Random Deviates

GalSim can produce random values according to a variety of probability distributions:

  • UniformDeviate implements \(p(x) = 1\) for \(0 \le x < 1\).

  • GaussianDeviate implements \(p(x) = \frac{1}{\sqrt{2\pi\sigma}} e^{-\frac{(x-\mu)^2}{2\sigma^2}}\).

  • PoissonDeviate implements \(p(x) = \frac{e^{-\mu}\mu^x}{x!}\) for integer \(x > 0\).

  • BinomialDeviate implements \(p(x) = {N \choose x}p^k(1-p)^{N-x}\) for integer \(0 \le x \le N\).

  • Chi2Deviate implements \(p(x) = \frac{x^{(n/2)-1}e^{-x/2}}{\Gamma(n/2)2^{n/2}}\) for \(x > 0\).

  • GammaDeviate implements \(p(x) = x^{k-1}\frac{e^{-x/\theta}}{\theta^k\Gamma(k)}\) for \(x > 0\).

  • WeibullDeviate implements \(p(x) = \frac{a}{b}\left(\frac{x}{b}\right)^{a-1}e^{-\left(\frac{x}{b}\right)^a}\) for \(x \ge 0\).

  • DistDeviate implements any arbitrary, user-supplied \(p(x)\).

These are all subclasses of the base class BaseDeviate, which implements the underlying pseudo-random number generator using the Boost libraries Mersenne twister.

We have fixed the implementation of this to Boost version 1.48.0, the relevant files of which are bundled with the GalSim distribution, so that random numbers produced by GalSim simulations are deterministic across different user platforms and operating systems. These Boost files are included with GalSim, so the user does not need to have Boost installed on their system.

There are ways to connect various different deviate objects to use the same underlying BaseDeviate, which is often important for producing deterministic simulations given a particular random number seed. See the docstring of BaseDeviate for details.

Note

We have put some care into the way we seed the random number generator such that it is safe to start several random number sequences seeded by sequential seeds. This is already supposed to be the case for the Boost Mersenne Twister implementation, but we add some extra (probably overly paranoid) steps to ensure this by seeding one pseudo-rng, skip a few values, and then use that to seed the actual pseudo-rng that we will use.

This means you can start the rngs for sequential images or even galaxies with sequential seed values and there will not be any measurable correlations in the results. This can greatly ease the ability to split work across multiple processes and still achieve deterministic results.

class galsim.BaseDeviate(seed=None)[source]

Base class for all the various random deviates.

This holds the essential random number generator that all the other classes use.

All deviates take an initial seed argument that is used to seed the underlying random number generator. It has three different kinds of behavior.

  1. An integer value can be provided to explicitly seed the random number generator with a particular value. This is useful to have deterministic behavior. If you seed with an integer value, the subsequent series of “random” values will be the same each time you run the program.

  2. A seed of 0 or None means to pick some arbitrary value that will be different each time you run the program. Currently, this tries to get a seed from /dev/urandom if possible. If that doesn’t work, then it creates a seed from the current time. You can also get this behavior by omitting the seed argument entirely. (i.e. the default is None.)

  3. Providing another BaseDeviate object as the seed will make the new BaseDeviate share the same underlying random number generator as the other BaseDeviate. So you can make one BaseDeviate (of any type), and seed it with a particular deterministic value. Then if you pass that BaseDeviate to any other one you make, they will both be using the same RNG and the series of “random” values will be deterministic.

Usage:

The only kind of random number you can obtain from a pure BaseDeviate (i.e. one that is not actually one of the variosu subclasses) is a “raw” value. This is an unsigned 32-bit integer that behind the scenes is used by all sub-classes to generate floating point values with various distributions.

>>> rng = galsim.BaseDeviate(215324)
>>> rng.raw()
3559052779

Most other usage is effected by creating sub-classes corresponding to particular random deviates with various distributions. E.g. UniformDeviate generates random values following a uniform distribution between 0 and 1.

>>> rng = galsim.BaseDeviate(215324)
>>> ud = galsim.UniformDeviate(rng)
>>> ud()
0.58736140513792634
>>> ud2 = galsim.UniformDeviate(215324)
>>> ud2()
0.58736140513792634
_seed(seed=0)[source]

Equivalent to seed, but without any type checking.

_reset(rng)[source]

Equivalent to reset, but rng must be a BaseDeviate (not an int), and there is no type checking.

add_generate(array)[source]

Generate many pseudo-random values, adding them to the values of a numpy array.

as_numpy_generator()[source]

Return a numpy.random.Generator object that uses the current BaseDeviate for the underlying bit generations.

This allows you to use the (probably) more familiar numpy functions, while maintaining GalSim’s guarantees about random number stability across platforms.

Example:

>>> rng = galsim.BaseDeviate(1234)
>>> gen = rng.as_numpy_generator()
>>> uniform = gen.uniform(1, 10, size=10)
>>> norm = gen.normal(0, 3, size=20)

There is also a shorthand syntax that may be convenient. The property np is equivalent to this method, so you can also write:

>>> uniform = rng.np.uniform(1, 10, size=10)
>>> norm = rng.np.normal(0, 3, size=20)
clearCache()[source]

Clear the internal cache of the BaseDeviate, if any. This is currently only relevant for GaussianDeviate, since it generates two values at a time, saving the second one to use for the next output value.

discard(n, suppress_warnings=False)[source]

Discard n values from the current sequence of pseudo-random numbers.

This is typically used to keep two random number generators in sync when one of them is used to generate some random values. The other can quickly discard the same number of random values to stay in sync with the first one.

Parameters:

  • n – The number of raw random numbers to discard.

  • suppress_warnings – Whether to suppress warnings related to detecting when this action is not likely to reliably keep two random number generators in sync. [default: False]

duplicate()[source]

Create a duplicate of the current BaseDeviate object.

The subsequent series from each copy of the BaseDeviate will produce identical values:

>>> u = galsim.UniformDeviate(31415926)
>>> u()
0.17100770119577646
>>> u2 = u.duplicate()
>>> u()
0.49095047544687986
>>> u()
0.10306670609861612
>>> u2()
0.49095047544687986
>>> u2()
0.10306670609861612
>>> u2()
0.13129289541393518
>>> u()
0.13129289541393518
generate(array)[source]

Generate many pseudo-random values, filling in the values of a numpy array.

property generates_in_pairs

Indicates whether the generator uses 2 rngs values per 2 returned values.

GaussianDeviate has a slight wrinkle to the BaseDeviate.has_reliable_discard story. It always uses two rng values to generate two Gaussian deviates. So if the number of generated values is even, then discard can keep things in sync. However, if an odd number of values are generated, then you to generate one more value (and discard it) to keep things synced up.

This is only True for GaussianDeviate.

property has_reliable_discard

Indicates whether the generator always uses 1 rng per value.

If it does, then discard can reliably be used to keep two generators in sync when one of them generated some values and the other didn’t.

This is False for PoissonDeviate, Chi2Deviate, and GammaDeviate.

See also BaseDeviate.generates_in_pairs.

property np

Shorthand for self.as_numpy_generator()

raw()[source]

Generate the next pseudo-random number and rather than return the appropriate kind of random deviate for this class, just return the raw integer value that would have been used to generate this value.

reset(seed=None)[source]

Reset the pseudo-random number generator, severing connections to any other deviates. Providing another BaseDeviate object as the seed connects this deviate with the other one, so they will both use the same underlying random number generator.

Parameters:

seed – Something that can seed a BaseDeviate: an integer seed or another BaseDeviate. Using None means to generate a seed from the system. [default: None]

seed(seed=0)[source]

Seed the pseudo-random number generator with a given integer value.

Parameters:

seed – An int value to be used to seed the random number generator. Using 0 means to generate a seed from the system. [default: 0]

class galsim.UniformDeviate(seed=None)[source]

Bases: BaseDeviate

Pseudo-random number generator with uniform distribution in interval [0.,1.).

Successive calls to u() generate pseudo-random values distributed uniformly in the interval [0., 1.):

>>> u = galsim.UniformDeviate(31415926)
>>> u()
0.17100770119577646
>>> u()
0.49095047544687986
Parameters:

seed – Something that can seed a BaseDeviate: an integer seed or another BaseDeviate. Using 0 means to generate a seed from the system. [default: None]

__call__()[source]

Draw a new random number from the distribution.

Returns a uniform deviate between 0 and 1.

class galsim.GaussianDeviate(seed=None, mean=0.0, sigma=1.0)[source]

Bases: BaseDeviate

Pseudo-random number generator with Gaussian distribution.

See http://en.wikipedia.org/wiki/Gaussian_distribution for further details.

Successive calls to g() generate pseudo-random values distributed according to a Gaussian distribution with the provided mean, sigma:

>>> g = galsim.GaussianDeviate(31415926)
>>> g()
0.5533754000847082
>>> g()
1.0218588970190354
Parameters:

  • seed – Something that can seed a BaseDeviate: an integer seed or another BaseDeviate. Using 0 means to generate a seed from the system. [default: None]

  • mean – Mean of Gaussian distribution. [default: 0.]

  • sigma – Sigma of Gaussian distribution. [default: 1.; Must be > 0]

__call__()[source]

Draw a new random number from the distribution.

Returns a Gaussian deviate with the given mean and sigma.

generate_from_variance(array)[source]

Generate many Gaussian deviate values using the existing array values as the variance for each.

property generates_in_pairs

Indicates whether the generator uses 2 rngs values per 2 returned values.

GaussianDeviate has a slight wrinkle to the BaseDeviate.has_reliable_discard story. It always uses two rng values to generate two Gaussian deviates. So if the number of generated values is even, then discard can keep things in sync. However, if an odd number of values are generated, then you to generate one more value (and discard it) to keep things synced up.

This is only True for GaussianDeviate.

property mean

The mean of the Gaussian distribution.

property sigma

The sigma of the Gaussian distribution.

class galsim.PoissonDeviate(seed=None, mean=1.0)[source]

Bases: BaseDeviate

Pseudo-random Poisson deviate with specified mean.

The input mean sets the mean and variance of the Poisson deviate. An integer deviate with this distribution is returned after each call. See http://en.wikipedia.org/wiki/Poisson_distribution for more details.

Successive calls to p() generate pseudo-random integer values distributed according to a Poisson distribution with the specified mean:

>>> p = galsim.PoissonDeviate(31415926, mean=100)
>>> p()
94
>>> p()
106
Parameters:

  • seed – Something that can seed a BaseDeviate: an integer seed or another BaseDeviate. Using 0 means to generate a seed from the system. [default: None]

  • mean – Mean of the distribution. [default: 1; Must be > 0]

__call__()[source]

Draw a new random number from the distribution.

Returns a Poisson deviate with the given mean.

generate_from_expectation(array)[source]

Generate many Poisson deviate values using the existing array values as the expectation value (aka mean) for each.

property has_reliable_discard

Indicates whether the generator always uses 1 rng per value.

If it does, then discard can reliably be used to keep two generators in sync when one of them generated some values and the other didn’t.

This is False for PoissonDeviate, Chi2Deviate, and GammaDeviate.

See also BaseDeviate.generates_in_pairs.

property mean

The mean of the distribution.

class galsim.BinomialDeviate(seed=None, N=1, p=0.5)[source]

Bases: BaseDeviate

Pseudo-random Binomial deviate for N trials each of probability p.

N is number of ‘coin flips,’ p is probability of ‘heads,’ and each call returns an integer value where 0 <= value <= N gives the number of heads. See http://en.wikipedia.org/wiki/Binomial_distribution for more information.

Successive calls to b() generate pseudo-random integer values distributed according to a binomial distribution with the provided N, p:

>>> b = galsim.BinomialDeviate(31415926, N=10, p=0.3)
>>> b()
2
>>> b()
3
Parameters:

  • seed – Something that can seed a BaseDeviate: an integer seed or another BaseDeviate. Using 0 means to generate a seed from the system. [default: None]

  • N – The number of ‘coin flips’ per trial. [default: 1; Must be > 0]

  • p – The probability of success per coin flip. [default: 0.5; Must be > 0]

__call__()[source]

Draw a new random number from the distribution.

Returns a Binomial deviate with the given n and p.

property n

The number of ‘coin flips’.

property p

The probability of success per ‘coin flip’.

class galsim.Chi2Deviate(seed=None, n=1.0)[source]

Bases: BaseDeviate

Pseudo-random Chi^2-distributed deviate for degrees-of-freedom parameter n.

See http://en.wikipedia.org/wiki/Chi-squared_distribution (note that k=n in the notation adopted in the Boost.Random routine called by this class). The Chi^2 distribution is a real-valued distribution producing deviates >= 0.

Successive calls to chi2() generate pseudo-random values distributed according to a chi-square distribution with the specified degrees of freedom, n:

>>> chi2 = galsim.Chi2Deviate(31415926, n=7)
>>> chi2()
7.9182211987712385
>>> chi2()
6.644121724269535
Parameters:

  • seed – Something that can seed a BaseDeviate: an integer seed or another BaseDeviate. Using 0 means to generate a seed from the system. [default: None]

  • n – Number of degrees of freedom for the output distribution. [default: 1; Must be > 0]

__call__()[source]

Draw a new random number from the distribution.

Returns a Chi2-distributed deviate with the given number of degrees of freedom.

property has_reliable_discard

Indicates whether the generator always uses 1 rng per value.

If it does, then discard can reliably be used to keep two generators in sync when one of them generated some values and the other didn’t.

This is False for PoissonDeviate, Chi2Deviate, and GammaDeviate.

See also BaseDeviate.generates_in_pairs.

property n

The number of degrees of freedom.

class galsim.GammaDeviate(seed=None, k=1.0, theta=1.0)[source]

Bases: BaseDeviate

A Gamma-distributed deviate with shape parameter k and scale parameter theta. See http://en.wikipedia.org/wiki/Gamma_distribution. (Note: we use the k, theta notation. If you prefer alpha, beta, use k=alpha, theta=1/beta.) The Gamma distribution is a real valued distribution producing deviates >= 0.

Successive calls to g() generate pseudo-random values distributed according to a gamma distribution with the specified shape and scale parameters k and theta:

>>> gam = galsim.GammaDeviate(31415926, k=1, theta=2)
>>> gam()
0.37508882726316
>>> gam()
1.3504199388358704
Parameters:

  • seed – Something that can seed a BaseDeviate: an integer seed or another BaseDeviate. Using 0 means to generate a seed from the system. [default: None]

  • k – Shape parameter of the distribution. [default: 1; Must be > 0]

  • theta – Scale parameter of the distribution. [default: 1; Must be > 0]

__call__()[source]

Draw a new random number from the distribution.

Returns a Gamma-distributed deviate with the given k and theta.

property has_reliable_discard

Indicates whether the generator always uses 1 rng per value.

If it does, then discard can reliably be used to keep two generators in sync when one of them generated some values and the other didn’t.

This is False for PoissonDeviate, Chi2Deviate, and GammaDeviate.

See also BaseDeviate.generates_in_pairs.

property k

The shape parameter, k.

property theta

The scale parameter, theta.

class galsim.WeibullDeviate(seed=None, a=1.0, b=1.0)[source]

Bases: BaseDeviate

Pseudo-random Weibull-distributed deviate for shape parameter a and scale parameter b.

The Weibull distribution is related to a number of other probability distributions; in particular, it interpolates between the exponential distribution (a=1) and the Rayleigh distribution (a=2). See http://en.wikipedia.org/wiki/Weibull_distribution (a=k and b=lambda in the notation adopted in the Wikipedia article) for more details. The Weibull distribution is real valued and produces deviates >= 0.

Successive calls to w() generate pseudo-random values distributed according to a Weibull distribution with the specified shape and scale parameters a and b:

>>> w = galsim.WeibullDeviate(31415926, a=1.3, b=4)
>>> w()
1.1038481241018219
>>> w()
2.957052966368049
Parameters:

  • seed – Something that can seed a BaseDeviate: an integer seed or another BaseDeviate. Using 0 means to generate a seed from the system. [default: None]

  • a – Shape parameter of the distribution. [default: 1; Must be > 0]

  • b – Scale parameter of the distribution. [default: 1; Must be > 0]

__call__()[source]

Draw a new random number from the distribution.

Returns a Weibull-distributed deviate with the given shape parameters a and b.

property a

The shape parameter, a.

property b

The scale parameter, b.

class galsim.DistDeviate(seed=None, function=None, x_min=None, x_max=None, interpolant=None, npoints=None, clip_neg=False)[source]

Bases: BaseDeviate

A class to draw random numbers from a user-defined probability distribution.

DistDeviate is a BaseDeviate class that can be used to draw from an arbitrary probability distribution. The probability distribution passed to DistDeviate can be given one of three ways: as the name of a file containing a 2d ASCII array of x and P(x), as a LookupTable mapping x to P(x), or as a callable function.

Once given a probability, DistDeviate creates a table of the cumulative probability and draws from it using a UniformDeviate. The precision of its outputs can be controlled with the keyword npoints, which sets the number of points DistDeviate creates for its internal table of CDF(x). To prevent errors due to non-monotonicity, the interpolant for this internal table is always linear.

Two keywords, x_min and x_max, define the support of the function. They must be passed if a callable function is given to DistDeviate, unless the function is a LookupTable, which has its own defined endpoints. If a filename or LookupTable is passed to DistDeviate, the use of x_min or x_max will result in an error.

If given a table in a file, DistDeviate will construct an interpolated LookupTable to obtain more finely gridded probabilities for generating the cumulative probability table. The default interpolant is linear, but any interpolant understood by LookupTable may be used. We caution against the use of splines because they can cause non-monotonic behavior. Passing the interpolant keyword next to anything but a table in a file will result in an error.

Examples:

Some sample initialization calls:

>>> d = galsim.DistDeviate(function=f, x_min=x_min, x_max=x_max)

Initializes d to be a DistDeviate instance with a distribution given by the callable function f(x) from x=x_min to x=x_max and seeds the PRNG using current time:

>>> d = galsim.DistDeviate(1062533, function=file_name, interpolant='floor')

Initializes d to be a DistDeviate instance with a distribution given by the data in file file_name, which must be a 2-column ASCII table, and seeds the PRNG using the integer seed 1062533. It generates probabilities from file_name using the interpolant ‘floor’:

>>> d = galsim.DistDeviate(rng, function=galsim.LookupTable(x,p))

Initializes d to be a DistDeviate instance with a distribution given by P(x), defined as two arrays x and p which are used to make a callable LookupTable, and links the DistDeviate PRNG to the already-existing random number generator rng.

Successive calls to d() generate pseudo-random values with the given probability distribution:

>>> d = galsim.DistDeviate(31415926, function=lambda x: 1-abs(x), x_min=-1, x_max=1)
>>> d()
-0.4151921102709466
>>> d()
-0.00909781188974034
Parameters:

  • seed – Something that can seed a BaseDeviate: an integer seed or another BaseDeviate. Using 0 means to generate a seed from the system. [default: None]

  • function – A callable function giving a probability distribution or the name of a file containing a probability distribution as a 2-column ASCII table. [required]

  • x_min – The minimum desired return value (required for non-LookupTable callable functions; will raise an error if not passed in that case, or if passed in any other case) [default: None]

  • x_max – The maximum desired return value (required for non-LookupTable callable functions; will raise an error if not passed in that case, or if passed in any other case) [default: None]

  • interpolant – Type of interpolation used for interpolating a file (causes an error if passed alongside a callable function). Options are given in the documentation for LookupTable. [default: ‘linear’]

  • npoints – Number of points DistDeviate should create for its internal interpolation tables. [default: 256, unless the function is a non-log LookupTable, in which case it uses the table’s x values]

  • clip_neg – Clip any negative input values to zero. [default: False; an error will be raised if any negative probabilities are found.]

__call__()[source]

Draw a new random number from the distribution.

add_generate(array)[source]

Generate many pseudo-random values, adding them to the values of a numpy array.

generate(array)[source]

Generate many pseudo-random values, filling in the values of a numpy array.

val(p)[source]

Return the value \(x\) of the input function to DistDeviate such that p = \(F(x)\), where \(F\) is the cumulattive probability distribution function:

\[F(x) = \int_{-\infty}^x \mathrm{pdf}(t) dt\]

This function is typically called by __call__, which generates a random p between 0 and 1 and calls self.val(p).

Parameters:

p – The desired cumulative probabilty p.

Returns:

the corresponding x such that \(p = F(x)\).