Edit File by line
/home/barbar84/www/wp-conte.../plugins/sujqvwi/ExeBy/exe_root.../opt/alt/python27/lib64/python2....
File: random.py
"""Random variable generators.
[0] Fix | Delete
[1] Fix | Delete
integers
[2] Fix | Delete
--------
[3] Fix | Delete
uniform within range
[4] Fix | Delete
[5] Fix | Delete
sequences
[6] Fix | Delete
---------
[7] Fix | Delete
pick random element
[8] Fix | Delete
pick random sample
[9] Fix | Delete
generate random permutation
[10] Fix | Delete
[11] Fix | Delete
distributions on the real line:
[12] Fix | Delete
------------------------------
[13] Fix | Delete
uniform
[14] Fix | Delete
triangular
[15] Fix | Delete
normal (Gaussian)
[16] Fix | Delete
lognormal
[17] Fix | Delete
negative exponential
[18] Fix | Delete
gamma
[19] Fix | Delete
beta
[20] Fix | Delete
pareto
[21] Fix | Delete
Weibull
[22] Fix | Delete
[23] Fix | Delete
distributions on the circle (angles 0 to 2pi)
[24] Fix | Delete
---------------------------------------------
[25] Fix | Delete
circular uniform
[26] Fix | Delete
von Mises
[27] Fix | Delete
[28] Fix | Delete
General notes on the underlying Mersenne Twister core generator:
[29] Fix | Delete
[30] Fix | Delete
* The period is 2**19937-1.
[31] Fix | Delete
* It is one of the most extensively tested generators in existence.
[32] Fix | Delete
* Without a direct way to compute N steps forward, the semantics of
[33] Fix | Delete
jumpahead(n) are weakened to simply jump to another distant state and rely
[34] Fix | Delete
on the large period to avoid overlapping sequences.
[35] Fix | Delete
* The random() method is implemented in C, executes in a single Python step,
[36] Fix | Delete
and is, therefore, threadsafe.
[37] Fix | Delete
[38] Fix | Delete
"""
[39] Fix | Delete
[40] Fix | Delete
from __future__ import division
[41] Fix | Delete
from warnings import warn as _warn
[42] Fix | Delete
from types import MethodType as _MethodType, BuiltinMethodType as _BuiltinMethodType
[43] Fix | Delete
from math import log as _log, exp as _exp, pi as _pi, e as _e, ceil as _ceil
[44] Fix | Delete
from math import sqrt as _sqrt, acos as _acos, cos as _cos, sin as _sin
[45] Fix | Delete
from os import urandom as _urandom
[46] Fix | Delete
from binascii import hexlify as _hexlify
[47] Fix | Delete
import hashlib as _hashlib
[48] Fix | Delete
[49] Fix | Delete
__all__ = ["Random","seed","random","uniform","randint","choice","sample",
[50] Fix | Delete
"randrange","shuffle","normalvariate","lognormvariate",
[51] Fix | Delete
"expovariate","vonmisesvariate","gammavariate","triangular",
[52] Fix | Delete
"gauss","betavariate","paretovariate","weibullvariate",
[53] Fix | Delete
"getstate","setstate","jumpahead", "WichmannHill", "getrandbits",
[54] Fix | Delete
"SystemRandom"]
[55] Fix | Delete
[56] Fix | Delete
NV_MAGICCONST = 4 * _exp(-0.5)/_sqrt(2.0)
[57] Fix | Delete
TWOPI = 2.0*_pi
[58] Fix | Delete
LOG4 = _log(4.0)
[59] Fix | Delete
SG_MAGICCONST = 1.0 + _log(4.5)
[60] Fix | Delete
BPF = 53 # Number of bits in a float
[61] Fix | Delete
RECIP_BPF = 2**-BPF
[62] Fix | Delete
[63] Fix | Delete
[64] Fix | Delete
# Translated by Guido van Rossum from C source provided by
[65] Fix | Delete
# Adrian Baddeley. Adapted by Raymond Hettinger for use with
[66] Fix | Delete
# the Mersenne Twister and os.urandom() core generators.
[67] Fix | Delete
[68] Fix | Delete
import _random
[69] Fix | Delete
[70] Fix | Delete
class Random(_random.Random):
[71] Fix | Delete
"""Random number generator base class used by bound module functions.
[72] Fix | Delete
[73] Fix | Delete
Used to instantiate instances of Random to get generators that don't
[74] Fix | Delete
share state. Especially useful for multi-threaded programs, creating
[75] Fix | Delete
a different instance of Random for each thread, and using the jumpahead()
[76] Fix | Delete
method to ensure that the generated sequences seen by each thread don't
[77] Fix | Delete
overlap.
[78] Fix | Delete
[79] Fix | Delete
Class Random can also be subclassed if you want to use a different basic
[80] Fix | Delete
generator of your own devising: in that case, override the following
[81] Fix | Delete
methods: random(), seed(), getstate(), setstate() and jumpahead().
[82] Fix | Delete
Optionally, implement a getrandbits() method so that randrange() can cover
[83] Fix | Delete
arbitrarily large ranges.
[84] Fix | Delete
[85] Fix | Delete
"""
[86] Fix | Delete
[87] Fix | Delete
VERSION = 3 # used by getstate/setstate
[88] Fix | Delete
[89] Fix | Delete
def __init__(self, x=None):
[90] Fix | Delete
"""Initialize an instance.
[91] Fix | Delete
[92] Fix | Delete
Optional argument x controls seeding, as for Random.seed().
[93] Fix | Delete
"""
[94] Fix | Delete
[95] Fix | Delete
self.seed(x)
[96] Fix | Delete
self.gauss_next = None
[97] Fix | Delete
[98] Fix | Delete
def seed(self, a=None):
[99] Fix | Delete
"""Initialize internal state of the random number generator.
[100] Fix | Delete
[101] Fix | Delete
None or no argument seeds from current time or from an operating
[102] Fix | Delete
system specific randomness source if available.
[103] Fix | Delete
[104] Fix | Delete
If a is not None or is an int or long, hash(a) is used instead.
[105] Fix | Delete
Hash values for some types are nondeterministic when the
[106] Fix | Delete
PYTHONHASHSEED environment variable is enabled.
[107] Fix | Delete
"""
[108] Fix | Delete
[109] Fix | Delete
if a is None:
[110] Fix | Delete
try:
[111] Fix | Delete
# Seed with enough bytes to span the 19937 bit
[112] Fix | Delete
# state space for the Mersenne Twister
[113] Fix | Delete
a = long(_hexlify(_urandom(2500)), 16)
[114] Fix | Delete
except NotImplementedError:
[115] Fix | Delete
import time
[116] Fix | Delete
a = long(time.time() * 256) # use fractional seconds
[117] Fix | Delete
[118] Fix | Delete
super(Random, self).seed(a)
[119] Fix | Delete
self.gauss_next = None
[120] Fix | Delete
[121] Fix | Delete
def getstate(self):
[122] Fix | Delete
"""Return internal state; can be passed to setstate() later."""
[123] Fix | Delete
return self.VERSION, super(Random, self).getstate(), self.gauss_next
[124] Fix | Delete
[125] Fix | Delete
def setstate(self, state):
[126] Fix | Delete
"""Restore internal state from object returned by getstate()."""
[127] Fix | Delete
version = state[0]
[128] Fix | Delete
if version == 3:
[129] Fix | Delete
version, internalstate, self.gauss_next = state
[130] Fix | Delete
super(Random, self).setstate(internalstate)
[131] Fix | Delete
elif version == 2:
[132] Fix | Delete
version, internalstate, self.gauss_next = state
[133] Fix | Delete
# In version 2, the state was saved as signed ints, which causes
[134] Fix | Delete
# inconsistencies between 32/64-bit systems. The state is
[135] Fix | Delete
# really unsigned 32-bit ints, so we convert negative ints from
[136] Fix | Delete
# version 2 to positive longs for version 3.
[137] Fix | Delete
try:
[138] Fix | Delete
internalstate = tuple( long(x) % (2**32) for x in internalstate )
[139] Fix | Delete
except ValueError, e:
[140] Fix | Delete
raise TypeError, e
[141] Fix | Delete
super(Random, self).setstate(internalstate)
[142] Fix | Delete
else:
[143] Fix | Delete
raise ValueError("state with version %s passed to "
[144] Fix | Delete
"Random.setstate() of version %s" %
[145] Fix | Delete
(version, self.VERSION))
[146] Fix | Delete
[147] Fix | Delete
def jumpahead(self, n):
[148] Fix | Delete
"""Change the internal state to one that is likely far away
[149] Fix | Delete
from the current state. This method will not be in Py3.x,
[150] Fix | Delete
so it is better to simply reseed.
[151] Fix | Delete
"""
[152] Fix | Delete
# The super.jumpahead() method uses shuffling to change state,
[153] Fix | Delete
# so it needs a large and "interesting" n to work with. Here,
[154] Fix | Delete
# we use hashing to create a large n for the shuffle.
[155] Fix | Delete
s = repr(n) + repr(self.getstate())
[156] Fix | Delete
n = int(_hashlib.new('sha512', s).hexdigest(), 16)
[157] Fix | Delete
super(Random, self).jumpahead(n)
[158] Fix | Delete
[159] Fix | Delete
## ---- Methods below this point do not need to be overridden when
[160] Fix | Delete
## ---- subclassing for the purpose of using a different core generator.
[161] Fix | Delete
[162] Fix | Delete
## -------------------- pickle support -------------------
[163] Fix | Delete
[164] Fix | Delete
def __getstate__(self): # for pickle
[165] Fix | Delete
return self.getstate()
[166] Fix | Delete
[167] Fix | Delete
def __setstate__(self, state): # for pickle
[168] Fix | Delete
self.setstate(state)
[169] Fix | Delete
[170] Fix | Delete
def __reduce__(self):
[171] Fix | Delete
return self.__class__, (), self.getstate()
[172] Fix | Delete
[173] Fix | Delete
## -------------------- integer methods -------------------
[174] Fix | Delete
[175] Fix | Delete
def randrange(self, start, stop=None, step=1, _int=int, _maxwidth=1L<<BPF):
[176] Fix | Delete
"""Choose a random item from range(start, stop[, step]).
[177] Fix | Delete
[178] Fix | Delete
This fixes the problem with randint() which includes the
[179] Fix | Delete
endpoint; in Python this is usually not what you want.
[180] Fix | Delete
[181] Fix | Delete
"""
[182] Fix | Delete
[183] Fix | Delete
# This code is a bit messy to make it fast for the
[184] Fix | Delete
# common case while still doing adequate error checking.
[185] Fix | Delete
istart = _int(start)
[186] Fix | Delete
if istart != start:
[187] Fix | Delete
raise ValueError, "non-integer arg 1 for randrange()"
[188] Fix | Delete
if stop is None:
[189] Fix | Delete
if istart > 0:
[190] Fix | Delete
if istart >= _maxwidth:
[191] Fix | Delete
return self._randbelow(istart)
[192] Fix | Delete
return _int(self.random() * istart)
[193] Fix | Delete
raise ValueError, "empty range for randrange()"
[194] Fix | Delete
[195] Fix | Delete
# stop argument supplied.
[196] Fix | Delete
istop = _int(stop)
[197] Fix | Delete
if istop != stop:
[198] Fix | Delete
raise ValueError, "non-integer stop for randrange()"
[199] Fix | Delete
width = istop - istart
[200] Fix | Delete
if step == 1 and width > 0:
[201] Fix | Delete
# Note that
[202] Fix | Delete
# int(istart + self.random()*width)
[203] Fix | Delete
# instead would be incorrect. For example, consider istart
[204] Fix | Delete
# = -2 and istop = 0. Then the guts would be in
[205] Fix | Delete
# -2.0 to 0.0 exclusive on both ends (ignoring that random()
[206] Fix | Delete
# might return 0.0), and because int() truncates toward 0, the
[207] Fix | Delete
# final result would be -1 or 0 (instead of -2 or -1).
[208] Fix | Delete
# istart + int(self.random()*width)
[209] Fix | Delete
# would also be incorrect, for a subtler reason: the RHS
[210] Fix | Delete
# can return a long, and then randrange() would also return
[211] Fix | Delete
# a long, but we're supposed to return an int (for backward
[212] Fix | Delete
# compatibility).
[213] Fix | Delete
[214] Fix | Delete
if width >= _maxwidth:
[215] Fix | Delete
return _int(istart + self._randbelow(width))
[216] Fix | Delete
return _int(istart + _int(self.random()*width))
[217] Fix | Delete
if step == 1:
[218] Fix | Delete
raise ValueError, "empty range for randrange() (%d,%d, %d)" % (istart, istop, width)
[219] Fix | Delete
[220] Fix | Delete
# Non-unit step argument supplied.
[221] Fix | Delete
istep = _int(step)
[222] Fix | Delete
if istep != step:
[223] Fix | Delete
raise ValueError, "non-integer step for randrange()"
[224] Fix | Delete
if istep > 0:
[225] Fix | Delete
n = (width + istep - 1) // istep
[226] Fix | Delete
elif istep < 0:
[227] Fix | Delete
n = (width + istep + 1) // istep
[228] Fix | Delete
else:
[229] Fix | Delete
raise ValueError, "zero step for randrange()"
[230] Fix | Delete
[231] Fix | Delete
if n <= 0:
[232] Fix | Delete
raise ValueError, "empty range for randrange()"
[233] Fix | Delete
[234] Fix | Delete
if n >= _maxwidth:
[235] Fix | Delete
return istart + istep*self._randbelow(n)
[236] Fix | Delete
return istart + istep*_int(self.random() * n)
[237] Fix | Delete
[238] Fix | Delete
def randint(self, a, b):
[239] Fix | Delete
"""Return random integer in range [a, b], including both end points.
[240] Fix | Delete
"""
[241] Fix | Delete
[242] Fix | Delete
return self.randrange(a, b+1)
[243] Fix | Delete
[244] Fix | Delete
def _randbelow(self, n, _log=_log, _int=int, _maxwidth=1L<<BPF,
[245] Fix | Delete
_Method=_MethodType, _BuiltinMethod=_BuiltinMethodType):
[246] Fix | Delete
"""Return a random int in the range [0,n)
[247] Fix | Delete
[248] Fix | Delete
Handles the case where n has more bits than returned
[249] Fix | Delete
by a single call to the underlying generator.
[250] Fix | Delete
"""
[251] Fix | Delete
[252] Fix | Delete
try:
[253] Fix | Delete
getrandbits = self.getrandbits
[254] Fix | Delete
except AttributeError:
[255] Fix | Delete
pass
[256] Fix | Delete
else:
[257] Fix | Delete
# Only call self.getrandbits if the original random() builtin method
[258] Fix | Delete
# has not been overridden or if a new getrandbits() was supplied.
[259] Fix | Delete
# This assures that the two methods correspond.
[260] Fix | Delete
if type(self.random) is _BuiltinMethod or type(getrandbits) is _Method:
[261] Fix | Delete
k = _int(1.00001 + _log(n-1, 2.0)) # 2**k > n-1 > 2**(k-2)
[262] Fix | Delete
r = getrandbits(k)
[263] Fix | Delete
while r >= n:
[264] Fix | Delete
r = getrandbits(k)
[265] Fix | Delete
return r
[266] Fix | Delete
if n >= _maxwidth:
[267] Fix | Delete
_warn("Underlying random() generator does not supply \n"
[268] Fix | Delete
"enough bits to choose from a population range this large")
[269] Fix | Delete
return _int(self.random() * n)
[270] Fix | Delete
[271] Fix | Delete
## -------------------- sequence methods -------------------
[272] Fix | Delete
[273] Fix | Delete
def choice(self, seq):
[274] Fix | Delete
"""Choose a random element from a non-empty sequence."""
[275] Fix | Delete
return seq[int(self.random() * len(seq))] # raises IndexError if seq is empty
[276] Fix | Delete
[277] Fix | Delete
def shuffle(self, x, random=None):
[278] Fix | Delete
"""x, random=random.random -> shuffle list x in place; return None.
[279] Fix | Delete
[280] Fix | Delete
Optional arg random is a 0-argument function returning a random
[281] Fix | Delete
float in [0.0, 1.0); by default, the standard random.random.
[282] Fix | Delete
[283] Fix | Delete
"""
[284] Fix | Delete
[285] Fix | Delete
if random is None:
[286] Fix | Delete
random = self.random
[287] Fix | Delete
_int = int
[288] Fix | Delete
for i in reversed(xrange(1, len(x))):
[289] Fix | Delete
# pick an element in x[:i+1] with which to exchange x[i]
[290] Fix | Delete
j = _int(random() * (i+1))
[291] Fix | Delete
x[i], x[j] = x[j], x[i]
[292] Fix | Delete
[293] Fix | Delete
def sample(self, population, k):
[294] Fix | Delete
"""Chooses k unique random elements from a population sequence.
[295] Fix | Delete
[296] Fix | Delete
Returns a new list containing elements from the population while
[297] Fix | Delete
leaving the original population unchanged. The resulting list is
[298] Fix | Delete
in selection order so that all sub-slices will also be valid random
[299] Fix | Delete
samples. This allows raffle winners (the sample) to be partitioned
[300] Fix | Delete
into grand prize and second place winners (the subslices).
[301] Fix | Delete
[302] Fix | Delete
Members of the population need not be hashable or unique. If the
[303] Fix | Delete
population contains repeats, then each occurrence is a possible
[304] Fix | Delete
selection in the sample.
[305] Fix | Delete
[306] Fix | Delete
To choose a sample in a range of integers, use xrange as an argument.
[307] Fix | Delete
This is especially fast and space efficient for sampling from a
[308] Fix | Delete
large population: sample(xrange(10000000), 60)
[309] Fix | Delete
"""
[310] Fix | Delete
[311] Fix | Delete
# Sampling without replacement entails tracking either potential
[312] Fix | Delete
# selections (the pool) in a list or previous selections in a set.
[313] Fix | Delete
[314] Fix | Delete
# When the number of selections is small compared to the
[315] Fix | Delete
# population, then tracking selections is efficient, requiring
[316] Fix | Delete
# only a small set and an occasional reselection. For
[317] Fix | Delete
# a larger number of selections, the pool tracking method is
[318] Fix | Delete
# preferred since the list takes less space than the
[319] Fix | Delete
# set and it doesn't suffer from frequent reselections.
[320] Fix | Delete
[321] Fix | Delete
n = len(population)
[322] Fix | Delete
if not 0 <= k <= n:
[323] Fix | Delete
raise ValueError("sample larger than population")
[324] Fix | Delete
random = self.random
[325] Fix | Delete
_int = int
[326] Fix | Delete
result = [None] * k
[327] Fix | Delete
setsize = 21 # size of a small set minus size of an empty list
[328] Fix | Delete
if k > 5:
[329] Fix | Delete
setsize += 4 ** _ceil(_log(k * 3, 4)) # table size for big sets
[330] Fix | Delete
if n <= setsize or hasattr(population, "keys"):
[331] Fix | Delete
# An n-length list is smaller than a k-length set, or this is a
[332] Fix | Delete
# mapping type so the other algorithm wouldn't work.
[333] Fix | Delete
pool = list(population)
[334] Fix | Delete
for i in xrange(k): # invariant: non-selected at [0,n-i)
[335] Fix | Delete
j = _int(random() * (n-i))
[336] Fix | Delete
result[i] = pool[j]
[337] Fix | Delete
pool[j] = pool[n-i-1] # move non-selected item into vacancy
[338] Fix | Delete
else:
[339] Fix | Delete
try:
[340] Fix | Delete
selected = set()
[341] Fix | Delete
selected_add = selected.add
[342] Fix | Delete
for i in xrange(k):
[343] Fix | Delete
j = _int(random() * n)
[344] Fix | Delete
while j in selected:
[345] Fix | Delete
j = _int(random() * n)
[346] Fix | Delete
selected_add(j)
[347] Fix | Delete
result[i] = population[j]
[348] Fix | Delete
except (TypeError, KeyError): # handle (at least) sets
[349] Fix | Delete
if isinstance(population, list):
[350] Fix | Delete
raise
[351] Fix | Delete
return self.sample(tuple(population), k)
[352] Fix | Delete
return result
[353] Fix | Delete
[354] Fix | Delete
## -------------------- real-valued distributions -------------------
[355] Fix | Delete
[356] Fix | Delete
## -------------------- uniform distribution -------------------
[357] Fix | Delete
[358] Fix | Delete
def uniform(self, a, b):
[359] Fix | Delete
"Get a random number in the range [a, b) or [a, b] depending on rounding."
[360] Fix | Delete
return a + (b-a) * self.random()
[361] Fix | Delete
[362] Fix | Delete
## -------------------- triangular --------------------
[363] Fix | Delete
[364] Fix | Delete
def triangular(self, low=0.0, high=1.0, mode=None):
[365] Fix | Delete
"""Triangular distribution.
[366] Fix | Delete
[367] Fix | Delete
Continuous distribution bounded by given lower and upper limits,
[368] Fix | Delete
and having a given mode value in-between.
[369] Fix | Delete
[370] Fix | Delete
http://en.wikipedia.org/wiki/Triangular_distribution
[371] Fix | Delete
[372] Fix | Delete
"""
[373] Fix | Delete
u = self.random()
[374] Fix | Delete
try:
[375] Fix | Delete
c = 0.5 if mode is None else (mode - low) / (high - low)
[376] Fix | Delete
except ZeroDivisionError:
[377] Fix | Delete
return low
[378] Fix | Delete
if u > c:
[379] Fix | Delete
u = 1.0 - u
[380] Fix | Delete
c = 1.0 - c
[381] Fix | Delete
low, high = high, low
[382] Fix | Delete
return low + (high - low) * (u * c) ** 0.5
[383] Fix | Delete
[384] Fix | Delete
## -------------------- normal distribution --------------------
[385] Fix | Delete
[386] Fix | Delete
def normalvariate(self, mu, sigma):
[387] Fix | Delete
"""Normal distribution.
[388] Fix | Delete
[389] Fix | Delete
mu is the mean, and sigma is the standard deviation.
[390] Fix | Delete
[391] Fix | Delete
"""
[392] Fix | Delete
# mu = mean, sigma = standard deviation
[393] Fix | Delete
[394] Fix | Delete
# Uses Kinderman and Monahan method. Reference: Kinderman,
[395] Fix | Delete
# A.J. and Monahan, J.F., "Computer generation of random
[396] Fix | Delete
# variables using the ratio of uniform deviates", ACM Trans
[397] Fix | Delete
# Math Software, 3, (1977), pp257-260.
[398] Fix | Delete
[399] Fix | Delete
random = self.random
[400] Fix | Delete
while 1:
[401] Fix | Delete
u1 = random()
[402] Fix | Delete
u2 = 1.0 - random()
[403] Fix | Delete
z = NV_MAGICCONST*(u1-0.5)/u2
[404] Fix | Delete
zz = z*z/4.0
[405] Fix | Delete
if zz <= -_log(u2):
[406] Fix | Delete
break
[407] Fix | Delete
return mu + z*sigma
[408] Fix | Delete
[409] Fix | Delete
## -------------------- lognormal distribution --------------------
[410] Fix | Delete
[411] Fix | Delete
def lognormvariate(self, mu, sigma):
[412] Fix | Delete
"""Log normal distribution.
[413] Fix | Delete
[414] Fix | Delete
If you take the natural logarithm of this distribution, you'll get a
[415] Fix | Delete
normal distribution with mean mu and standard deviation sigma.
[416] Fix | Delete
mu can have any value, and sigma must be greater than zero.
[417] Fix | Delete
[418] Fix | Delete
"""
[419] Fix | Delete
return _exp(self.normalvariate(mu, sigma))
[420] Fix | Delete
[421] Fix | Delete
## -------------------- exponential distribution --------------------
[422] Fix | Delete
[423] Fix | Delete
def expovariate(self, lambd):
[424] Fix | Delete
"""Exponential distribution.
[425] Fix | Delete
[426] Fix | Delete
lambd is 1.0 divided by the desired mean. It should be
[427] Fix | Delete
nonzero. (The parameter would be called "lambda", but that is
[428] Fix | Delete
a reserved word in Python.) Returned values range from 0 to
[429] Fix | Delete
positive infinity if lambd is positive, and from negative
[430] Fix | Delete
infinity to 0 if lambd is negative.
[431] Fix | Delete
[432] Fix | Delete
"""
[433] Fix | Delete
# lambd: rate lambd = 1/mean
[434] Fix | Delete
# ('lambda' is a Python reserved word)
[435] Fix | Delete
[436] Fix | Delete
# we use 1-random() instead of random() to preclude the
[437] Fix | Delete
# possibility of taking the log of zero.
[438] Fix | Delete
return -_log(1.0 - self.random())/lambd
[439] Fix | Delete
[440] Fix | Delete
## -------------------- von Mises distribution --------------------
[441] Fix | Delete
[442] Fix | Delete
def vonmisesvariate(self, mu, kappa):
[443] Fix | Delete
"""Circular data distribution.
[444] Fix | Delete
[445] Fix | Delete
mu is the mean angle, expressed in radians between 0 and 2*pi, and
[446] Fix | Delete
kappa is the concentration parameter, which must be greater than or
[447] Fix | Delete
equal to zero. If kappa is equal to zero, this distribution reduces
[448] Fix | Delete
to a uniform random angle over the range 0 to 2*pi.
[449] Fix | Delete
[450] Fix | Delete
"""
[451] Fix | Delete
# mu: mean angle (in radians between 0 and 2*pi)
[452] Fix | Delete
# kappa: concentration parameter kappa (>= 0)
[453] Fix | Delete
# if kappa = 0 generate uniform random angle
[454] Fix | Delete
[455] Fix | Delete
# Based upon an algorithm published in: Fisher, N.I.,
[456] Fix | Delete
# "Statistical Analysis of Circular Data", Cambridge
[457] Fix | Delete
# University Press, 1993.
[458] Fix | Delete
[459] Fix | Delete
# Thanks to Magnus Kessler for a correction to the
[460] Fix | Delete
# implementation of step 4.
[461] Fix | Delete
[462] Fix | Delete
random = self.random
[463] Fix | Delete
if kappa <= 1e-6:
[464] Fix | Delete
return TWOPI * random()
[465] Fix | Delete
[466] Fix | Delete
s = 0.5 / kappa
[467] Fix | Delete
r = s + _sqrt(1.0 + s * s)
[468] Fix | Delete
[469] Fix | Delete
while 1:
[470] Fix | Delete
u1 = random()
[471] Fix | Delete
z = _cos(_pi * u1)
[472] Fix | Delete
[473] Fix | Delete
d = z / (r + z)
[474] Fix | Delete
u2 = random()
[475] Fix | Delete
if u2 < 1.0 - d * d or u2 <= (1.0 - d) * _exp(d):
[476] Fix | Delete
break
[477] Fix | Delete
[478] Fix | Delete
q = 1.0 / r
[479] Fix | Delete
f = (q + z) / (1.0 + q * z)
[480] Fix | Delete
u3 = random()
[481] Fix | Delete
if u3 > 0.5:
[482] Fix | Delete
theta = (mu + _acos(f)) % TWOPI
[483] Fix | Delete
else:
[484] Fix | Delete
theta = (mu - _acos(f)) % TWOPI
[485] Fix | Delete
[486] Fix | Delete
return theta
[487] Fix | Delete
[488] Fix | Delete
## -------------------- gamma distribution --------------------
[489] Fix | Delete
[490] Fix | Delete
def gammavariate(self, alpha, beta):
[491] Fix | Delete
"""Gamma distribution. Not the gamma function!
[492] Fix | Delete
[493] Fix | Delete
Conditions on the parameters are alpha > 0 and beta > 0.
[494] Fix | Delete
[495] Fix | Delete
The probability distribution function is:
[496] Fix | Delete
[497] Fix | Delete
x ** (alpha - 1) * math.exp(-x / beta)
[498] Fix | Delete
pdf(x) = --------------------------------------
[499] Fix | Delete
12
It is recommended that you Edit text format, this type of Fix handles quite a lot in one request
Function