Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame^] | 1 | # R A N D O M V A R I A B L E G E N E R A T O R S |
| 2 | # |
| 3 | # distributions on the real line: |
| 4 | # ------------------------------ |
| 5 | # normal (Gaussian) |
| 6 | # lognormal |
| 7 | # negative exponential |
| 8 | # gamma |
| 9 | # |
| 10 | # distributions on the circle (angles 0 to 2pi) |
| 11 | # --------------------------------------------- |
| 12 | # circular uniform |
| 13 | # von Mises |
| 14 | |
| 15 | # Translated from anonymously contributed C/C++ source. |
| 16 | |
| 17 | from whrandom import random, uniform, randint, choice # Also for export! |
| 18 | from math import log, exp, pi, e, sqrt, acos, cos |
| 19 | |
| 20 | # Housekeeping function to verify that magic constants have been |
| 21 | # computed correctly |
| 22 | |
| 23 | def verify(name, expected): |
| 24 | computed = eval(name) |
| 25 | if abs(computed - expected) > 1e-7: |
| 26 | raise ValueError, \ |
| 27 | 'computed value for %s deviates too much (computed %g, expected %g)' % \ |
| 28 | (name, computed, expected) |
| 29 | |
| 30 | # -------------------- normal distribution -------------------- |
| 31 | |
| 32 | NV_MAGICCONST = 4*exp(-0.5)/sqrt(2) |
| 33 | verify('NV_MAGICCONST', 1.71552776992141) |
| 34 | def normalvariate(mu, sigma): |
| 35 | # mu = mean, sigma = standard deviation |
| 36 | |
| 37 | # Uses Kinderman and Monahan method. Reference: Kinderman, |
| 38 | # A.J. and Monahan, J.F., "Computer generation of random |
| 39 | # variables using the ratio of uniform deviates", ACM Trans |
| 40 | # Math Software, 3, (1977), pp257-260. |
| 41 | |
| 42 | while 1: |
| 43 | u1 = random() |
| 44 | u2 = random() |
| 45 | z = NV_MAGICCONST*(u1-0.5)/u2 |
| 46 | zz = z*z/4 |
| 47 | if zz <= -log(u2): |
| 48 | break |
| 49 | return mu+z*sigma |
| 50 | |
| 51 | # -------------------- lognormal distribution -------------------- |
| 52 | |
| 53 | def lognormvariate(mu, sigma): |
| 54 | return exp(normalvariate(mu, sigma)) |
| 55 | |
| 56 | # -------------------- circular uniform -------------------- |
| 57 | |
| 58 | def cunifvariate(mean, arc): |
| 59 | # mean: mean angle (in radians between 0 and pi) |
| 60 | # arc: range of distribution (in radians between 0 and pi) |
| 61 | |
| 62 | return (mean + arc * (random() - 0.5)) % pi |
| 63 | |
| 64 | # -------------------- exponential distribution -------------------- |
| 65 | |
| 66 | def expovariate(lambd): |
| 67 | # lambd: rate lambd = 1/mean |
| 68 | # ('lambda' is a Python reserved word) |
| 69 | |
| 70 | u = random() |
| 71 | while u <= 1e-7: |
| 72 | u = random() |
| 73 | return -log(u)/lambd |
| 74 | |
| 75 | # -------------------- von Mises distribution -------------------- |
| 76 | |
| 77 | TWOPI = 2*pi |
| 78 | verify('TWOPI', 6.28318530718) |
| 79 | |
| 80 | def vonmisesvariate(mu, kappa): |
| 81 | # mu: mean angle (in radians between 0 and 180 degrees) |
| 82 | # kappa: concentration parameter kappa (>= 0) |
| 83 | |
| 84 | # if kappa = 0 generate uniform random angle |
| 85 | if kappa <= 1e-6: |
| 86 | return TWOPI * random() |
| 87 | |
| 88 | a = 1.0 + sqrt(1 + 4 * kappa * kappa) |
| 89 | b = (a - sqrt(2 * a))/(2 * kappa) |
| 90 | r = (1 + b * b)/(2 * b) |
| 91 | |
| 92 | while 1: |
| 93 | u1 = random() |
| 94 | |
| 95 | z = cos(pi * u1) |
| 96 | f = (1 + r * z)/(r + z) |
| 97 | c = kappa * (r - f) |
| 98 | |
| 99 | u2 = random() |
| 100 | |
| 101 | if not (u2 >= c * (2.0 - c) and u2 > c * exp(1.0 - c)): |
| 102 | break |
| 103 | |
| 104 | u3 = random() |
| 105 | if u3 > 0.5: |
| 106 | theta = mu + 0.5*acos(f) |
| 107 | else: |
| 108 | theta = mu - 0.5*acos(f) |
| 109 | |
| 110 | return theta % pi |
| 111 | |
| 112 | # -------------------- gamma distribution -------------------- |
| 113 | |
| 114 | LOG4 = log(4) |
| 115 | verify('LOG4', 1.38629436111989) |
| 116 | |
| 117 | def gammavariate(alpha, beta): |
| 118 | # beta times standard gamma |
| 119 | ainv = sqrt(2 * alpha - 1) |
| 120 | return beta * stdgamma(alpha, ainv, alpha - LOG4, alpha + ainv) |
| 121 | |
| 122 | SG_MAGICCONST = 1+log(4.5) |
| 123 | verify('SG_MAGICCONST', 2.50407739677627) |
| 124 | |
| 125 | def stdgamma(alpha, ainv, bbb, ccc): |
| 126 | # ainv = sqrt(2 * alpha - 1) |
| 127 | # bbb = alpha - log(4) |
| 128 | # ccc = alpha + ainv |
| 129 | |
| 130 | if alpha <= 0.0: |
| 131 | raise ValueError, 'stdgamma: alpha must be > 0.0' |
| 132 | |
| 133 | if alpha > 1.0: |
| 134 | |
| 135 | # Uses R.C.H. Cheng, "The generation of Gamma |
| 136 | # variables with non-integral shape parameters", |
| 137 | # Applied Statistics, (1977), 26, No. 1, p71-74 |
| 138 | |
| 139 | while 1: |
| 140 | u1 = random() |
| 141 | u2 = random() |
| 142 | v = log(u1/(1-u1))/ainv |
| 143 | x = alpha*exp(v) |
| 144 | z = u1*u1*u2 |
| 145 | r = bbb+ccc*v-x |
| 146 | if r + SG_MAGICCONST - 4.5*z >= 0 or r >= log(z): |
| 147 | return x |
| 148 | |
| 149 | elif alpha == 1.0: |
| 150 | # expovariate(1) |
| 151 | u = random() |
| 152 | while u <= 1e-7: |
| 153 | u = random() |
| 154 | return -log(u) |
| 155 | |
| 156 | else: # alpha is between 0 and 1 (exclusive) |
| 157 | |
| 158 | # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle |
| 159 | |
| 160 | while 1: |
| 161 | u = random() |
| 162 | b = (e + alpha)/e |
| 163 | p = b*u |
| 164 | if p <= 1.0: |
| 165 | x = pow(p, 1.0/alpha) |
| 166 | else: |
| 167 | # p > 1 |
| 168 | x = -log((b-p)/alpha) |
| 169 | u1 = random() |
| 170 | if not (((p <= 1.0) and (u1 > exp(-x))) or |
| 171 | ((p > 1) and (u1 > pow(x, alpha - 1.0)))): |
| 172 | break |
| 173 | return x |
| 174 | |
| 175 | # -------------------- test program -------------------- |
| 176 | |
| 177 | def test(): |
| 178 | print 'TWOPI =', TWOPI |
| 179 | print 'LOG4 =', LOG4 |
| 180 | print 'NV_MAGICCONST =', NV_MAGICCONST |
| 181 | print 'SG_MAGICCONST =', SG_MAGICCONST |
| 182 | N = 100 |
| 183 | test_generator(N, 'random()') |
| 184 | test_generator(N, 'normalvariate(0.0, 1.0)') |
| 185 | test_generator(N, 'lognormvariate(0.0, 1.0)') |
| 186 | test_generator(N, 'cunifvariate(0.0, 1.0)') |
| 187 | test_generator(N, 'expovariate(1.0)') |
| 188 | test_generator(N, 'vonmisesvariate(0.0, 1.0)') |
| 189 | test_generator(N, 'gammavariate(0.5, 1.0)') |
| 190 | test_generator(N, 'gammavariate(0.9, 1.0)') |
| 191 | test_generator(N, 'gammavariate(1.0, 1.0)') |
| 192 | test_generator(N, 'gammavariate(2.0, 1.0)') |
| 193 | test_generator(N, 'gammavariate(20.0, 1.0)') |
| 194 | test_generator(N, 'gammavariate(200.0, 1.0)') |
| 195 | |
| 196 | def test_generator(n, funccall): |
| 197 | import sys |
| 198 | print '%d calls to %s:' % (n, funccall), |
| 199 | sys.stdout.flush() |
| 200 | code = compile(funccall, funccall, 'eval') |
| 201 | sum = 0.0 |
| 202 | sqsum = 0.0 |
| 203 | for i in range(n): |
| 204 | x = eval(code) |
| 205 | sum = sum + x |
| 206 | sqsum = sqsum + x*x |
| 207 | avg = sum/n |
| 208 | stddev = sqrt(sqsum/n - avg*avg) |
| 209 | print 'avg %g, stddev %g' % (avg, stddev) |
| 210 | |
| 211 | if __name__ == '__main__': |
| 212 | test() |