A sequence of nonrecurring events can be simulated via RandomSample. For instance, the probability of randomly sampling the integers 1 through in order might be simulated.

This estimates the probability of getting elements in order for from 2 to 8:

Random number generation is at the heart of Monte Carlo estimates. An estimate of an expected value of a function can be obtained by generating values from the desired distribution and finding the mean of applied to those values.

This provides no evidence that and are different for real values:

This provides evidence that and differ for at least some complex values:

Random Numbers

RandomInteger[{imin,imax}]	give a pseudorandom integer in the range {imin,…,imax}
RandomInteger[imax]	give a pseudorandom integer in the range {0,…,imax}
RandomInteger[]	pseudorandomly give 0 or 1 with probability
RandomInteger[domain,n]	give a list of n pseudorandom integers
RandomInteger[domain,{n1,n2,…}]	give an n1×n2×… array of pseudorandom integers

For multidimensional arrays with dimensions through , the total number of required pseudorandom numbers is generated and then partitioned. This makes the multidimensional array generation as efficient as possible because the total number of random values is generated as efficiently as possible and the time required for partitioning is negligible.

The following converts the first million digits of to a string of integers:

This gives the positions where the string of five digits appears in the first million digits of :

The following example will construct a Gibbs sampler for an example given by Casella and George [2]. The distribution of interest is bivariate. The conditional distribution of given is a binomial, and the conditional distribution of given is a beta. As Casella and George mention, various strategies for detecting convergence and sampling using the Gibbs sampler have been suggested. For simplicity, assume that convergence will occur within 1000 iterations. A sample of size from the distribution will be taken as the values following the 1000 ^th iteration. It should be noted that these values will, however, be dependent.

This compares the empirical and theoretical distributions of for :

This compares the empirical and theoretical distributions of for :

Arbitrary-Precision Reals and Complexes

This gives a precision-50 ‐distributed number:

This estimates the worst precision loss in computing on the interval [0,1000]:

Random Elements

If the first argument to RandomChoice or RandomSample is a list, elements are selected with equal probability. The weight specification defines a distribution on the set of the e_i. The weights must be positive, but need not sum to 1. For weights {w₁,…,w_n} the probability of e_i in the initial distribution is . Since RandomSample samples without replacement, weights are updated internally based on the total remaining weight after each selection.

SeedRandom and BlockRandom in Parallel Computations

This defines a CompiledFunction that approximates the area of a quarter circle using samples that will run in parallel when given a list:

"Congruential"	linear congruential generator (low ‐ quality randomness)
"ExtendedCA"	extended cellular automaton generator (default)
"Legacy"	default generators prior to version 6.0
"MersenneTwister"	Mersenne Twister shift register generator
"MKL"	Intel MKL generator (Intel ‐ based systems)
"ParallelGenerator"	used for initializing and seeding generators for parallel computations.
"ParallelMersenneTwister"	set of 1024 Mersenne Twister generators of period
"Rule30CA"	Wolfram rule 30 generator

Congruential

"Congruential" uses a linear congruential generator. This is one of the simplest types of pseudorandom number generators, with pseudorandom numbers between 0 and 1 obtained from , where is given by the modular recurrence relation

for some fixed integers , , and called the multiplier, increment, and modulus, respectively. If the increment is 0, the generator is a multiplicative congruential generator. The values of , , and can be set via options to the "Congruential" method.

Linear congruential generators are periodic and tend to give a lower quality of randomness, especially when a large number of random values is needed. If reals are generated directly from the congruence relation, the period is less than or equal to .

The actual period can be determined by finding the smallest integer such that mod .

The "Bits" option can be Automatic, a nonzero integer, or a list of two nonzero integers specifying the range of bits in the modulus used for constructing numbers from bits. Automatic uses {2,-1} unless is a power of 2, in which case {1,-1} is used.

ExtendedCA

The length of the state vectors used is by default set to cells. The multiple of 64 can be controlled by the "Size" option. Once a state vector is computed by evolving the cellular automaton using the five-neighbor rule, bits are selected for random numbers from bits {start,start+skip,…}.

Legacy

MersenneTwister

"MersenneTwister" uses the Mersenne Twister generator due to Matsumoto and Nishimura [3][4]. The Mersenne Twister is a generalized feedback shift register generator with period .

Reference implementation

genrand64_int64 ()

#include <stdio.h>
#include "mt64.h"

int main(void) {
int i;
unsigned long long init[1], length=1;

init[0]=1ULL; /*SeedRandom[1, Method -> "MersenneTwister"];*/
init_by_array64(init, length);

for (i=0; i<5; i++) {
   printf("%20llu\n", genrand64_int64());
}
    
return 0;
}

/*Output*/ 
7259937129391483703
 7973299316636211948
16865006314979686608
 5442441613857606270
14480929463982189498

Wolfram Language

Reals

genrand64_real1()

#include <stdio.h>
#include "mt64.h"

int main(void) {
int i;
unsigned long long init[1], length=1;

init[0]=1ULL; /*SeedRandom[1, Method -> "MersenneTwister"];*/
init_by_array64(init, length);

for (i=0; i<3; i++) {
        printf("%16.15f\n", genrand64_real1());
    }
    
return 0;
}

/*Output*/ 
0.393561980389721
0.432233422048709
0.914253824284137

Wolfram Language

MKL

Rule30CA

where is the value of cell at time .

The length of the state vectors used is by default set to cells. The multiplier for 29 can be controlled by the "Size" option.

ParallelMersenneTwister

"ParallelMersenneTwister" uses a set of Mersenne Twister generators due to Matsumoto and Nishimura [3][4] with parameters chosen using their "Dynamic Creator" program dcmt [19]. The program computes parameters for the Mersenne Twister generator that are relatively prime and so should produce independent results. The parameters were computed to produce Mersenne Twister generalized feedback shift register generators with period .

ParallelGenerator

"ParallelMersenneTwister"	parallel Mersenne twister generators with period
"ExtendedCA"	extended CA generators with different starting positions
f	generator f[i] used for the i th thread
"Default"	restores the default method

Defining Your Own Generator

If random integers are supported, then gobj["GenerateIntegers"[n,{a,b}]] is expected to return a list of n random integers in the range . A warning message will be issued when results are out of range.

If random reals are supported, then gobj["GenerateReals"[n,{a,b},prec]] is expected to return a list of n random reals with precision prec in the range . A warning message will be issued when results are out of range or of the wrong precision.

Example: Multiplicative Congruential Generator

Example: Blum–Blum–Shub Generator

The Blum–Blum–Shub generator is a quadratic congruential method for generating pseudorandom bits for cryptographic purposes [5]. The congruence is mod × for specified primes and .

Continuous Distributions

NoncentralChiSquareDistribution[ν,λ], , variate generation uses additive properties of distributions to avoid expensive inverse CDF computations for nonintegral . The additive properties are given in, for instance, Johnson, Kotz, and Balakrishnan [7]. For a noncentral variate can be generated as the square of a normal variate with mean and variance 1. For noncentral variates are obtained as the sum of a central and a noncentral random variable. For , is distributed if and . This relationship cannot be used for . In that case the construction is with and , where is the limiting noncentral distribution as goes to 0. The limiting distribution is a mixture of Poisson and variables, which has a nonzero probability mass at 0 and a continuous density for positive values. NoncentralFRatioDistribution variates are obtained from one central and one noncentral variate.

For GammaDistribution[α,β] exponential variates are generated when . Otherwise, methods due to Cheng and Feast [9] and Ahrens and Dieter [10] are used.

Beta variates are constructed by switching between multiple methods, depending on the values of the beta parameters and . If both parameters are 1, uniform random variates will be generated. If one of the beta parameters is 1, then a closed-form inverse CDF evaluation is used. Otherwise, RandomVariate switches between acceptance‐rejection methods due to Jöhnk [11], Cheng [12], and Atkinson [13]. An example of the advantage of using an acceptance‐rejection method over construction from two gammas can be seen in the following. The direct acceptance‐rejection method is nearly twice as fast as the gamma‐pair construction.

For StudentTDistribution the method used by RandomVariate is a polar rejection method due to Bailey [14]. This method is more efficient than direct construction from normal and variates, as can be seen in the following. The direct construction takes roughly 1.5 times as long as the polar method for a million Student variates.

This shows a comparison of direct construction and Bailey's polar rejection method for Student :

Discrete Distributions

GeometricDistribution, BetaBinomialDistribution, and BetaNegativeBinomialDistribution use direct construction. GeometricDistribution variates are generated as where follows UniformDistribution[0,1]. BetaBinomialDistribution and BetaNegativeBinomialDistribution are constructed from BinomialDistribution and NegativeBinomialDistribution variates with probability parameters taken as random BetaDistribution variates.

When used, table lookups for random integer generation are implemented via RandomChoice using the distribution's probability mass function for the weights. Most discrete distributions switch to other methods whenever construction of the list of weights is expected to be expensive given the desired sampled size. For example, as approaches 1 LogSeriesDistribution[p] switches to the direct construction , where and are uniformly distributed on the interval [15]. Depending on parameters and sample size, NegativeBinomialDistribution[n,p] may switch to construction as a Poisson-gamma mixture, which is a Poisson variable with mean following a gamma distribution [6].

The acceptance‐rejection method used by PoissonDistribution is due to Ahrens and Dieter [18]. The acceptance and rejection is carried out using discrete normal variates, taking advantage of the tendency of PoissonDistribution[μ] toward NormalDistribution[μ,]as μ increases.

Defining Distributional Generators

Rules for generating random values from distributions are generally defined via a TagSet on the head of the distribution. The distribution itself may contain parameters. As a simple example, the following defines rules for NegativeOfUniform[a,b], which represents a uniform distribution on the interval .

While the previous examples show the basic framework for defining distributional generators, the distributions themselves are not particularly interesting. In fact, it would have been easier in these two cases to just generate values from RandomVariate and multiply the end result by instead of attaching definitions to a new symbol. The following examples will demonstrate slightly more complicated distributions, in which case attaching definitions to a symbol will be more useful.

Example: Normal Distribution by Inversion

• generate a uniform random number on the interval

• compute the inverse cumulative distribution function of for the distribution of interest
To demonstrate the process, use the normal distribution. To generate random normal variates using this method, start with the Quantile for the normal distribution at a point and replace with a random uniform between 0 and 1.
Here is the Quantile function for a normal with mean and standard deviation :
A new distribution object can now be used to define a normal random number generator that uses inversion of the cumulative distribution function.
This defines generation of normals by inversion:
Here are 10 random normals generated by inversion:
Here is a sample of 10⁴ random normals along with the sample mean and standard deviation:
The normal distribution is one example where methods other than direct inversion are generally preferred. While inversion of the CDF is a perfectly valid method for generating pseudorandom normal numbers, it is not particularly efficient. Numeric evaluation of InverseErf is computationally much more expensive than the sinusoid and logarithmic evaluations required by the Box–Müller method used by NormalDistribution.
The built-in method takes almost no time for the same number of values:

Example: Uniform Distribution on a Disk

Random`DistributionVector can also be used to define generators for multidimensional distributions. For instance, suppose a random point from a uniform distribution on the unit disk, the set of real points with , is desired. Such a random point can be constructed as follows:

• generate a random angle uniformly distributed on

• generate a random vector uniformly distributed on

• return
The returned ordered pair can be multiplied by to generate points uniformly distributed on a disk of radius .
The following defines a generator for a uniform disk of radius :
Here is one random ordered pair from a disk of radius 2:
The following visualizes the distribution of 10⁴ generated points on this disk:

Example: Gibbs Sampler

Here are 5 vectors from the sampler with and :

[1] Geman, S. and D. Geman. "Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images." IEEE Transactions on Pattern Analysis and Machine Intelligence 6, no. 6 (1984): 721–741.

[2] Casella, G. and E. I. George. "Explaining the Gibbs Sampler." The American Statistician 46, no. 3 (1992): 167–174.

[3] Matsumoto, M. and T. Nishimura. "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudorandom Number Generator." ACM Transactions on Modeling and Computer Simulation 8, no. 1 (1998): 3–30.

[4] Nishimura, T. "Tables of 64-Bit Mersenne Twisters." ACM Transactions on Modeling and Computer Simulation 10, no. 4 (2000): 348–357.

[5] Junod, P. "Cryptographic Secure Pseudo-Random Bits Generation: The Blum–Blum–Shub Generator." August 1999. http://crypto.junod.info/bbs.pdf

[6] Gentle, J. E. Random Number Generation and Monte Carlo Methods, 2nd ed. Springer-Verlag, 2003.

[7] Johnson, N. L., S. Kotz, and N. Balakrishnan. Continuous Univariate Distributions, Volume 2, 2nd ed. John Wiley & Sons, 1995.

[8] Smith, W. B. and R. R. Hocking. "Algorithm AS 53: Wishart Variate Generator." Applied Statistics 21, no. 3 (1972): 341–345.

[9] Cheng, R. C. H. and G. M. Feast. "Some Simple Gamma Variate Generators." Applied Statistics 28, no. 3 (1979): 290–295.

[10] Johnson, M. E. Multivariate Statistical Simulation. John Wiley & Sons, 1987.

[11] Jöhnk, M. D. "Erzeugung von Betaverteilten und Gammaverteilten Zufallszahlen." Metrika 8 (1964): 5–15.

[12] Cheng, R. C. H. "Generating Beta Variables with Nonintegral Shape Parameters." Communications of the ACM 21, no. 4 (1978): 317–322.

[13] Atkinson, A. C. "A Family of Switching Algorithms for the Computer Generation of Beta Random Variables." Biometrika 66, no. 1 (1979): 141–145.

[14] Bailey, R. W. "Polar Generation of Random Variates with the t-Distribution." Mathematics of Computation 62, no. 206 (1994): 779–781.

[15] Devroye, L. Non-Uniform Random Variate Generation. Springer-Verlag, 1986.

[16] Kachitvichyanukul, V. and B. W. Schmeiser. "Binomial Random Variate Generation." Communications of the ACM 31, no. 2 (1988): 216–223.

[17] Kachitvichyanukul, V. and B. W. Schmeiser. "Computer Generation of Hypergeometric Random Variates." Journal of Statistical Computation and Simulation 22, no. 2 (1985): 127–145.

[18] Ahrens, J. H. and U. Dieter. "Computer Generation of Poisson Deviates from Modified Normal Distributions." ACM Transactions on Mathematical Software 8, no. 2 (1982): 163–179.

[19] Matsumoto, M. and T. Nishimura. "Dynamic Creation of Pseudorandom Number Generators." In Proceedings of the Third International Conference on Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing: Monte Carlo and Quasi‐Monte Carlo Methods 1998, 56–69, 2000.