1
from __future__ import nested_scopes
3
import string, time, array
12
special = scipy.special
13
from scipy_base.fastumath import *
14
from scipy_base import vectorize
17
SequenceType = [types.TupleType, types.ListType, array.ArrayType, Num.ArrayType]
20
if type(sh) not in SequenceType:
23
if not isinstance(val,types.IntType):
24
raise ValueError, "Each element of the shape parameter must be an integer."
25
prod = Num.product(sh)
26
return tuple(sh), prod
29
ArgumentError = "ArgumentError"
31
def multivariate_normal(mean, cov, size=None):
32
"""returns an array containing multivariate normally distributed random
33
numbers with specified mean and covariance.
35
mean must be a 1 dimensional array. cov must be a square two dimensional
36
array with the same number of rows and columns as mean has elements.
38
The first form returns a single 1-D array containing a multivariate
41
The second form returns an array of shape (m, n, ..., cov.shape[0]).
42
In this case, output[i,j,...,:] is a 1-D array containing a multivariate
44
if isinstance(size, IntType):
50
output = rand.multivariate_normal(mean, cov, n)
52
final_shape = list(size[:])
53
final_shape.append(mean.shape[0])
54
output.shape = final_shape
57
#####################################
58
# General purpose continuous
59
######################################
61
def randwppf(ppf, args=(), size=None):
62
"""returns an array of randomly distributed integers of a distribution
63
whose percent point function (inverse of the CDF) is given.
65
args is a tuple of extra arguments to the ppf function (i.e. shape,
66
location, scale), and size is the size of the output. Note the ppf
67
function must accept an array of q values to compute over.
70
return apply(ppf, (U,)+args)
72
def randwcdf(cdf, mean=1.0, args=(), size=None):
73
"""returns an array of randomly distributed integers of a distribution
74
whose cumulative distribution function (CDF) is given.
76
mean is the mean of the distribution (helps the solver).
77
args is a tuple of extra arguments to the cdf function (i.e. shape,
78
location, scale), and size is the size of the output. Note the
79
cdf function needs to accept a single value to compute over.
81
import scipy.optimize as optimize
82
def _ppfopt(x, q, *nargs):
84
return cdf(*newargs) - q
87
return optimize.fsolve(_ppfopt, mean, args=(q,)+nargs)
89
_vppf = vectorize(_ppf)
91
return apply(_vppf,(U,)+args)
94
#################################################
96
##################################################
98
def multinom(trials, probs, size=None):
99
"""returns array of multinomial distributed integer vectors.
101
trials is the number of trials in each multinomial distribution.
102
probs is a one dimensional array. There are len(prob)+1 events.
103
prob[i] is the probability of the i-th event, 0<=i<len(prob).
104
The probability of event len(prob) is 1.-Numeric.sum(prob).
106
The first form returns a single 1-D array containing one multinomially
109
The second form returns an array of size (m, n, ..., len(probs)).
110
In this case, output[i,j,...,:] is a 1-D array containing a multinomially
111
distributed integer 1-D array."""
112
# Check preconditions on arguments
113
probs = Num.array(probs)
114
if len(probs.shape) != 1:
115
raise ArgumentError, "probs must be 1 dimensional."
116
# Compute shape of output
117
if type(size) == type(0): size = [size]
118
final_shape = size[:]
119
final_shape.append(probs.shape[0]+1)
120
x = rand.multinomial(trials, probs.astype(Num.Float32), Num.multiply.reduce(size))
121
# Change its shape to the desire one
122
x.shape = final_shape
126
##############################
127
## Functions from old rv2.py
128
##############################
132
"""Encapsulating class for all random number generator methods and data
134
All <_pranv> data and method attributes are private. The single
135
instance <_inst> is created and public aliases to external methods are
136
provided at the end of the module. To get the alias, delete the leading
137
underscore from the name of the corresponding method in class <_pranv>.
138
There are no public data attributes. This class should have only one
139
instance, <_inst>."""
141
# Internal (private) data attributes:
144
# _algorithm_name one-line description of uniform(0,1) algorithm used
145
# _count number of delivered random()s; Python long integer
146
# _index integer index to next uniform(0,1) in ranbuf[]
147
# _iterator a tuple of tuples used by some uniform generators
148
# _ln_fac a double array containing ln(n!) for n = 0, ..., 129
149
# _ranbuf a double array to store batched uniform(0,1) randoms
150
# _ranbuf_size integer, length of _ranbuf
151
# _fillbuf a function variable, aliasing a uniform(0,1) generator
152
# _seed initial Python long integer from user or clock
153
# _series array of pseudo-random series values, or single value
154
# _second_nrv Normal pseudo-randoms are produced in pairs; 2nd here.
156
# Internal (private) utility methods:
158
# __init__() Initialize <_pranv> instance _inst; called once.
159
# _build_iterator() Construct _iterator; allocate _series, _ranbuf.
160
# _cmrg(size=None) Fill ranbuf with random #s from algorithm 'cmrg'.
161
# _fill_ln_fac() Calculate and store ln(n!) for n = 0 through 129.
162
# _flip(size=None) Fill ranbuf with randoms from algorithm 'flip'.
163
# _ln_factorial(n) Return ln(n!), from ln_fac[n] or by calculation.
164
# _smrg(size=None) Fill ranbuf with random #s from algorithm 'sMRG'.
165
# _twister(size=None) Fill ranbuf with randoms from algorithm 'twister'.
167
# External (private, but aliased public) utility methods:
169
# _initial_seed() Return the seed used to start this random series.
170
# _initialize(file_string=None, algorithm='cmrg', seed=0L) Restart.
171
# _random_algorithm() Return a 1-line description of current generator.
172
# _random_count() Return number of randoms generated in current series.
173
# _random(size=None) Return next pseudo-random uniform(0,1).
174
# _save_state(file_string='prvstate.dat') Save <_inst> state to disk.
176
# External (private, but aliased public) non-uniform random variate methods:
178
# _choice(seq=(0,1), size=None) <seq> item
179
# _geom(pr_failure=0.5,
180
# size=None) integer non-negative
181
# _hypergeom(tot=35, good=25,
182
# sample=10, size=None) integer >= max(0, sample-good)
183
# <= min(sample, bad)
184
# _logser(p=0.5, size=None) integer positive
185
# _von_Mises(mode=0.0, shape=1.0,
186
# size=None) double > -pi, < +pi
187
# _Wald(mean=1.0, scale=1.0,
188
# size=None) double positive
189
# _Zipf(a=4.0, size=None) integer positive
191
# External (private, but aliased public) geometrical,
192
# permutation, and subsampling routines:
194
# _in_simplex(mseq=5*[0.0], bound=1.0) Return point in simplex.
195
# _in_sphere(center=5*[0.0], radius=1.0) Return point in sphere.
196
# _on_simplex(mseq=5*[0.0], bound=1.0) Return point on simplex.
197
# _on_sphere(center=5*[0.0], radius=1.0) Return point on sphere.
198
# _sample(inlist=[0,1,2], sample_size=2) Return simple random sample.
199
# _smart_sample(inlist=[0,1,2], sample_size=2) Return 'no-dups' sample.
203
"""Initialize the single class <_pranv> instance.
206
# Declaration of class _pranv data attributes, all introduced here mostly for
207
# documentation purposes. Most are allocated or set in _initialize().
209
self._algorithm_name = '' # one-line aLgorithm description
210
self._count = 0L # number of pseudorandoms generated in this series
211
self._index = 0 # index to output buffer array ranbuf[]
212
self._iterator = ((),) # tuple of tuples used by 'flip' and 'twister'.
213
self._ln_fac = array.array('d', 130*[0.0]) # ln(n!); n = 0 to 129;
214
# double array used by some discrete generators.
215
self._ranbuf = [] # uniform(0,1) output buffer array; allocated in
216
# _build_iterator(), called by _initialize().
217
self._ranbuf_size = 0 # This is just len(_ranbuf), after allocation.
218
self._fillbuf = None # function variable pointing to uniform(0,1) gen.
219
self._seed = 0L # starting seed from user or system clock
220
self._series = [] # rotating array of random integers, longs or
221
# doubles; allocated in _build_iterator.
222
self._second_nrv = None # 2nd of generated pair of normal random vars.
223
# End of class <_pranv> data attribute declarations
225
self._fill_ln_fac() # Calculate values of _ln_fac[n] = ln(n!)
226
# for n = 0 to 129; only need to do this once.
227
self._initialize() # Do detailed initialization. <_initialize()>
228
# may also be called by the user (via its global
229
# alias that drops the leading '_') to change
232
def _build_iterator(self):
233
"""Construct <_iterator[]> if necessary; allocate _ranbuf and _series.
237
'flip' and 'twister', pseudo-random integer generators which manipulate
238
long series of integers, utilize a tuple of tuples to simplify and speed
239
up the logic. This function constructs the (global) tuple of tuples,
240
<_iterator>, if it's not already there. The array <_ranbuf> is
241
allocated if needed and its length, <_ranbuf_size>, is set. Also, the
242
array <_series> is allocated if necessary. <_build_iterator> is an
243
internal (private) <_pranv> method, called whenever the selected
244
uniform(0,1) pseudo-random generator changes, by <_initialize()>."""
246
if self._fillbuf == self._cmrg:
247
if self._ranbuf_size != 660: # 'cmrg' does not use _iterator[].
249
self._ranbuf = array.array('d', 660*[0.0]) # output buffer
250
self._ranbuf_size = 660
251
if len(self._series) != 6:
253
self._series = array.array('d', 6*[0.0])
255
elif self._fillbuf == self._flip:
256
if len(self._iterator) != 55: # May need to recreate _ranbuf ,
257
del self._series # _iterator and _series.
258
self._series = array.array('l',56*[0])
260
self._iterator = 55 * [()]
261
i = 1 # The 'flip' _iterator is:
262
j = 32 # ( (1,32), (2,33), ..., (24,55), (25,1),
263
while j < 56: # (26,2), ..., (55,31) ). Note that
264
self._iterator[i-1] = (i, j) # _iterator is global. A list
265
i = i + 1 # here, _iterator converts to tuple below.
270
self._iterator[i-1] = (i, j)
274
self._iterator = tuple(self._iterator)
276
if self._ranbuf_size != 660:
278
self._ranbuf = array.array('d', 660*[0.0]) # output buffer
279
self._ranbuf_size = 660
281
elif self._fillbuf == self._smrg:
282
if self._ranbuf_size != 660: # 'SMRG' does not use _iterator[].
284
self._ranbuf = array.array('d', 660*[0.0]) # output buffer
285
self._ranbuf_size = 660
288
self._series = 0L # _series is a single Python long.
290
elif self._fillbuf == self._twister:
291
if self._ranbuf_size != 624: # May need to reload _ranbuf,
292
del self._ranbuf # _series and _iterator.
293
self._ranbuf = array.array('d', 624*[0.0])
294
self._ranbuf_size = 624
296
self._series = array.array('L', 624*[0L])
298
self._iterator = 624 * [()]
300
k = 0 # ( (0,1,397),(1,2,398),...,(226,227,623),
301
while k < 227: # (227,228,0),228,229,1),...,(622,623,395),
302
self._iterator[k] = (k, k + 1, k + 397) # (623,0,386) ).A list
303
k = k + 1 # here, _iterator converts to tuple below.
306
self._iterator[k] = (k, k + 1, k - 227)
309
self._iterator[623] = (623, 0, 396)
310
self._iterator = tuple(self._iterator)
313
pass # Can't get here.
316
"""Produce batch of uniform(0,1) RVs from L'Ecuyer's MRG32k3a generator
320
See L'Ecuyer, Pierre., "Good parameters and implementations for combined
321
multiple recursive random number generators," May, 1998, to appear in
322
"Operations Research." To download a Postscript copy of the article:
323
http://www.iro.umontreal.ca/~lecuyer/myftp/papers/combmrg2.ps
324
The period is about 2**191, (or 10**57), and it is well-behaved in all
325
dimensions through 45. The generator has been tested extensively by
326
L'Ecuyer; no statistical faults were found."""
328
# This algorithm is much simpler than the implementation appearing below,
329
# which has been hand-optimized to minimize array references. The savings
330
# in execution time was about 25%. Here is the underlying algorithm:
332
# s = self.series # local alias
334
# p1 = 1403580.0 * s[1] - 810728.0 * s[0]
335
# k = floor(p1 / 4294967087.0) # first modulus, 2**32 - 209
336
# p1 = p1 - k * 4294967087.0
337
# if p1 < 0.0: p1 = p1 + 4294967087.0
339
# p2 = 527612.0 * s[5] - 1370589.0 * s[3]
340
# k = floor(p2 / 4294944443.0) # 2nd modulus, 2**32 - 22853
341
# p2 = p2 - k * 4294944443.0
342
# if p2 < 0.0: p2 = p2 + 4294944443.0
347
# if dif <= 0.0: return (dif + 4294967087 * 2.328306549295728e-10
348
# else: return dif * 2.328306549295728e-10
351
ranlim = self._ranbuf_size # local alias
352
buf = self._ranbuf # local alias
353
s = self._series # local alias
354
(s0, s1, s2, s3, s4, s5) = s # Unpack to local aliases.
357
p11 = 1403580.0 * s1 - 810728.0 * s0
358
p11 = p11 - floor(p11 / 4294967087.0) * 4294967087.0
359
if p11 < 0.0: p11 = p11 + 4294967087.0 # first mod, 2**32-209
361
p2 = 527612.0 * s5 - 1370589.0 * s3
362
p2 = p2 - floor(p2 / 4294944443.0) * 4294944443.0
363
if p2 < 0.0: p2 = p2 + 4294944443.0 # 2nd mod, 2**32-22853
368
if dif < 0: dif = dif + 4294967087.0
369
buf[i] = dif * 2.328306549295728e-10
372
p1 = 1403580.0 * s2 - 810728.0 * s1
373
p1 = p1 - floor(p1 / 4294967087.0) * 4294967087.0
374
if p1 < 0.0: p1 = p1 + 4294967087.0 # first mod, 2**32-209
376
p2 = 527612.0 * p2 - 1370589.0 * s4
377
p2 = p2 - floor(p2 / 4294944443.0) * 4294944443.0
378
if p2 < 0.0: p2 = p2 + 4294944443.0 # 2nd mod, 2**32-22853
383
if dif < 0: dif = dif + 4294967087.0
384
buf[i] = dif * 2.328306549295728e-10
387
p1 = 1403580.0 * p11 - 810728.0 * s2
388
p1 = p1 - floor(p1 / 4294967087.0) * 4294967087.0
389
if p1 < 0.0: p1 = p1 + 4294967087.0 # first mod, 2**32-209
391
p2 = 527612.0 * p2 - 1370589.0 * s5
392
p2 = p2 - floor(p2 / 4294944443.0) * 4294944443.0
393
if p2 < 0.0: p2 = p2 + 4294944443.0 # 2nd mod, 2**32-22853
398
if dif < 0: dif = dif + 4294967087.0
399
buf[i] = dif * 2.328306549295728e-10
402
s[0] = s0 # Save the coefficients
403
s[1] = s1 # for the next batch of
404
s[2] = s2 # 660 uniform randoms.
405
s[3] = s3 # s = [s0,s1,s2,s3,s4,s5]
406
s[4] = s4 # doesn't work here because
407
s[5] = s5 # s is an array, not a list.
409
self._count = self._count + ranlim
412
def _fill_ln_fac(self):
413
"""Calculate and store array values of <_ln_fac[n]> = ln(n!).
417
This internal routine is called only once, by __init__() in
418
class <_pranv> during the instantiation of <_inst>."""
419
self._ln_fac[0] = 0.0
424
sum = sum + log( float(i) )
425
self._ln_fac[i] = sum
429
"""Produce a batch of uniform(0,1) RVs with Knuth's 1993 generator.
433
See Knuth, D.E., "The Stanford GraphBase: A Platform for Combinatorial
434
Computing," ACM Press, Addison-Wesley, 1993, pp 216-221. The period for
435
the underlying pseudo-random integer series is at least 3.6e16, and may
436
be nearly 4e25. The low-order bits of the integers are just as random
437
as the high-order bits. The requirements are 32-bit integers and two's
438
complement arithmetic. The recurrence is quite fast, but for sensitive
439
applications, autocorrelations or other statistical defects may dictate
440
the use of a more complex generator, a longer series or both."""
441
ranlim = self._ranbuf_size # local alias
442
buf = self._ranbuf # local alias
443
s = self._series # local alias
444
# Here is the algorithm, C-language style:
445
# ii = 1 # Refresh series of 55 random integers.
447
# while jj < 56: # Calculate (s[ii] - s[ii+31]) % 2**31.
448
# s[ii] = (s[ii] - s[jj]) & 0x7fffffff
452
# while jj < 32: # Calculate (s[ii] - s[ii-24]) % 2**31.
453
# s[ii] = (s[ii] - s[jj]) & 0x7fffffff
457
# And here it is using <_iterator>:
461
for (ii, jj) in self._iterator:
462
s[ii] = (s[ii] - s[jj]) & 0x7fffffff
463
buf[i] = ( float(s[ii]) + 1.0 ) * 4.6566128709089882e-10
464
# The constant is 1/(2**31+1), so stored value is in (0,1).
468
self._count = self._count + ranlim
471
def _ln_factorial(self, n=0.0):
472
"""Return ln(n!)from the array <_ln_fac[n]> or use Stirling's formula.
476
For 0 <= n <= 129, the value pre-calculated and stored in <_ln_fac[]>
477
by internal routine <_fill_ln_fac()> is returned. For n >= 130, the
478
Stirling asymptotic approximation with 3 terms is used.
479
Approximation error is negligible. NOTE: n must be non-negative
480
and should be integral. There is no error-checking in this function
481
because it is internal, and we don't make mistakes inside this class."""
484
return self._ln_fac[int(n)]
489
return (((( 0.793650793650793650e-3) * yi2 # 1 / (1260 * (n+1)**5)
490
- 0.277777777777777778e-2) * yi2 # -1 / ( 360 * (n+1)**3)
491
+ 0.833333333333333333e-1) * yi # +1 / ( 8 * (n+1) )
492
+ 0.918938533204672742 # +(log(sqrt(2*math.pi))
493
+ (y - 0.5) * log(y) - y )
496
"""Produce a batch of uniform(0,1) RVs from Wu's mod 2**61 - 1 generator
500
See Wu, Pei-Chi, "Multiplicative, Congruential Random-Number Generators
501
with multiplier (+ or -) 2**k1 (+ or -) 2**k2 and Modulus 2**p - 1, ACM
502
Transactions on mathematical Software, June, 1997, Vol. 23, No. 2,
503
pp 255 - 265. The generator has modulus 2**61 - 1, which is the
504
Mersenne prime immediately following 2**31 - 1, and multiplier
505
37**458191 % (2**61 - 1). Because 37 is the minimal primitive root of
506
2**61 -1, the generator has full period (2**61 -2, about 2.3e18).
507
It was found by Wu to have the best performance on spectral tests
508
(through dimension 8) of any multiplier of the type 37**k % (2**61 - 1)
509
with k in the range from 1 to 1,000,000.
510
Generated pseudo-random numbers range from 4.337e-19 to 1.0 - 2**53
511
-- all bits are random, an advantage over 31- or 32-bit generators. """
512
buf = self._ranbuf # local alias
513
s = self._series # local alias
514
ranlim = self._ranbuf_size # local alias
515
for i in xrange(ranlim):
516
s = s * 2137866620694229420L % 0x1fffffffffffffffL # The first
517
# constant is 37**458191 % (2**61 - 1),
518
# and the second is 2**61-1.
519
buf[i] = float(s) * 4.3368086899420168e-19 # For s = 2**61 - 2,
520
# this is 1 - 2**-53.
522
self._series = s # save the current seed in <_series>
523
self._count = self._count + ranlim
526
"""Produce a batch of uniform(0,1) RVs from the Mersenne Twister MT19937
530
See M. matsumoto and T. Nishamura, "Mersenne Twister," ACM Transactions
531
on Modeling and Computer Simulation, Jan. 1998, vol. 8, no. 1, pp 3-30.
532
The period is 2**19937 - 1,(> 1e6000); the generator has a 623 dim-
533
ensional equi-distributional property. This means that every sequence
534
of less than 624 pseudo-random 32-bit integers occurs equally often. It
535
has passed a number of statistical tests, including those in Marsaglia's
536
"Diehard" package. """
537
buf = self._ranbuf # local alias
538
s = self._series # local alias
539
# Here is the first part of the algorithm, C-language style:
540
# if i == 624: # self._index points to long ints in _series[].
541
# k = 0 # Generate 624 new pseudo-random long integers.
542
# while k < 227: # 227 = 624 - 397
543
# y = (s[ k ] & 0x80000000L) | (s[k + 1] & 0x7fffffffL)
544
# s[k] = s[k + 397] ^ (y >> 1) ^ (y & 0x1) * 0x9908b0dfL
548
# y = (s[ k ] & 0x80000000L) | (s[k + 1] & 0x7fffffffL)
549
# s[k] = s[k - 227] ^ (y >> 1) ^ (y & 0x1) * 0x9908b0dfL
552
# y = (s[623] & 0x80000000L) | (s[ 0 ] & 0x7fffffffL)
553
# s[623] = s[ 396 ] ^ (y >> 1) ^ (y & 0x1) * 0x9908b0dfL
555
# And here is is in Python, using the global iterator, just three lines:
556
for (k, kp1, koff) in self._iterator:
557
y = (s[k] & 0x80000000L) | (s[kp1] & 0x7fffffffL)
558
y = s[koff] ^ (y >>1) ^ (y & 0x1) * 0x9908b0dfL
561
y = y ^ (y << 7) & 0x9d2c5680L
562
y = y ^ (y << 15) & 0xefc60000L
564
buf[k] = (float(y) + 1.0) * 2.3283064359965952e-10
565
# The constant is 1.0 / (2.0**32 +1.0).
567
self._count = self._count + self._ranbuf_size
569
# External (private; aliased public) utility methods:
571
def _initialize(self, file_string=None, algorithm='cmrg', seed=0L):
572
"""Set uniform algorithm and seed, or restore state from disk.
574
Called by <__init__()> or by user, via the alias <initialize()>.
575
if <file_string> is specified, it can be a null string, indicating
576
that class _pranv state is not to be retrieved from a file; a single
577
blank or the string 'default', indicating that _pranv state is to be
578
restored from the file with the default file name ('prvstate.dat'); or
579
the path/file name of a file containing a _pranv state. <algorithm> may
580
be 'cmrg', 'flip', 'smrg', or 'twister'; see the documentation for
581
guidance. <seed> is by default 0, indicating the system clock is to
582
be used to determine algorithm starting value(s), or a Python long
583
integer. For some of the algorithms, a negative seed has a special
584
meaning; see the code below. <algorithm> and <seed are ignored if
585
<file_string is non-null."""
586
if file_string: # User wants to restore generator state.
587
if (file_string == ' ' or # a space (not the null string)
588
string.strip( string.lower(file_string) ) == 'default'):
589
file_string = 'prvstate.dat'
590
f = open(file_string, 'r') # Get saved state of uniform generator.
591
inlist = f.readlines()
593
self._algorithm_name = inlist[1] # Skip the header (inlist[0]).
594
self._algorithm_name = self._algorithm_name[:-1] # Drop the \n.
595
self._seed = eval(inlist[2])
596
self._count = eval(inlist[3])
597
self._second_nrv = eval(inlist[4])
598
self._index = eval(inlist[5])
600
pos = string.find(self._algorithm_name, ":")
601
if algorithm == 'cmrg':
602
self._fillbuf = self._cmrg # Set generator function pointer.
604
elif algorithm == 'flip':
605
self._fillbuf = self._flip
607
elif algorithm == 'smrg':
608
self._fillbuf = self._smrg
610
elif algorithm == 'twister':
611
self._fillbuf = self._twister
614
pass # can't get here
616
self._build_iterator() # Construct _iterator, and allocate
617
# _series and _ranbuf if necessary.
618
if type(self._series) == array.ArrayType:
620
for i in xrange( len(self._series) ): # 'CMRG', 'Flip', or
621
self._series[i] = eval(inlist[offset + i]) # 'Twister'
623
offset = len(self._series) + 6
624
for i in xrange(self._ranbuf_size):
625
self._ranbuf[i] = eval(inlist[offset + i])
628
self._series = eval( inlist[6] ) # 'SMRG'
630
for i in xrange(self._ranbuf_size):
631
self._ranbuf[i] = eval(inlist[offset + i])
633
else: # User wants a new generator.
634
self._second_nrv = None # Initialize storage for 2nd normal RV.
635
self._count = 0L # Initialize random() delivery count.
636
algorithm = string.strip( string.lower(algorithm) )
637
if algorithm == 'cmrg':
638
self._algorithm_name = \
639
"'CMRG': Combined Multiplicative Recursion MRG32k3a (L'Ecuyer, 1998)"
640
self._fillbuf = self._cmrg # Set generator function pointer.
641
self._build_iterator() # Allocate _series and _ranbuf.
643
self._seed = long(seed) # User wants all initial
644
seed = float((-seed) & 0x7fffffffL) # values equal.
645
for i in range(6): self._series[i] = seed
648
if seed == 0: # Generate initial seed from epoch time.
649
t = (long(time.time()*128) * 1048576L) & 0x7fffffffffffffffL
650
seed = ((t & 0x7fffffffL)^(t >> 32)) & 0x7fffffffL # 31 bits
651
seed = max(seed, 1L) # <seed> must be positive in CMRG
652
seed = min(seed, 4294944442L) # 2**32 - 22853 - 1
654
self._seed = long(seed) # Save for user inquiry.
655
self._series[0] = float(seed)
656
for j in range(1,6): # Use a standard linear multiplicative
657
k = seed / 127773 # congruential generator to get initial
658
seed = 16807 * (seed - k * 127773) - 2836 * k # values for
659
if seed < 0: seed = seed + 0x7fffffff # other 5 multipliers.
660
self._series[j] = float(seed)
662
self._index = self._ranbuf_size-1 # Request new batch of randoms.
664
elif algorithm == 'flip':
665
self._algorithm_name = \
666
"'Flip': Subtractive Series Mod 2**31 (Knuth,1993)"
667
self._fillbuf = self._flip # Set generator function pointer.
668
self._build_iterator() # define_iterator;_ranbuf and _series
669
if seed == 0L: # Get seed from epoch time.
670
t = (long(time.time() * 128) * 1048576L) & 0x7fffffffffffffffL
671
seed = int( (t & 0x7fffffffL)^(t >> 32) ) & 0x7fffffffL # 31 b.
672
seed = max(seed, 1) # Insure seed is positive.
674
self._seed = long(seed) # Save for possible user inquiry.
675
seed = seed & 0x7fffffff # Need positive 4-byte integer here.
676
s = self._series # local alias
677
s[0] = -1 # <_series[0]> is a sentinel unused here.
678
prev = seed # <_series[]> contains signed integers.
684
next = (prev - next) & 0x7fffffff # difference mod 2**31
685
if seed & 0x1: seed = 0x40000000 + (seed >> 1)
686
else: seed = seed >> 1 # cyclic shift right
687
next = (next - seed) & 0x7fffffff # difference mod 2**31
689
i = (i + 21) % 55 # Jump arround in array.
691
self._index = self._ranbuf_size-1 # Request new batch of uniforms.
692
self._random() # Exercize the generator by
693
self._index = self._index + 219 # "using" 4 sets of 55 integers.
695
elif algorithm == 'smrg':
696
self._algorithm_name = \
697
"'SMRG': Single Multiplicative Recursion m61-p3019 (Wu, 1997)"
698
self._fillbuf = self._smrg # Set generator function pointer.
699
self._build_iterator() # Allocate _ranbuf and _series
701
t = (long( time.time() * 128 ) * 544510892281524561475242L)
702
t = t & 0x3ffffffffffffffffffffffffffffffL # Keep 122 bits.
703
seed = ( t ^ (t >> 61) ) & 0x1fffffffffffffffL # 61 bits
704
# The time() multiplier above is arbitrary
705
# so long as t >= (but close to) 2**122.
706
self._seed = seed # Save for user inquiry.
707
seed = abs(seed) # Insure seed is positive.
708
seed = min(seed, 2L**61 - 2) # 2**61 - 1 is the modulus
709
self._series = seed #_series is not an array here; a Python long.
710
self._index = self._ranbuf_size-1 # Request new batch of uniforms.
711
self._random() # Exercize the generator by "using" 16
712
self._index = self._index + 15 # pseudo-random numbers.
714
elif algorithm == 'twister':
715
self._algorithm_name = \
716
"'Twister': Mersenne Twister MT19937 (matsumoto and Nishamura, 1998)"
717
self._fillbuf = self._twister # Set generator function pointer.
718
self._build_iterator() # Define _iterator;_ranbuf, _series.
720
# Construct a seed from epoch time.
721
# 32-bit unsigned integers needed.
722
t = (long(time.time() * 128) * 2097152L) & 0xffffffffffffffffL
723
seed = ( (t & 0xffffffffL)^(t >> 32) ) & 0xffffffffL # 32 bits
725
self._series[0] = seed # <_series[]> contains unsigned integers.
727
self._seed = long(seed) # Save the seed for later user inquiry.
728
seed = seed & 0xffffffffL
730
# Use an algorithm from Knuth(1981) to continue filling <_series>
731
# with 623 more unsigned 32-bit pseudo-random integers.
733
for k in xrange(1, 624):
734
seed = (69069L * seed) & 0xffffffffL # Retain lower 32 bits.
735
self._series[k] = seed
737
self._index = self._ranbuf_size - 1
738
# Force generation of first batch of 624
739
# pseudo-randoms. <_index> points to the
740
# pseudo-random uniform just delivered from
741
# <_ranbuf>; none remain when <_index> is 623.
743
else: raise ValueError, \
744
"algorithm? --must be 'CMRG', 'Flip', 'SMRG', or 'Twister'"
747
def _initial_seed(self):
748
"""Return seed for this series; was either user-supplied or from clock.
752
The result is a Python long integer."""
755
def _random_algorithm(self):
756
"""Return 1-line description of uniform random number algorithm used.
760
The result is a string."""
761
return self._algorithm_name
763
def _random_count(self):
764
"""Return number of uniform(0,1) RVs used since the seed was last reset.
768
A Python long integer is returned. This the number of uniform(0,1)
769
random numbers delivered to the user or consumed by _pranv methods."""
770
return self._count - (self._ranbuf_size - 1 - self._index)
773
def _save_state(self, file_string='prvstate.dat'):
774
"""Save <_pranv> state to disk for later recall and continuation.
776
save_state(file_string='prvstate.dat')
778
A backup file ('filename.bak' or 'prvstate.bak') is created if
779
the file <file_string> can be opened, where <filename*> is
780
the leading portion of <file_string> before the '.', if present."""
783
f = open(file_string, 'r') # If file exists, create backup file.
785
dot_pos = string.rfind(file_string,'.')
786
if dot_pos == -1: bak_file_string = file_string + '.bak'
787
else: bak_file_string = file_string[:dot_pos] + '.bak'
789
f = open(bak_file_string, 'w') # Write backup file.
792
except IOError: # No file was found.
794
f = open(file_string, 'w') # Open file for current state.
795
# Save state. Arrays can't be
796
outlist = [] # pickled, so just write out the data.
797
outlist.append('Python module rv 1.1 random numbers save_state: '
798
+ time.ctime( time.time() ) + '\n') # save_state file header
799
outlist.append(self._algorithm_name + '\n')
800
outlist.append(`self._seed` + '\n')
801
outlist.append(`self._count` + '\n')
802
outlist.append(`self._second_nrv` + '\n')
803
outlist.append(`self._index` + '\n')
804
if type(self._series) is array.ArrayType:
805
for i in xrange( len(self._series) ): # 'CMRG', 'Flip', and 'Twister'
806
outlist.append(`self._series[i]` + '\n')
808
else: # In 'SMRG', _series is a
809
outlist.append(`self._series` + '\n') # single Python long.
811
for i in xrange (self._ranbuf_size): # buffer of uniform(0,1) RVs
812
outlist.append(`self._ranbuf[i]` + '\n')
814
f.writelines(outlist)
817
# Uniform(0,1) Pseudo-random Variate Generator
819
def _random(self, size=None):
820
"""Return a single random uniform(0,1), or <None> and fill buffer.
823
i = self._index # local alias
825
size, Ns = _check_shape(size)
826
buffer = Num.zeros(Ns,Num.Float)
827
buf = self._ranbuf # local alias
828
ranlim = self._ranbuf_size # local alias
829
for j in xrange( Ns ):
831
if i == ranlim: self._fillbuf(); i = 0
833
return Num.reshape(buffer, size)
837
if i == self._ranbuf_size: self._fillbuf(); i = 0
839
return self._ranbuf[i]
841
self._index = i # Restore _ranbuf index and return (buffer filled).
843
def _choice(self, seq=(0,1), size=None):
844
"""Return element k with probability 1/len(seq) from non-empty sequence.
846
choice(seq=(0,1), size=None)
848
Default is a 0, 1 coin flip. <seq> must not be empty. If <buffer> is
849
specified, it must be a mutable sequence. It is filled with random
850
<seq> elements and <None> is returned. Otherwise, a single <seq>
851
element is returned."""
854
raise ValueError, '<seq> must not be empty'
857
i = self._index # local alias
859
size, Ns = _check_shape(size)
860
buffer = Num.zeros(Ns,Num.Float)
861
buflim = self._ranbuf_size # local alias
862
buf = self._ranbuf # local alias
863
for j in xrange( Ns ):
865
if i == buflim: self._fillbuf(); i = 0
867
buffer[j] = seq[int( floor(x * lenseq) )]
868
return Num.reshape(buffer, size)
872
if i == self._ranbuf_size: self._fillbuf(); i = 0
874
self._index = i # Restore updated value of _ranbuf index.
875
return seq[int( floor(x * lenseq) )]
877
self._index = i # Update _ranbuf index and return (buffer[] is filled).
879
def _geom(self, pr=0.5, size=None):
880
"""Return non-negative random integers from a geometric distribution.
882
geom(pr=0.5, size=None)
884
Z has the geometric distribution if it is the number of successes
885
before the first failure in a series of independent Bernouli trials
886
with probability of success 1 - <pr>. 0 <= pr < 1.
887
The result is a non-negative integer less than or equal to 2**31 -1.
888
If <buffer> is specified, it must be a mutable sequence.
889
<buffer> is filled with geometric(pr) pseudo-randoms.
890
Otherwise, a single geometric pseudo-random variate is returned."""
892
pr = 1.0-pr # Added to be consistent with distributions.py
893
if not 0.0 <= pr < 1.0:
894
raise ValueError, '0.0 <= <pr> < 1.0'
896
i = self._index # local alias
897
buf = self._ranbuf # local alias
898
buflim = self._ranbuf_size # local alias
900
size, Ns = _check_shape(size)
901
buffer = Num.zeros(Ns,Num.Float)
909
if pr <= 0.9: # <pr> is small or moderate;
910
while j < out_len: # use inverse transformation.
914
if i == buflim: self._fillbuf(); i = 0
917
successes = successes + 1
922
buffer[j] = successes
926
self._index = i # Restore updated value of _ranbuf index.
929
else: # <pr> larger than 0.9.
932
if i == buflim: self._fillbuf(); i = 0
934
successes = floor( log(u) / log(pr) )
935
successes = int( min(successes, 2147483647.0) ) # 2**31 - 1
937
buffer[j] = successes
941
self._index = i # Restore updated value of _ranbuf index.
944
self._index = i # Restore _ranbuf index and return (buffer[] filled).
945
return Num.reshape(buffer, size)
947
def _hypergeom(self, tot=35, good=25, N=10, size=None):
948
"""Return hypergeometric pseudorandom variates: #"bad" in <sample>
950
hypergeom(tot=35, good=25, N=10)
952
Z has a hypergeometric distribution if it is the number of "bad"
953
(tot-good) items in a sample of size <N> from a population of
954
size tot items. <tot>, <good> and <N> must be positive integers. Also,
955
<tot> >= <sample> >= 1. The result Z satisfies:
956
max(0, <sample> - <good>) <= Z <= min(<sample>, <tot>-<good>).
957
if size is not None is
958
supplied, it must be a mutable sequence (list or array). It is filled
959
with hypergeometric pseudo-random integers and <None> is returned.
960
Otherwise, a single hypergeometric pseudo-random integer is returned.
961
See Fishman, "Monte Carlo," pp 218-221. Algorithms HYP and HRUA* were
964
alpha = int(bad); beta=int(good); n = int(N)
965
if (alpha < 1) or (beta < 1) or (n < 1) or (alpha + beta) < n:
966
raise ValueError, '<bad>, <good>, or <sample> out of range'
969
i = self._index # local alias
970
buf = self._ranbuf # local alias
971
buflim = self._ranbuf_size # local alias
973
size, Ns = _check_shape(size)
974
buffer = Num.zeros(Ns,Num.Float)
982
if n <= 10: # 10 is an empirical estimate.
983
d1 = alpha + beta - n # Use Fishman's HYP algorithm.
984
d2 = float( min(alpha, beta) )
990
if i == buflim: self._fillbuf(); i = 0
992
y = y - floor( u + y / (d1 + k) )
998
if alpha > beta: z = n - z
1007
else: # Use Fishman's HRUA* algorithm.
1008
minalbe = min(alpha, beta)
1009
popsize = alpha + beta
1010
maxalbe = popsize - minalbe
1011
lnfac = self._ln_factorial # function alias
1012
m = min(n, popsize - n)
1013
d1 = 1.715527770 # 2 * sqrt(2 / math.e)
1014
d2 = 0.898916162 # 3 - 2 * sqrt(3 / math.e)
1015
d4 = float(minalbe) / popsize
1018
d7 = sqrt((popsize - m) * n * d4 * d5 / (popsize - 1) + 0.5)
1020
d9 = floor( float((m + 1) * (minalbe + 1)) / (popsize + 2) )
1021
d10 = ( lnfac(d9) + lnfac(minalbe - d9) +
1022
lnfac(m - d9) + lnfac(maxalbe - m + d9) )
1023
d11 = min( min(m, minalbe) + 1.0, floor(d6 + 7 * d7) ) # 7 for 9-
1024
while j < out_len: # digit precision
1025
while 1: # in d1 and d2
1027
if i == buflim: self._fillbuf(); i = 0
1030
if i == buflim: self._fillbuf(); i = 0
1032
w = d6 + d8 * (y - 0.5) / x
1033
if w < 0.0 or w >= d11:
1034
continue # fast rejection; try another x, y
1038
t = d10 - ( lnfac(z) + lnfac(minalbe - z) +
1039
lnfac(m - z) + lnfac(maxalbe - m + z) )
1040
if x * (4.0 - x) - 3.0 <= t:
1041
break # fast acceptance
1044
if x * (x - t) >= 1:
1045
continue # fast rejection, try another x, y
1048
if 2.0 * log(x) <= t:
1051
if alpha > beta: # Error in HRUA*, this is correct.
1054
if m < n: # This fix allows n to exceed
1055
z = alpha - z # popsize / 2.
1065
self._index = i # Restore ranbuf index and return (buffer is filled).
1066
return Num.reshape(buffer, size)
1069
def _logser(self, pr=0.5, size=None):
1070
"""Return logarithmic (log-series) positive pseudo-random integers;
1073
logser(pr=0.5, size=None)
1075
Z has the logseries(logarithmic)distribution if Pr{Z=i} = (a/i) * p ** i,
1076
for i = 1, 2, ... and a = -1.0 / log(1.0 - p). If <buffer> is specified,
1077
it must be a mutable sequence (list or array). It is filled with
1078
logseries (positive) pseudo-random integers and <None> is returned.
1079
Otherwise, a single logarithmic series pseudo-random is returned.
1080
The algorithm is by A.W. Kemp; see Devroye, 1986, pp 547-548."""
1082
if not(0.0 < p < 1.0):
1083
raise ValueError, '<p> must be in (0.0, 1.0)'
1087
i = self._index # local alias
1088
buf = self._ranbuf # local alias
1089
buflim = self._ranbuf_size # local alias
1090
if size is not None:
1091
size, Ns = _check_shape(size)
1092
buffer = Num.zeros(Ns,Num.Float)
1105
if i == buflim: self._fillbuf(); i = 0
1111
sum = sum * p * (k - 1) / k
1121
else: # p is >= 0.95.
1125
if i == buflim: self._fillbuf(); i = 0
1129
if i == buflim: self._fillbuf(); i = 0
1131
q = 1.0 - exp(r * u)
1134
k = int( floor(1.0 + log(v) / log(q)) )
1150
self._index = i # Restore ranbuf index and return (<buffer> filled).
1151
return Num.reshape(buffer, size)
1153
def _von_Mises(self, b, loc=0.0, size=None):
1154
"""Return von Mises distribution pseudo-random variates on [-pi, +pi].
1156
von_Mises(b, loc=0.0, size=None)
1158
<mean> must be in the open interval (-math.pi, +math.pi). If <buffer>
1159
is specified, it must be a mutable sequence (list or array). It is
1160
filled with von Mises(mean, shape) pseudo-random variates and <None> is
1161
returned. Otherwise, a single von Mises RV is returned. The method is
1162
an algorithm of Best and Fisher, 1979; see Fisher, N. I., "Statistical
1163
Analysis of Circular Data," Cambridge University Press, 1995, p. 49."""
1165
shape, mean = b, loc
1167
mean = arctan2(z.imag, z.real)
1168
if not (-3.1415926535897931 < mean < +3.1415926535897931):
1170
'<loc> must be in the open interval (-math.pi, math.pi)'
1173
a = 1.0 + sqrt( 1.0 + 4.0 * shape*shape )
1174
b = ( a - sqrt(a + a) ) / (shape + shape)
1175
r = (1.0 + b * b) / (b + b)
1176
i = self._index # local alias
1177
buf = self._ranbuf # local alias
1178
buflim = self._ranbuf_size # local alias
1179
if size is not None:
1180
size, Ns = _check_shape(size)
1181
buffer = Num.zeros(Ns,Num.Float)
1192
if i == buflim: self._fillbuf(); i = 0
1193
z = cos(3.1415926535897931 * buf[i])
1194
f = (1.0 + r * z) / (r + z)
1197
if i == buflim: self._fillbuf(); i = 0
1199
if (c * (2.0 - c) - u > 0.0):
1200
break # quick acceptance
1202
if log(c / u) + 1.0 - c < 0.0:
1203
continue # quick rejection
1209
if i == buflim: self._fillbuf(); i = 0
1212
buffer[j] = fmod(acos(f) + mean, 6.2831853071795862)
1215
buffer[j] = -fmod(acos(f) + mean, 6.2831853071795862)
1220
self._index = i # Restore updated value of _ranbuf index.
1221
if buf[i] > 0.5: return fmod(acos(f)+mean, 6.2831853071795862)
1222
else: return -fmod(acos(f)+mean, 6.2831853071795862)
1224
self._index = i # Restore _ranbuf index and return (buffer[] filled).
1225
return Num.reshape(buffer, size)
1227
def _Wald(self, mu, loc=0.0, scale=1.0, size=None):
1228
"""Return Inverse gaussian pseudo-random variates.
1230
invnorm(mu, loc=0.0, scale=1.0, size=None)
1232
Both <mean> and <scale> must be positive. If <buffer> is specified, it
1233
must be a mutable sequence (list or array). It is filled with inverse
1234
gaussian (wald if mu=1.0) pseudo-random variates, and <None> is returned.
1235
Otherwise, a single Wald pseudo-random variate is returned."""
1239
if (mean <= 0.0) or (scale <= 0.0):
1240
raise ValueError, 'both <mean> and <scale> must be positive'
1243
i = self._index # local alias
1244
buf = self._ranbuf # local alias
1245
buflim = self._ranbuf_size # local alias
1246
if size is not None:
1247
size, Ns = _check_shape(size)
1248
buffer = Num.zeros(Ns,Num.Float)
1256
normal = self._normal # local alias
1257
r = 0.5 * mean / scale
1259
meansq = mean * mean
1263
x1 = mean + r * ( muy - sqrt(muy * (con + muy)) )
1265
if i == buflim: self._fillbuf(); i = 0
1267
if buf[i] <= mean / (mean + x1):
1271
buffer[j] = meansq / x1
1276
self._index = i # Restore updated value of _ranbuf index.
1277
if buf[i] <= mean / (mean + x1):
1283
self._index = i # Restore _ranbuf index and return (buffer[] filled).
1284
return Num.reshape(buffer*oscale+loc, size)
1286
def _Zipf(self, a=4.0, size=None):
1287
"""Return positive pseudo-random integers from the Zipf distribution
1289
Zipf(a=4.0, size=None)
1291
Z has the Zipf ditribution with parameter a if
1292
Pr{Z=i} = 1.0 / (zeta(a) * i ** a) for i = 1, 2, ... and a > 1.0.
1293
Here zeta(a) is the Riemann zeta function, the sum from 1 to infinity of
1294
1.0 / i **a. If <buffer> is specified, it must be a mutable sequence
1295
(list or array). It is filled with Zipf pseudo-random integers and
1296
<None> is returned. Otherwise, a single (positive) Zipf pseudo-random
1297
integer is returned; see Devroye, 1986, pp 550-551 for the algorithm."""
1299
raise ValueError, '<a> must be larger than 1.0'
1302
i = self._index # local alias
1303
buf = self._ranbuf # local alias
1304
buflim = self._ranbuf_size # local alias
1307
if size is not None:
1308
size, Ns = _check_shape(size)
1309
buffer = Num.zeros(Ns,Num.Float)
1320
if i == buflim: self._fillbuf(); i = 0
1321
x = floor( buf[i] ** (-1 / am1) )
1322
t = (1.0 + 1.0 / x) ** am1
1324
if i == buflim: self._fillbuf(); i = 0
1325
if buf[i] * x * (t - 1.0) / (b - 1.0) <= t / b:
1336
self._index = i # Restore ranbuf index and return (buffer[] filled).
1337
return Num.reshape(buffer, size)
1339
# End of Non-uniform(0,1) Generators
1341
# Geometric and Subsampling Routines
1343
def _in_simplex(self, mseq=5*[0.0], bound=1.0):
1344
"""Return a pseudorandom point in a simplex.
1346
in_simplex(mseq=2*[0.0], bound=1.0)
1348
The simplex of dimension d and bound b > 0 is defined by all points
1349
z = (z1, z2, ..., zd) with z1 + z2 + ... + zd <= b and all zi >= 0.
1350
<bound> must be positive.
1351
A two-dimensional simplex with bound b is the region enclosed by the
1352
horizontal (x) and vertical (y) axes and the line y = b - x. A
1353
one dimensional simplex with bound b is the single point b.
1354
<mseq> must be a mutable sequence (list or array); if it is an array, it
1355
must be typed double or single float. A sequence of the same type and
1356
length as <mseq> containing the coordinates of a pseudo-random point in
1357
the simplex is returned.
1358
See Fishman, George, "Monte Carlo," Springer-Verlag, 1996, pp. 232-233
1359
and 226-229. Algorithms EUNIFORM and OSIMP were used."""
1363
raise ValueError, '<bound> must be positive.'
1369
i = self._index # local alias
1370
buf = self._ranbuf # local alias
1371
buflim = self._ranbuf_size # local alias
1374
if i == buflim: self._fillbuf(); i = 0
1375
point[0] = -log(buf[i])
1376
for j in range_d[1:]: # Accumulate running sum of exponential(1) RVs.
1378
if i == buflim: self._fillbuf(); i = 0
1379
point[j] = point[j-1] - log(buf[i])
1382
if i == buflim: self._fillbuf(); i = 0
1383
total = point[d-1] - log(buf[i])
1385
point[j] = point[j] / total
1387
point[0] = point[0] * bound
1390
point[j] = bound * (point[j] - point[j-1])
1393
self._index = i # Restore updated value of _ranbuf index and return.
1397
def _in_sphere(self, center=5*[0.0], radius=1.0):
1398
"""Return a pseudo-random point in a sphere centered at <center>.
1400
in_sphere(center=5*[0.0], radius=1.0)
1402
A one-dimensional sphere centered at 0 is the two points <+radius> and
1403
<-radius>. A two-dimensional sphere centered at 0 is the circle with
1404
radius <radius>. <center> must be a sequence containing double
1405
coordinates for the sphere's center. <radius> must be positive.
1406
pseudo-random point in the "sphere" with center 0 and radius <radius>.
1407
See Fishman, George, "Monte Carlo," Springer-Verlag, 1996, pp 234-235.
1408
Algorithm OSPHERE was used."""
1411
raise ValueError, '<radius> must be positive.'
1414
return self._on_sphere( center, radius * self._beta(float(d), 1.0) )
1416
def _on_simplex(self, mseq=5*[0.0], bound=1.0):
1417
"""Return a pseudo_random point on the boundary of a simplex.
1419
on_simplex(mseq=5*[0.0], bound=1.0)
1421
The simplex of dimension d and bound b > 0 is defined by all vectors
1422
z = (z1, z2, ..., zd) with z1 + z2 + ... + zd <= b and all zi >= 0.
1423
<mseq> must be a mutable sequence (list or array) and <bound> must be
1424
positive. A boundary point on the simplex has the sum of its
1425
coefficients equal to b.
1426
The coordinates of a pseudo-random point on the boundary of the
1427
simplex of dimension len(seq) and bound <bound> are returned as a
1428
sequence of the same type and length as <mseq>. The input values in
1429
<mseq> are ignored, and undisturbed.
1430
See Fishman, George, "Monte Carlo," Springer-Verlag, 1996,
1431
pp 232-233 and 226-229. Algorithms EUNIFORM and OSIMP were used."""
1435
raise ValueError, '<bound> must be positive.'
1437
elif d == 1: # The point is just <bound>.
1440
i = self._index # local alias
1441
buf = self._ranbuf # local alias
1442
buflim = self._ranbuf_size # local alias
1443
range_d = range(d) # [0,1,...,d-1]
1444
range_dm1 = range_d[:(d - 1)] # [0,1,...,d-2]
1445
range_1_dm1 = range_dm1[1:] # [1,2,...,d-2]
1447
if i == buflim: self._fillbuf(); i = 0
1448
point[0] = -log(buf[i])
1449
for j in range_1_dm1: # range_1_dm1 is empty for d = 2.
1451
if i == buflim: self._fillbuf(); i = 0
1452
point[j] = point[j - 1] - log(buf[i])
1455
if i == buflim: self._fillbuf(); i = 0
1456
total = point[d - 2] - log(buf[i])
1458
point[j] = point[j] / total
1461
point[0] = bound * point[0]
1463
while j > 0: # not executed for d = 2.
1464
point[j] = bound * (point[j] - point[j - 1])
1467
point[d - 1] = bound - last
1469
self._index = i # Restore updated value of _ranbuf index and return.
1472
def _on_sphere(self, center=5*[0.0], radius=1.0):
1473
"""Return a pseudo-random point on the sphere centered at point[].
1475
on_sphere(center=5*[0.0], radius=1.0)
1477
The sphere of dimension d centered at 0 is all points (x1,x2,...,xd)
1478
with x1**2 + x2**2 + ... + xd**2 <= <radius>**2. The surface of this
1479
sphere is any point with the sum of its components squared equal to
1480
the radius squared. A 2-dimensional sphere is a circle, and a
1481
1-dimensional sphere is the line from <-radius> to <+radius>. <center>
1482
must be a sequence and <radius> must be positive. A sequence of the
1483
same type as <center[]> containing the coordinates of a pseudo-random
1484
point on the "sphere" centered at <center[]> is returned.
1485
See Fishman, George, "Monte Carlo," Springer-Verlag, 1996, pp 234-235.
1486
Algorithm OSPHERE was used. For dimensions 2, 3, and 4, rejection
1487
algorithms were used. See Marsaglia, G., "Choosing a point from the
1488
surface of a sphere," Annals of mathematical Statistics, vol. 43,
1489
pp. 645-646, 1972."""
1491
buflim = self._ranbuf_size # local alias
1492
buf = self._ranbuf # local alias
1493
i = self._index # local alias
1496
raise ValueError, '<radius> must be positive.'
1498
elif d == 0: # Just return center[:].
1503
if i == buflim: self._fillbuf(); i = 0
1505
point[0] = radius + center[0]
1508
point[0] = -radius + center[0]
1510
elif d == 2: # rejection from circumscribed square.
1514
if i == buflim: self._fillbuf(); i = 0
1517
if i == buflim: self._fillbuf(); i = 0
1519
x = x + x - 1.0 # (x, y) is in the square circum-
1520
y = y + y - 1.0 # scribed on the unit circle
1521
ss = x * x + y * y # centered at (0, 0).
1523
ss = radius / sqrt(ss) # (x, y) now restricted to unit
1524
point[0] = x * ss + center[0] # circle.
1525
point[1] = y * ss + center[1]
1526
self._index = i # Restore ranbuf index and return.
1528
elif d == 3: # algorithm by Marsaglia, 1972
1532
if i == buflim: self._fillbuf(); i = 0
1535
if i == buflim: self._fillbuf(); i = 0
1541
con = sqrt(1.0 - ss) * radius
1542
point[0] = 2.0 * u * con + center[0]
1543
point[1] = 2.0 * v * con + center[1]
1544
point[2] = (1.0 - ss - ss) * radius + center[2]
1545
self._index = i # Restore ranbuf index and return.
1548
ssuv = 2.0 # algorithm by Marsaglia, 1972
1551
if i == buflim: self._fillbuf(); i = 0
1554
if i == buflim: self._fillbuf(); i = 0
1558
ssuv = u * u + v * v
1563
if i == buflim: self._fillbuf(); i = 0
1566
if i == buflim: self._fillbuf(); i = 0
1570
ssrs = r * r + s * s
1572
con = sqrt((1.0 - ssuv) / ssrs) * radius
1573
point[0] = u * radius + center[0]
1574
point[1] = v * radius + center[1]
1575
point[2] = r * con + center[2]
1576
point[3] = s * con + center[3]
1577
self._index = i # Restore ranbuf index and return.
1579
else: # Use radially symmetric normals
1580
range_d = range(d) # (Fishmans algorithm OSPHERE).
1582
normal = self._normal
1588
constant = radius / sqrt(sum)
1590
point[j] = constant * point[j] + center[j]
1594
def _sample(self, inlist=[0,1,2], sample_size=2):
1595
"""Return simple random sample of <sample_size> items from <inlist>.
1597
sample(inlist=[0,1,2], sample_size=2)
1599
<inlist> must be a mutable sequence (list or array). The
1600
<sample_size> <= len(inlist) elements of the returned mutable sequence
1601
(which has the same type as <inlist>) are a simple random sample from
1602
all the elements in <inlist>.
1603
Sampling is with replacement; duplicates may occur in the returned
1604
sequence. The length of the returned list or array is <sample_size>,
1605
which may be any positive integer."""
1606
pop_size = len(inlist)
1607
n = int(sample_size)
1612
outlist = 0 * inlist[:1]
1615
i = self._index # local alias
1616
buf = self._ranbuf # local alias
1617
buflim = self._ranbuf_size # local alias
1618
outlist = n*inlist[:1] # <sample_size> >= 1
1621
if i == buflim: self._fillbuf(); i = 0
1622
outlist[j] = inlist[int( floor(buf[i] * pop_size) )]
1624
self._index = i # Restore _ranbuf index arnd return.
1627
def _smart_sample(self, inlist=[0,1,2], sample_size=2):
1628
"""Return a random sample without replacement of elements in <inlist>.
1630
smart_sample(inlist=[0,1,2], sample_size=2)
1632
<inlist> must be a mutable sequence (list or array). The
1633
<sample_size> <= len(inlist) elements of the returned mutable sequence
1634
(which has the same type as <inlist>) are a sample without replacement
1635
of all the elements in <inlist>. Since sampling is without replacement,
1636
there is no possibility of duplicates in the returned list or array
1637
unless there are duplicates in <inlist>. The length of the
1638
returned mutable sequence is the smaller of <sample_size> and
1639
len(inlist), and its type is the type of <inlist>."""
1640
pop_size = len(inlist)
1641
n = int( min(sample_size, pop_size) ) # can't have sample > population.
1646
outlist = 0*inlist[:1]
1649
outlist = inlist[:] # <outlist> has same type as <inlist>.
1650
i = self._index # local alias
1651
buf = self._ranbuf # local alias
1652
buflim = self._ranbuf_size # local alias
1654
temp = outlist[j] # Save item in ith position.
1656
if i == buflim: self._fillbuf(); i = 0
1657
k = j + int( floor(buf[i] * (pop_size - j)) )
1658
outlist[j] = outlist[k] # Move sampled item to ith position.
1659
outlist[k] = temp # copy to sampled item's position.
1662
del outlist[n:pop_size] # Remove any unsampled elements.
1664
self._index = i # Restore _ranbuf index and return.
1667
_inst = _pranv() # Initialize uniform(0,1) generator to default; clock seed.
1669
initial_seed = _inst._initial_seed
1670
initialize = _inst._initialize
1671
random_algorithm = _inst._random_algorithm
1672
random_count = _inst._random_count
1673
save_state = _inst._save_state
1675
choice = _inst._choice
1677
# Geometrical Point Generators, Permutations and Subsampling Routines
1678
in_simplex = _inst._in_simplex
1679
in_sphere = _inst._in_sphere
1680
on_simplex = _inst._on_simplex
1681
on_sphere = _inst._on_sphere
1682
sample = _inst._sample
1683
smart_sample = _inst._smart_sample
1686
##def mean_var_test(x, type, mean, var, skew=[]):
1688
## x_mean = Num.sum(x)/n
1689
## x_minus_mean = x - x_mean
1690
## x_var = Num.sum(x_minus_mean*x_minus_mean)/(n-1.0)
1691
## print "\nAverage of ", len(x), type
1692
## print "(should be about ", mean, "):", x_mean
1693
## print "Variance of those random numbers (should be about ", var, "):", x_var
1695
## x_skew = (Num.sum(x_minus_mean*x_minus_mean*x_minus_mean)/9998.)/x_var**(3./2.)
1696
## print "Skewness of those random numbers (should be about ", skew, "):", x_skew
1699
## x, y = get_seed()
1700
## print "Initial seed", x, y
1702
## x1, y1 = get_seed()
1703
## if x1 != x or y1 != y:
1704
## raise SystemExit, "Failed seed test."
1705
## print "First random number is", random()
1706
## print "Average of 10000 random numbers is", Num.sum(random(10000))/10000.
1707
## x = random([10,1000])
1708
## if len(x.shape) != 2 or x.shape[0] != 10 or x.shape[1] != 1000:
1709
## raise SystemExit, "random returned wrong shape"
1710
## x.shape = (10000,)
1711
## print "Average of 100 by 100 random numbers is", Num.sum(x)/10000.
1712
## y = uniform(0.5,0.6, (1000,10))
1713
## if len(y.shape) !=2 or y.shape[0] != 1000 or y.shape[1] != 10:
1714
## raise SystemExit, "uniform returned wrong shape"
1715
## y.shape = (10000,)
1716
## if Num.minimum.reduce(y) <= 0.5 or Num.maximum.reduce(y) >= 0.6:
1717
## raise SystemExit, "uniform returned out of desired range"
1718
## print "randint(1, 10, size=[50])"
1719
## print randint(1, 10, size=[50])
1720
## print "permutation(10)", permutation(10)
1721
## print "randint(3,9)", randint(3,9)
1722
## print "random_integers(10, size=[20])"
1723
## print random_integers(10, size=[20])
1725
## x = norm(2.0, s, [10, 1000])
1726
## if len(x.shape) != 2 or x.shape[0] != 10 or x.shape[1] != 1000:
1727
## raise SystemExit, "standard_normal returned wrong shape"
1728
## x.shape = (10000,)
1729
## mean_var_test(x, "normally distributed numbers with mean 2 and variance %f"%(s**2,), 2, s**2, 0)
1730
## x = exponential(3, 10000)
1731
## mean_var_test(x, "random numbers exponentially distributed with mean %f"%(s,), s, s**2, 2)
1732
## x = multivariate_normal(Num.array([10,20]), Num.array(([1,2],[2,4])))
1733
## print "\nA multivariate normal", x
1734
## if x.shape != (2,): raise SystemExit, "multivariate_normal returned wrong shape"
1735
## x = multivariate_normal(Num.array([10,20]), Num.array([[1,2],[2,4]]), [4,3])
1736
## print "A 4x3x2 array containing multivariate normals"
1738
## if x.shape != (4,3,2): raise SystemExit, "multivariate_normal returned wrong shape"
1739
## x = multivariate_normal(Num.array([-100,0,100]), Num.array([[3,2,1],[2,2,1],[1,1,1]]), 10000)
1740
## x_mean = Num.sum(x)/10000.
1741
## print "Average of 10000 multivariate normals with mean [-100,0,100]"
1743
## x_minus_mean = x - x_mean
1744
## print "Estimated covariance of 10000 multivariate normals with covariance [[3,2,1],[2,2,1],[1,1,1]]"
1745
## print Num.matrixmultiply(Num.transpose(x_minus_mean),x_minus_mean)/9999.
1746
## x = beta(5.0, 10.0, 10000)
1747
## mean_var_test(x, "beta(5.,10.) random numbers", 0.333, 0.014)
1748
## x = gamma(.01, 2., 10000)
1749
## mean_var_test(x, "gamma(.01,2.) random numbers", 2*100, 2*100*100)
1750
## x = chi_square(11., 10000)
1751
## mean_var_test(x, "chi squared random numbers with 11 degrees of freedom", 11, 22, 2*Num.sqrt(2./11.))
1752
## x = F(5., 10., 10000)
1753
## mean_var_test(x, "F random numbers with 5 and 10 degrees of freedom", 1.25, 1.35)
1754
## x = poisson(50., 10000)
1755
## mean_var_test(x, "poisson random numbers with mean 50", 50, 50, 0.14)
1756
## print "\nEach element is the result of 16 binomial trials with probability 0.5:"
1757
## print binomial(16, 0.5, 16)
1758
## print "\nEach element is the result of 16 negative binomial trials with probability 0.5:"
1759
## print negative_binomial(16, 0.5, [16,])
1760
## print "\nEach row is the result of 16 multinomial trials with probabilities [0.1, 0.5, 0.1 0.3]:"
1761
## x = multinomial(16, [0.1, 0.5, 0.1], 8)
1763
## print "Mean = ", Num.sum(x)/8.
1765
##if __name__ == '__main__':
added
removed
Lib/stats/rv.py
