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 |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 9 | # beta |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 10 | # |
| 11 | # distributions on the circle (angles 0 to 2pi) |
| 12 | # --------------------------------------------- |
| 13 | # circular uniform |
| 14 | # von Mises |
| 15 | |
| 16 | # Translated from anonymously contributed C/C++ source. |
| 17 | |
| 18 | from whrandom import random, uniform, randint, choice # Also for export! |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 19 | from math import log, exp, pi, e, sqrt, acos, cos, sin |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 20 | |
| 21 | # Housekeeping function to verify that magic constants have been |
| 22 | # computed correctly |
| 23 | |
| 24 | def verify(name, expected): |
| 25 | computed = eval(name) |
| 26 | if abs(computed - expected) > 1e-7: |
| 27 | raise ValueError, \ |
| 28 | 'computed value for %s deviates too much (computed %g, expected %g)' % \ |
| 29 | (name, computed, expected) |
| 30 | |
| 31 | # -------------------- normal distribution -------------------- |
| 32 | |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 33 | NV_MAGICCONST = 4*exp(-0.5)/sqrt(2.0) |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 34 | verify('NV_MAGICCONST', 1.71552776992141) |
| 35 | def normalvariate(mu, sigma): |
| 36 | # mu = mean, sigma = standard deviation |
| 37 | |
| 38 | # Uses Kinderman and Monahan method. Reference: Kinderman, |
| 39 | # A.J. and Monahan, J.F., "Computer generation of random |
| 40 | # variables using the ratio of uniform deviates", ACM Trans |
| 41 | # Math Software, 3, (1977), pp257-260. |
| 42 | |
| 43 | while 1: |
| 44 | u1 = random() |
| 45 | u2 = random() |
| 46 | z = NV_MAGICCONST*(u1-0.5)/u2 |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 47 | zz = z*z/4.0 |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 48 | if zz <= -log(u2): |
| 49 | break |
| 50 | return mu+z*sigma |
| 51 | |
| 52 | # -------------------- lognormal distribution -------------------- |
| 53 | |
| 54 | def lognormvariate(mu, sigma): |
| 55 | return exp(normalvariate(mu, sigma)) |
| 56 | |
| 57 | # -------------------- circular uniform -------------------- |
| 58 | |
| 59 | def cunifvariate(mean, arc): |
| 60 | # mean: mean angle (in radians between 0 and pi) |
| 61 | # arc: range of distribution (in radians between 0 and pi) |
| 62 | |
| 63 | return (mean + arc * (random() - 0.5)) % pi |
| 64 | |
| 65 | # -------------------- exponential distribution -------------------- |
| 66 | |
| 67 | def expovariate(lambd): |
| 68 | # lambd: rate lambd = 1/mean |
| 69 | # ('lambda' is a Python reserved word) |
| 70 | |
| 71 | u = random() |
| 72 | while u <= 1e-7: |
| 73 | u = random() |
| 74 | return -log(u)/lambd |
| 75 | |
| 76 | # -------------------- von Mises distribution -------------------- |
| 77 | |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 78 | TWOPI = 2.0*pi |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 79 | verify('TWOPI', 6.28318530718) |
| 80 | |
| 81 | def vonmisesvariate(mu, kappa): |
Guido van Rossum | 5810297 | 1998-04-06 14:12:13 +0000 | [diff] [blame] | 82 | # mu: mean angle (in radians between 0 and 2*pi) |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 83 | # kappa: concentration parameter kappa (>= 0) |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 84 | # if kappa = 0 generate uniform random angle |
Guido van Rossum | 5810297 | 1998-04-06 14:12:13 +0000 | [diff] [blame] | 85 | |
| 86 | # Based upon an algorithm published in: Fisher, N.I., |
| 87 | # "Statistical Analysis of Circular Data", Cambridge |
| 88 | # University Press, 1993. |
| 89 | |
| 90 | # Thanks to Magnus Kessler for a correction to the |
| 91 | # implementation of step 4. |
| 92 | |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 93 | if kappa <= 1e-6: |
| 94 | return TWOPI * random() |
| 95 | |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 96 | a = 1.0 + sqrt(1.0 + 4.0 * kappa * kappa) |
| 97 | b = (a - sqrt(2.0 * a))/(2.0 * kappa) |
| 98 | r = (1.0 + b * b)/(2.0 * b) |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 99 | |
| 100 | while 1: |
| 101 | u1 = random() |
| 102 | |
| 103 | z = cos(pi * u1) |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 104 | f = (1.0 + r * z)/(r + z) |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 105 | c = kappa * (r - f) |
| 106 | |
| 107 | u2 = random() |
| 108 | |
| 109 | if not (u2 >= c * (2.0 - c) and u2 > c * exp(1.0 - c)): |
| 110 | break |
| 111 | |
| 112 | u3 = random() |
| 113 | if u3 > 0.5: |
Guido van Rossum | 5810297 | 1998-04-06 14:12:13 +0000 | [diff] [blame] | 114 | theta = (mu % TWOPI) + acos(f) |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 115 | else: |
Guido van Rossum | 5810297 | 1998-04-06 14:12:13 +0000 | [diff] [blame] | 116 | theta = (mu % TWOPI) - acos(f) |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 117 | |
Guido van Rossum | 5810297 | 1998-04-06 14:12:13 +0000 | [diff] [blame] | 118 | return theta |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 119 | |
| 120 | # -------------------- gamma distribution -------------------- |
| 121 | |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 122 | LOG4 = log(4.0) |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 123 | verify('LOG4', 1.38629436111989) |
| 124 | |
| 125 | def gammavariate(alpha, beta): |
| 126 | # beta times standard gamma |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 127 | ainv = sqrt(2.0 * alpha - 1.0) |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 128 | return beta * stdgamma(alpha, ainv, alpha - LOG4, alpha + ainv) |
| 129 | |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 130 | SG_MAGICCONST = 1.0 + log(4.5) |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 131 | verify('SG_MAGICCONST', 2.50407739677627) |
| 132 | |
| 133 | def stdgamma(alpha, ainv, bbb, ccc): |
| 134 | # ainv = sqrt(2 * alpha - 1) |
| 135 | # bbb = alpha - log(4) |
| 136 | # ccc = alpha + ainv |
| 137 | |
| 138 | if alpha <= 0.0: |
| 139 | raise ValueError, 'stdgamma: alpha must be > 0.0' |
| 140 | |
| 141 | if alpha > 1.0: |
| 142 | |
| 143 | # Uses R.C.H. Cheng, "The generation of Gamma |
| 144 | # variables with non-integral shape parameters", |
| 145 | # Applied Statistics, (1977), 26, No. 1, p71-74 |
| 146 | |
| 147 | while 1: |
| 148 | u1 = random() |
| 149 | u2 = random() |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 150 | v = log(u1/(1.0-u1))/ainv |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 151 | x = alpha*exp(v) |
| 152 | z = u1*u1*u2 |
| 153 | r = bbb+ccc*v-x |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 154 | if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= log(z): |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 155 | return x |
| 156 | |
| 157 | elif alpha == 1.0: |
| 158 | # expovariate(1) |
| 159 | u = random() |
| 160 | while u <= 1e-7: |
| 161 | u = random() |
| 162 | return -log(u) |
| 163 | |
| 164 | else: # alpha is between 0 and 1 (exclusive) |
| 165 | |
| 166 | # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle |
| 167 | |
| 168 | while 1: |
| 169 | u = random() |
| 170 | b = (e + alpha)/e |
| 171 | p = b*u |
| 172 | if p <= 1.0: |
| 173 | x = pow(p, 1.0/alpha) |
| 174 | else: |
| 175 | # p > 1 |
| 176 | x = -log((b-p)/alpha) |
| 177 | u1 = random() |
| 178 | if not (((p <= 1.0) and (u1 > exp(-x))) or |
| 179 | ((p > 1) and (u1 > pow(x, alpha - 1.0)))): |
| 180 | break |
| 181 | return x |
| 182 | |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 183 | |
| 184 | # -------------------- Gauss (faster alternative) -------------------- |
| 185 | |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 186 | gauss_next = None |
| 187 | def gauss(mu, sigma): |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 188 | |
| 189 | # When x and y are two variables from [0, 1), uniformly |
| 190 | # distributed, then |
| 191 | # |
Guido van Rossum | 72c2e1b | 1998-02-19 21:17:42 +0000 | [diff] [blame] | 192 | # cos(2*pi*x)*sqrt(-2*log(1-y)) |
| 193 | # sin(2*pi*x)*sqrt(-2*log(1-y)) |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 194 | # |
| 195 | # are two *independent* variables with normal distribution |
| 196 | # (mu = 0, sigma = 1). |
| 197 | # (Lambert Meertens) |
Guido van Rossum | 72c2e1b | 1998-02-19 21:17:42 +0000 | [diff] [blame] | 198 | # (corrected version; bug discovered by Mike Miller, fixed by LM) |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 199 | |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 200 | global gauss_next |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 201 | |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 202 | if gauss_next != None: |
| 203 | z = gauss_next |
| 204 | gauss_next = None |
| 205 | else: |
| 206 | x2pi = random() * TWOPI |
Guido van Rossum | 72c2e1b | 1998-02-19 21:17:42 +0000 | [diff] [blame] | 207 | g2rad = sqrt(-2.0 * log(1.0 - random())) |
| 208 | z = cos(x2pi) * g2rad |
| 209 | gauss_next = sin(x2pi) * g2rad |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 210 | |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 211 | return mu + z*sigma |
| 212 | |
| 213 | # -------------------- beta -------------------- |
| 214 | |
| 215 | def betavariate(alpha, beta): |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 216 | |
| 217 | # Discrete Event Simulation in C, pp 87-88. |
| 218 | |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 219 | y = expovariate(alpha) |
| 220 | z = expovariate(1.0/beta) |
| 221 | return z/(y+z) |
| 222 | |
Guido van Rossum | 5bdea89 | 1997-12-09 19:43:18 +0000 | [diff] [blame] | 223 | # -------------------- Pareto -------------------- |
Guido van Rossum | cf4559a | 1997-12-02 02:47:39 +0000 | [diff] [blame] | 224 | |
| 225 | def paretovariate(alpha): |
| 226 | # Jain, pg. 495 |
| 227 | |
| 228 | u = random() |
| 229 | return 1.0 / pow(u, 1.0/alpha) |
| 230 | |
Guido van Rossum | 5bdea89 | 1997-12-09 19:43:18 +0000 | [diff] [blame] | 231 | # -------------------- Weibull -------------------- |
Guido van Rossum | cf4559a | 1997-12-02 02:47:39 +0000 | [diff] [blame] | 232 | |
| 233 | def weibullvariate(alpha, beta): |
| 234 | # Jain, pg. 499; bug fix courtesy Bill Arms |
| 235 | |
| 236 | u = random() |
| 237 | return alpha * pow(-log(u), 1.0/beta) |
| 238 | |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 239 | # -------------------- test program -------------------- |
| 240 | |
Guido van Rossum | 2922c6d | 1994-05-06 14:28:19 +0000 | [diff] [blame] | 241 | def test(N = 200): |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 242 | print 'TWOPI =', TWOPI |
| 243 | print 'LOG4 =', LOG4 |
| 244 | print 'NV_MAGICCONST =', NV_MAGICCONST |
| 245 | print 'SG_MAGICCONST =', SG_MAGICCONST |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 246 | test_generator(N, 'random()') |
| 247 | test_generator(N, 'normalvariate(0.0, 1.0)') |
| 248 | test_generator(N, 'lognormvariate(0.0, 1.0)') |
| 249 | test_generator(N, 'cunifvariate(0.0, 1.0)') |
| 250 | test_generator(N, 'expovariate(1.0)') |
| 251 | test_generator(N, 'vonmisesvariate(0.0, 1.0)') |
| 252 | test_generator(N, 'gammavariate(0.5, 1.0)') |
| 253 | test_generator(N, 'gammavariate(0.9, 1.0)') |
| 254 | test_generator(N, 'gammavariate(1.0, 1.0)') |
| 255 | test_generator(N, 'gammavariate(2.0, 1.0)') |
| 256 | test_generator(N, 'gammavariate(20.0, 1.0)') |
| 257 | test_generator(N, 'gammavariate(200.0, 1.0)') |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 258 | test_generator(N, 'gauss(0.0, 1.0)') |
| 259 | test_generator(N, 'betavariate(3.0, 3.0)') |
Guido van Rossum | cf4559a | 1997-12-02 02:47:39 +0000 | [diff] [blame] | 260 | test_generator(N, 'paretovariate(1.0)') |
| 261 | test_generator(N, 'weibullvariate(1.0, 1.0)') |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 262 | |
| 263 | def test_generator(n, funccall): |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 264 | import time |
| 265 | print n, 'times', funccall |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 266 | code = compile(funccall, funccall, 'eval') |
| 267 | sum = 0.0 |
| 268 | sqsum = 0.0 |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 269 | smallest = 1e10 |
Guido van Rossum | cc32ac9 | 1994-03-15 16:10:24 +0000 | [diff] [blame] | 270 | largest = -1e10 |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 271 | t0 = time.time() |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 272 | for i in range(n): |
| 273 | x = eval(code) |
| 274 | sum = sum + x |
| 275 | sqsum = sqsum + x*x |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 276 | smallest = min(x, smallest) |
| 277 | largest = max(x, largest) |
| 278 | t1 = time.time() |
| 279 | print round(t1-t0, 3), 'sec,', |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 280 | avg = sum/n |
| 281 | stddev = sqrt(sqsum/n - avg*avg) |
Guido van Rossum | 95bfcda | 1994-03-09 14:21:05 +0000 | [diff] [blame] | 282 | print 'avg %g, stddev %g, min %g, max %g' % \ |
| 283 | (avg, stddev, smallest, largest) |
Guido van Rossum | ff03b1a | 1994-03-09 12:55:02 +0000 | [diff] [blame] | 284 | |
| 285 | if __name__ == '__main__': |
| 286 | test() |