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