Gamma distribution. Not the gamma function! Conditions on the parameters are alpha > 0 and beta > 0. The probability distribution function is: x ** (alpha - 1) * math.exp(-x / beta) pdf(x) = --------------------------------------
(self, alpha, beta)
| 665 | return theta |
| 666 | |
| 667 | def gammavariate(self, alpha, beta): |
| 668 | """Gamma distribution. Not the gamma function! |
| 669 | |
| 670 | Conditions on the parameters are alpha > 0 and beta > 0. |
| 671 | |
| 672 | The probability distribution function is: |
| 673 | |
| 674 | x ** (alpha - 1) * math.exp(-x / beta) |
| 675 | pdf(x) = -------------------------------------- |
| 676 | math.gamma(alpha) * beta ** alpha |
| 677 | |
| 678 | The mean (expected value) and variance of the random variable are: |
| 679 | |
| 680 | E[X] = alpha * beta |
| 681 | Var[X] = alpha * beta ** 2 |
| 682 | |
| 683 | """ |
| 684 | |
| 685 | # Warning: a few older sources define the gamma distribution in terms |
| 686 | # of alpha > -1.0 |
| 687 | if alpha <= 0.0 or beta <= 0.0: |
| 688 | raise ValueError('gammavariate: alpha and beta must be > 0.0') |
| 689 | |
| 690 | random = self.random |
| 691 | if alpha > 1.0: |
| 692 | |
| 693 | # Uses R.C.H. Cheng, "The generation of Gamma |
| 694 | # variables with non-integral shape parameters", |
| 695 | # Applied Statistics, (1977), 26, No. 1, p71-74 |
| 696 | |
| 697 | ainv = _sqrt(2.0 * alpha - 1.0) |
| 698 | bbb = alpha - LOG4 |
| 699 | ccc = alpha + ainv |
| 700 | |
| 701 | while True: |
| 702 | u1 = random() |
| 703 | if not 1e-7 < u1 < 0.9999999: |
| 704 | continue |
| 705 | u2 = 1.0 - random() |
| 706 | v = _log(u1 / (1.0 - u1)) / ainv |
| 707 | x = alpha * _exp(v) |
| 708 | z = u1 * u1 * u2 |
| 709 | r = bbb + ccc * v - x |
| 710 | if r + SG_MAGICCONST - 4.5 * z >= 0.0 or r >= _log(z): |
| 711 | return x * beta |
| 712 | |
| 713 | elif alpha == 1.0: |
| 714 | # expovariate(1/beta) |
| 715 | return -_log(1.0 - random()) * beta |
| 716 | |
| 717 | else: |
| 718 | # alpha is between 0 and 1 (exclusive) |
| 719 | # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle |
| 720 | while True: |
| 721 | u = random() |
| 722 | b = (_e + alpha) / _e |
| 723 | p = b * u |
| 724 | if p <= 1.0: |
no outgoing calls