| | |
- abcrand.ABCRand()
-
- InverseRandomStream
class InverseRandomStream(abcrand.ABCRand) |
| |
GENERAL CLASS FOR INITIATING RANDOM NUMBER STREAMS AND CREATING RANDOM
VARIATES FROM DIFFERENT PROBABILITY DISTRIBUTIONS
NB. All methods return a single random number on each call.
The class inherits the methods explicitly available from the abstract base
class ABCRand. The methods of Python's built-in Random class inherited via
ABCRand can not be guaranteed to be inverse-based. They are killed off
here and NOT available in InverseRandomStream! Please consult the docstrings
of the methods of the ABCRand class for information on how to use the
methods that ARE inherited.
The class is normally used like this:
rstream = InverseRandomStream()
meanx = 2.5
x = rstream.rexpo(meanx)
muy = -1.5
sigmay = 3.0
y = rstream.rnormal(muy, sigmay) # or whatever....
If another seed than the default seed is desired, then the class can be
instantiated using rstream = InverseRandomStream(nseed) where nseed is
an integer (preferably a large one).
A generating method of this class is in many cases much slower than the
corresponding method provided by the GeneralRandomStream class. The Latin
Hypercube and antithetic variance reduction methods provided by the
RandomStructure class only work with the methods in InverseRandomStream,
on the other hand. So speed is based on the generating methods as well as
the possibility of using variance reduction. Tests should be made when
speed is crucial.
InverseRandomStream makes it possible to use the RandomStructure class
to enter rank correlations between random parameters regardless of
distribution type, as well as correlated multinormal distributions.
It might be desirable at times to bound the output from a variate
generator for practical or physical reasons, for instance (some output
values might be unrealistic). A majority of the methods in this class
allow that. The actual bounds can be given (xmin and xmax), or the
corresponding cdf values of the unbound distribution can be specified
(pmin and pmax), a faster alternative since the method will otherwise
compute the corresponding cdf values every time the method is used.
Prior cdf value bounds may be computed by calling the corresponding
function in the statlib.cdf module.
BUT:
prescribing limits alters the distribution so the outputs will not truly
adhere to the theoretical one. The input parameters (mean, standard
deviation etc) are, anyhow, for the corresponding full-span theoretical
distribution!
NB The reason for using a generator not based on inversion is that a
particular generator not using inversion is faster. Bounding the range
of the output variate may alter this however, particularly for tight
bounds, so "ten mucho cuidado!". AND prescribing limits alters the
distribution so the outputs will not truly adhere to the theoretical
one. The input parameters (mean, std dev etc) are, however, for the
corresponding full-span theoretical distribution!
An externally generated input [0.0, 1.0] stream or "feed" can also be
used instead of letting the object instance pick the random variates
on its own. Cf. the docstrings of the ABCRand class for details!
NB. Some methods may return float('inf') or float('-inf') !
|
| |
- Method resolution order:
- InverseRandomStream
- abcrand.ABCRand
Methods defined here:
- __init__(self, nseed=2147483647)
- The seed 'nseed' must be a positive integer or a feed (a list or
a tuple) of numbers in [0.0, 1.0]!
- rNexpo(self, means, xmax=inf, pmax=1.0)
- Generator of random variates from a distribution of the sum of
exponentially distributed random variables. means[k] = 1.0/lambda[k].
NB Numbers in the means vector can be exactly equal! This generator
is slow, however...
- rNexpo2(self, means, xmax=inf, pmax=1.0)
- Generator of random variates from a distribution of the sum of
exponentially distributed random variables. means[k] = 1.0/lambda[k].
NB No two numbers in the means vector can be exactly equal! If this
is desired, use rNexpo !
- rbeta(self, a, b, x1, x2, betaab=False)
- The beta distribution f = x**(a-1) * (1-x)**(b-1) / beta(a, b)
The cdf is the integral = the incomplete beta or the incomplete
beta/complete beta depending on how the incomplete beta function
is defined.
x, a, b >= 0; x2 > x1
NB It is possible to provide the value of the complete beta
function beta(a, b) as a pre-computed input (may be computed
using numlib.specfunc.beta) instead of the default "False", a
feature that will make rbeta 30 % faster!
- rbootstrap(self, values)
- Picks elements from the input sequence (list, tuple etc) at random
(could be any sequence). The input sequence is sorted to assure
inverse properties.
- rconst(self, const)
- Returns the input constant. Use it to keep up the synchronization
even when a parameter has no spread (a dummy random number is sampled
each time the method is called)!
- rcoxian(self, means, probs, xmax=inf, pmax=1.0)
- Generates a random number from the Coxian phased distribution, which is
based on the exponential.
probs is a list of probabilities for GOING ON TO THE NEXT PHASE rather
than reaching the absorbing state prematurely. The number of means must
(of course) be one more than the number of probabilities!
NB means are allowed to be equal!
NB It is better to use rNexpo when all probs=1.0 !
- rcoxian2(self, means, probs, xmax=inf, pmax=1.0)
- Generates a random number from the Coxian phased distribution, which is
based on the exponential.
probs is a list of probabilities for GOING ON TO THE NEXT PHASE rather
than
reaching the absorbing state prematurely. The number of means must (of
course) be one more than the number of probabilities!
NB No two means must be equal - if this is desired, use rcoxian instead!
NB It is better to use rNexpo2 when all probs=1.0 !
- rdiscrete(self, values, qumul)
- Generates random variates from a user-defined discrete cdf.
'values' is a list/tuple with numbers, and 'qumul' are the corresponding
CUMULATIVE FREQUENCIES such that qumul[k] = P(x<=values[k]). The number
of values must be equal to the number of cumulative frequencies.
The cumulative frequencies must of course obey qumul[k+1] >= qumul[k],
otherwise an exception will be raised!
The 'values' sequence is sorted to assure inverse properties.
- rerlang(self, nshape, phasemean, xmax=inf, pmax=1.0)
- Generator of Erlang-distributed random variates.
Represents the sum of nshape exponentially distributed random variables,
each having the same mean value = phasemean. For nshape = 1 it works as
a generator of exponentially distributed random numbers.
- rerlang_gen(self, nshapes, qumul, phasemean, xmax=inf, pmax=1.0)
- The generalized Erlang distribution - the Erlang equivalent of the
rhyperexpo generator
f = sumk pk * ferlang(m, nk), F = sumk pk * Ferlang(m, nk), the same
mean for all phases.
NB Input to the function is the list of CUMULATIVE FREQUENCIES !
- rexpo(self, mean, xmax=inf, pmax=1.0)
- Generator of exponentially distributed random variates with
mean = 1.0/lambda
F = 1 - exp(-x/mean)
mean >= 0.0
- rexppower(self, loc, scale, alpha, lngam1oalpha=False, xmin=-inf, xmax=inf, pmin=0.0, pmax=1.0, tolf=8.8817841970012523e-16, itmax=128)
- The exponential power distribution
f = (a/s) * exp(-abs([x-l]/s)**a) / [2*gamma(1/a)]
F = 1/2 * [1 + sgn(x-l) * Fgamma(1/a, abs([x-l]/s)**a)], x in R
s, a > 0
where Fgamma is the gamma distribution cdf.
NB It is possible to gain efficiency by providing the value of the
natural logarithm of the complete gamma function ln(gamma(1.0/alpha))
as a pre-computed input (may be computed using numlib.specfunc.lngamma)
instead of the default 'False'.
tolf and itmax are the numerical control parameters of iexppower
and cexppower.
- rfoldednormal(self, muunfold, sigmaunfold, xmax=inf, pmax=1.0)
- The distribution of a random variable that is the absolute value
of a variate drawn from the normal distribution (i. e. the distribution
of a variate that is the absolute value of a normal variate, the latter
having muunfold as its mean and sigmaunfold as its standard deviation).
sigmaunfold >= 0.0
- rgamma(self, alpha, lam, lngamalpha=False, xmax=inf, pmax=1.0, tolf=8.8817841970012523e-16, itmax=128)
- The gamma distribution
f = lam * exp(-lam*x) * (lam*x)**(alpha-1) / gamma(alpha)
The cdf is the integral = the incomplete gamma ratio.
x, lam, alpha >= 0
NB It is possible to provide the value of the natural logarithm of the
complete gamma function lngamma(alpha) as a pre-computed input (may be
computed using numlib.specfunc.lngamma) instead of the default "False",
a feature that will make rgamma 50 % faster!
tolf and itmax are the numerical control parameters of igamma
and cgamma.
- rhyperexpo(self, means, qumul, xmax=inf, pmax=1.0)
- Generates a random number from the hyperexponential distribution
f = sumk pk * exp(x/mk) / mk, F = sumk pk * (1-exp(x/mk))
NB Input to the function is the list of CUMULATIVE FREQUENCIES !
- rlevy(self, scale, xmax=inf, pmax=1.0)
- The Levy distribution: f = sqrt(s/2pi) * (1/x)**(3/2) * exp(-s/2x)
F = erfc(sqrt(s/2x)); x >= 0
(stable distribution with alpha = 1/2, and beta = 1, aka the Cournot
distribution or the right-skewed Levy).
- rlognormal(self, mulg, sigmalg, xmax=inf, pmax=1.0)
- Generator of lognormally distributed random variates, inverse variant.
The log10-converted form is assumed for mulg and sigmalg:
mulg is the mean of the log10 (and the log10 of the median) of
the random variate, NOT the log10 of the mean of the non-logged
variate!, and sigmalg is the standard deviation of the log10 of
the random variate, NOT the log10 of the standard deviation of
the non-logged variate!!
sigmalg > 0.0
- rnormal(self, mu, sigma, xmin=-inf, xmax=inf, pmin=0.0, pmax=1.0)
- Generator of normally distributed random variates using an inverse
method, claimed to give a maximum relative error less than 2.6e-9
(cf. statlib.invcdf.inormal for details).
- rpareto(self, lam, xm, xmax=inf, pmax=1.0)
- The Pareto distribution:
f = lam * xm**lam / x**(lam+1)
F = 1 - (xm/x)**lam
x >= xm, lam > 0
For lam < 1 all moments are infinite
For lam < 2 all moments are infinite except for the mean
- rpoisson(self, lam, tspan, nmax=False, pmax=1.0)
- The Poisson distribution: p(N=n) = exp(-lam*tspan) * (lam*tspan)**n / n!
n = 0, 1, ...., infinity
A maximum number for the output may be given in nmax - then it must be
a positive integer.
- rstable_sym(self, alpha, location, scale, xmin=-inf, xmax=inf, pmin=0.0, pmax=1.0)
- The SYMMETRICAL stable distribution where alpha is the tail exponent.
For numerical reasons alpha is restricted to [0.25, 0.9] and
[1.125, 1.9] - but alpha = 1.0 (the Cauchy) and alpha = 2.0 (scaled
normal) are also allowed!
Numerics are somewhat crude but the fractional error is mostly < 0.001 -
sometimes much less - and the absolute error is almost always < 0.001 -
sometimes much less...
NB This generator is slow - particularly for small alpha !!!!!
- ruser_defined(self, ifunc, *args, pmin=0.0, pmax=1.0)
- Random deviate generation based on a user-defined inverse cdf placed in
a function ('ifunc') with arguments p in [0.0, 1.0] and sequence *args.
NB 1 *args is optional in that it is OK for ifunc to have no arguments.
NB 2 THE VALUES OF pmin AND pmax MUST - AND CAN ONLY - BE PASSED
BY KEYWORD: pmin=??, pmax=?? MUST ALWAYS BE WRITTEN OUT EXPLICITLY
UNLESS DEFAULTS ARE USED ('pmin' may be entered without explicit
keywording of pmax if the default is accepted for pmax).
- rweibull(self, c, scale, xmax=inf, pmax=1.0)
- Generator of random variates from the Weibull distribution.
F = 1 - exp[-(x/s)**c]
s >= 0.0, c >= 1.0
For c = 1.0 it is equivalent to a generator of exponential variates
with mean = scale
- rwiener(self, tau, xmin=-inf, xmax=inf, pmin=0.0, pmax=1.0)
- Generates random numbers corrsponding to a Wiener process (the
intergral of white noise (Langevin's function)), inverse variant.
The Wiener process is W(t+tau) - W (t) = N(0, sqrt(tau))
where tau is the time increment and N(0, sqrt(tau)) is a
normally distributed random variable having zero mean and
sqrt(tau) as its standard deviation.
This method returns W(t+tau) - W(t) given tau and allows
tau to be negative.
Data and other attributes defined here:
- __abstractmethods__ = frozenset([])
Methods inherited from abcrand.ABCRand:
- rcauchy(self, location, scale, xmin=-inf, xmax=inf, pmin=0.0, pmax=1.0)
- Generator of random variates from the Cauchy distribution:
f = 1 / [s*pi*(1 + [(x-l)/s]**2)]
F = (1/pi)*arctan((x-l)/s) + 1/2
(also known as the Lorentzian or Lorentz distribution)
scale must be >= 0
- rchistogram(self, values, qumul)
- Generates random variates from an input CUMULATIVE histogram.
'values' is a list/tuple with FLOATS in ascending order - A MUST!
These values represent bin end points and must be one more than
the number of cumulative frequencies, and where...
...'qumul' are the corresponding CUMULATIVE FREQUENCIES such that
qumul[k] = P(x<=values[k+1]).
The cumulative frequencies must of course obey qumul[k+1] >= qumul[k],
otherwise an exception will be raised!
The values of the random variate are assumed to be uniformly
distributed within each bin.
- rchistogram_int(self, values, qumul)
- Generates random variates from an input CUMULATIVE histogram.
'values' is a list/tuple with INTEGERS in ascending order - A MUST!
These values represent bin end points and must be one more than
the number of cumulative frequencies, and where...
...'qumul' are the corresponding CUMULATIVE FREQUENCIES such that
qumul[k] = P(x<=values[k+1]).
NB The first element of the values list is will never be returned!
The first integer to be returned is values[0] + 1 !!!!
The cumulative frequencies must of course obey qumul[k+1] >= qumul[k],
otherwise an exception will be raised!
The integer values of the random variate are assumed to be uniformly
distributed within each bin.
- remp_exp(self, values, npexp=0, ordered=False, xmax=inf, pmax=1.0)
- The mixed expirical/exponential distribution from Bratley, Fox and
Schrage. A polygon (piecewise linearly interpolated cdf) is used
together with a (shifted) exponential for the tail. The procedure
is designed so as to preserve the mean of the input sample.
The input is a set of observed points (vector) and an integer
representing the npexp largest points that will be used to formulate
the exponential tail.
NB it is assumed that x is in [0, inf) (with the usual cutoff
provisions) !!!!!
The function may also be used for a piecewise linear cdf without the
exponential tail - corrections are made to preserve the mean in this
case as well !!!
- rexpo_gen(self, a, b, c, xmin=-inf, xmax=inf, pmin=0.0, pmax=1.0)
- The generalized continuous double-sided exponential
distribution (x in R):
x <= c: f = [a*b/(a+b)] * exp(+a*[x-c])
F = [b/(a+b)] * exp(+a*[x-c])
x >= c: f = [a*b/(a+b)] * exp(-b*[x-c])
F = 1 - [a/(a+b)]*exp(-b*[x-c])
a > 0, b > 0
NB The symmetrical double-sided exponential sits in rlaplace!
- rextreme_I(self, type, mu, scale, xmin=-inf, xmax=inf, pmin=0.0, pmax=1.0)
- Extreme value distribution type I (aka the Gumbel distribution or
Gumbel distribution type I):
F = exp{-exp[-(x-mu)/scale]} (max variant)
f = exp[-(x-mu)/scale] * exp{-exp[-(x-mu)/scale]} / scale
F = 1 - exp{-exp[+(x-mu)/scale]} (min variant)
f = exp[+(x-mu)/scale] * exp{-exp[+(x-mu)/scale]} / scale
type must be 'max' or 'min'
scale must be > 0.0
- rextreme_gen(self, type, shape, mu, scale, xmin=-inf, xmax=inf, pmin=0.0, pmax=1.0)
- Generalized extreme value distribution:
F = exp{-[1-shape*(x-mu)/scale]**(1/shape)} (max version)
f = [1-shape*(x-mu)/scale]**(1/shape-1) *
exp{-[1-shape*(x-mu)/scale]**(1/shape)} / scale
F = 1 - exp{-[1+shape*(x-mu)/scale]**(1/shape)} (min version)
f = [1+shape*(x-mu)/scale]**(1/shape-1) *
exp{-[1+shape*(x-mu)/scale]**(1/shape)} / scale
shape < 0 => Type II
shape > 0 => Type III
shape -> 0 => Type I - Gumbel
type must be 'max' or 'min'
scale must be > 0.0
A REASONABLE SCHEME SEEMS TO BE mu = scale WHICH SEEMS TO LIMIT THE
DISTRIBUTION TO EITHER SIDE OF THE Y-AXIS!
- rgeometric(self, phi)
- The geometric distribution with p(K=k) = phi * (1-phi)**(k-1) and
P(K>=k) = sum phi * (1-phi)**k = 1 - q**k, where q = 1 - phi and
0 < phi <= 1 is the success frequency or "Bernoulli probability"
and K >= 1 is the number of trials to the first success in a series
of Bernoulli trials. It is easy to prove that P(k) = 1 - (1-phi)**k:
let q = 1 - phi. p(k) = (1-q) * q**(k-1) = q**(k-1) - q**k. Then P(1) =
p(1) = 1 - q. P(2) = p(1) + p(2) = 1 - q + q - q**2 = 1 - q**2.
Induction can be used to show that P(k) = 1 - q**k = 1 - (1-phi)**k
- rkodlin(self, gam, eta, xmax=inf, pmax=1.0)
- The Kodlin distribution, aka the linear hazard rate distribution:
f = (gam + eta*x) * exp{-[gam*x + (1/2)*eta*x**2]},
F = 1 - exp{-[gam*x + (1/2)*eta*x**2]}; x, gam, eta >= 0
- rkumaraswamy(self, a, b, x1, x2)
- The Kumaraswamy distribution f = a*b*x**(a-1) * (1-x**a)**(b-1)
F = 1 - (1-x**a)**b
a, b >= 0; 0 <= x <= 1
The Kumaraswamy is very similar to the beta distribution !!!
x2 >= x1 !!!!
- rlaplace(self, loc, scale, xmin=-inf, xmax=inf, pmin=0.0, pmax=1.0)
- The Laplace aka the symmetrical double-sided exponential distribution
f = ((1/2)/s)) * exp(-abs([x-l]/s))
F = (1/2)*exp([x-l]/s) {x <= 0}, F = 1 - (1/2)*exp(-[x-l]/s) {x >= 0}
s >= 0
- rlogistic(self, mu, scale, xmin=-inf, xmax=inf, pmin=0.0, pmax=1.0)
- The logistic distribution: F = 1 / {1 + exp[-(x-m)/s]}; x on R
m is the mean and mode, and s is a scale parameter (s >= 0)
- rpareto_zero(self, lam, xm, xmax=inf, pmax=1.0)
- The Pareto distribution with the support shifted to [0, inf):
f = lam * xm**lam / (x+xm)**(lam+1)
F = 1 - [xm/(x+xm)]**lam
x in [0, inf)
lam > 0
For lam < 1 all moments are infinite
For lam < 2 all moments are infinite except for the mean
- rrayleigh(self, sigma, xmax=inf, pmax=1.0)
- The Rayleigh distribution:
f = (x/s**2) * exp[-x**2/(2*s**2)]
F = 1 - exp[-x**2/(2*s**2)]
x, s >= 0
- rsign(self)
- Returns -1.0 or 1.0 with probability 0.5 for each.
- rsinus(self, left, right)
- The "sinus distribution".
- rtri_unif_tri(self, a, b, c, d)
- Triangular-uniform-triangular distribution with support on [a, d] and
with breakpoints in b and c.
- rtriang(self, left, mode, right)
- Generator of triangularly distributed random numbers on [left, right]
with the peak of the pdf at mode.
- rtukeylambda_gen(self, lam1, lam2, lam3, lam4, pmin=0.0, pmax=1.0)
- The Friemer-Mudholkar-Kollia-Lin generalized Tukey lambda distribution.
lam1 is a location parameter and lam2 a scale parameter. lam3 and lam4
are associated with the shape of the distribution.
- runif_int0N(self, number)
- Generator of uniformly distributed integers in [0, n) (also the basis
of some other procedures for generating random variates).
Numbers returned are 0 through n-1.
- runifab(self, left, right)
- Generator of uniformly distributed floats between 'left' and 'right'.
Data descriptors inherited from abcrand.ABCRand:
- __dict__
- dictionary for instance variables (if defined)
- __weakref__
- list of weak references to the object (if defined)
| |