Random Number Generation in Mathematica 6

Rob Knapp

Wolfram Research, Inc.

Abstract

The framework for pseudorandom number generation in Mathematica has been
completely redesigned for version 6.  There are new functions optimized
for generating large arrays of random numbers as well as for making
random choices and samples from sets.  A mechanism is provided to extend
these functions for distributions beyond the several built-in distributions.
There are several new random generators along with a way to switch between
generators.

Goals

Functions for random numbers

Three new functions were introduced for getting random numbers in Mathematica, RandomInteger, RandomReal, and RandomComplex.  In addition to the range argument that Random has, they all accept a second argument that gives the dimensions of the returned array or random numbers.

With a range, the functions above all generate numbers that are uniformly distributed.  For RandomInteger and RandomReal you can give a distribution.

RandomInteger

The default is to generate a single number equally probable to be 0 or 1.

RandomInteger[]

A 4x4 matrix of 0's and 1's

MatrixForm[RandomInteger[1, {4,4}]]

Estimate the probability that a 4x4 matrix of ones and zeros will be singular.

Timing[
Block[{sing = 0},
Do[If[Det[RandomInteger[1, {4,4}]] == 0, sing++],{10^5}];
N[sing/10^5]
]
]

The same but a little bit faster.

Timing[
cf = Compile[{{n, _Integer}}, Block[{sing = 0},
Do[If[Det[RandomInteger[1, {4,4}]] == 0, sing++],{n}];
N[sing/n]
]];
cf[10^5]
]

Generate 5 random decimal digits.

digits = RandomInteger[{0,9}, 5]

Convert it to a string.

sd = Apply[StringJoin, Map[ToString, digits]]

The positions the digits appear in the first million digits of Pi.

pistring = ToString[FractionalPart[N[Pi, 10^6]]];
StringPosition[pistring,sd]

Generate a million random bits.

RandomInteger[{0, 2^(10^6 - 1)}]

A list of numbers distributed according to the Poisson distribution. The probability for integer value x in a Poisson distribution is (e^(-μ))(μ^x)/x! for x≥0.

RandomInteger[PoissonDistribution[10], 10]

Show convergence towards the PDF as more numbers are generated.

Table[n = 10^k;
p = RandomInteger[PoissonDistribution[10], n];
max = Max[p];
Show[{ListPlot[Transpose[{Range[0,max],BinCounts[p,{0,max + 1}]/n}]],
Plot[Exp[-10] 10^x/Gamma[x + 1],{x,0,max}]}, PlotLabel->ToString[10^k]],
{k,1,6}]

Generating arrays is much faster than individual numbers.

Timing[Table[RandomInteger[], {10^7}];]

Timing[RandomInteger[1, {10^7}];]

RandomReal

The default is to give a real uniformly distributed between 0 and 1.

RandomReal[]

Generate a 20 points in the unit rectangle 2 dimensions.

pts = RandomReal[1,{20,2}]

Generate 10 values between -1 and 1.

vals = RandomReal[{-1,1}, 20]

Plot the polynomial that interpolates those points.

data = Transpose[{pts, vals}];
poly[x_, y_] = InterpolatingPolynomial[data, {x,y}];
Show[{Plot3D[poly[x,y],{x,0,1},{y,0,1}, PlotRange->All],Graphics3D[{PointSize[0.025], Yellow, Map[Point[Flatten[#]]&, data]}]}]

Generate a 5x5 matrix of values with the standard normal distribution.

MatrixForm[rm = RandomReal[NormalDistribution[0,1], {5,5}]]

Estimate the probability that such a 5x5 matrix will have real eigenvalues.

Block[{count = 0}, 
Do[ev = Eigenvalues[RandomReal[NormalDistribution[0,1], {5,5}]]; If[Re[ev] == ev,count++],{10^5}];
N[count/10^5]]

Generate 10 numbers with precision 100 between 1 and 2.

RandomReal[{1,2}, 10, WorkingPrecision->100]

Estimate the worst precision loss from computing (Subscript[J, 1]) on the interval [0, 1000].

100 - Precision[BesselJ[1, RandomReal[{0,1000}, 1000, WorkingPrecision->100]]]

Generate 9 numbers with precision 27 distributed with the t-distribution.

RandomReal[StudentTDistribution[2], 9, WorkingPrecision->27]

Generating arrays is much faster than generating individual values.

Timing[Table[RandomReal[],{10^7}];]

Timing[RandomReal[1, {10^7}];]

The default is to give a number in the unit rectangle in the complex plane.

RandomComplex[]

Generate a list of 10 numbers on the imaginary axis:

RandomComplex[{-I,I}, 10]

Generate a list of numbers in a square centered at the origin and show them on a plot.

list = RandomComplex[{-1-I, 1 + I}, 10^4];
ListPlot[Transpose[{Re[list], Im[list]}], AspectRatio->1]

Test an identity for complex values.

Max[Abs[Map[Function[{x}, Sqrt[x^2] - Abs[x]],list]]]

Test another identity for complex values.

Max[Abs[Map[Function[{x},Sin[x]^2 + Cos[x]^2 - 1],list]]]

Generating arrays is much faster than generating individual values.

Timing[Table[RandomComplex[],{10^7}];]

Timing[RandomComplex[1, {10^7}];]

Functions for random selection

RandomChoice and RandomSample are functions that make random selections from lists with and without replacement, respectively.

Both functions allow giving weights that specify the relative probability an element will be selected.

RandomChoice

Four letter "words" formed by selecting at random from the lower case alphabet.

alphabet = Table[FromCharacterCode[i], {i,97,122}]

RandomChoice[alphabet, {10,4}]

Map[StringJoin, %]

Four letter lower case words in English.

lowercaseQ[w_] := Block[{cc = ToCharacterCode[w]},(Min[cc] >= 97) && (Max[cc] <= 122)];Length[fourletterwords =Select[DictionaryLookup["*"], ((StringLength[#] == 4) && lowercaseQ[#])&]]

Table[FourLetterWordQ[w] = True, {w, fourletterwords}];
FourLetterWordQ[_] = False;

Estimate the probability of getting an actual English four letter word.

Block[{count = 0},Table[count += Boole[FourLetterWordQ[w]], {w, Map[StringJoin, RandomChoice[alphabet, {10^5, 4}]]}]; N[count/10^5]]

Weights for the letters according to how often they appear in the four letter words.

chars = Apply[StringJoin, fourletterwords];
weights= Map[StringCount[chars, #]&, alphabet]

Estimate the probability of getting an actual English four letter word with letters chosen according to the weights.

Block[{count = 0},Table[count += Boole[FourLetterWordQ[w]], {w, Map[StringJoin, RandomChoice[weights->alphabet, {10^5, 4}]]}]; N[count/10^5]]

RandomSample

With just a list, the default is to give a random permutation.

RandomSample[Range[5]]

Estimate probability of getting n elements in order

Block[{trials = 10^5, count, rn},
Table[
count = 0;
rn = Range[n];
Do[If[RandomSample[rn] == rn, count++], {trials}]; 
{n, N[count/trials]},
{n, 2, 8}]]

Estimate probability of getting n elements in order when weighted

Block[{trials = 10^5, count, rn, w},
Table[
count = 0;
rn = Range[n];
w = Reverse[rn];
Do[If[RandomSample[w->rn] == rn, count++], {trials}]; 
{n, N[count/trials]},
{n, 2, 8}]]

A big list of points in the unit rectangle in 2 dimensions.

r = RandomReal[1, {10^6, 2}];

Compute the mean distance from the origin of a set of points.

dmean[pts_] :=Mean[Sqrt[Total[pts^2, {2}]]];

μ = dmean[r]

The standard deviation of the distance from the origin.

σ = StandardDeviation[Sqrt[Total[r^2, {2}]]]

A demonstration of the central limit theorem.

Block[{m = 10000, x}, 
Table[
n = 10^k;
pdf =PDF[ NormalDistribution[μ, σ/Sqrt[n]]][x];
samples = Table[RandomSample[r, n], {m}];
smeans = Map[dmean, samples];
{max, min} = {Max[smeans], Min[smeans]};
dx = (max - min)/100;
bins = Range[min, max,dx];
Show[{ListPlot[Map[{Mean[#], Length[#]/(m dx)}&, BinLists[smeans, {min, max, dx}]], PlotRange->All],
Plot[Evaluate[pdf],{x,min, max}]}, PlotLabel->n],
{k,0,3}]
]

Controlling the random number generator

The function SeedRandom is similar to before, but it now also has a method option.

The function BlockRandom allows you to localize the effects of using random number generators.

BlockRandom

Get a list of 1 bit numbers.

BlockRandom[RandomInteger[{0,1}, 8]]

Get a list of two bit numbers.

BlockRandom[RandomInteger[{0,3}, 4]]

Map[PadLeft[IntegerDigits[#,2],2,0]&, %]

Get a list of 1 bit numbers.

RandomInteger[{0,1}, 8]

Get a list of two bit numbers.

RandomInteger[{0,3}, 4]

Map[PadLeft[IntegerDigits[#,2],2,0]&, %]

SeedRandom

The default is to seed with a number derived from the time.

SeedRandom[];
s128bit = RandomInteger[2^128 - 1]

If you want to be able to reproduce calculations, seed with a particular value.

SeedRandom[192673033350623974216981134157370125810];
{RandomReal[{-1,1},5], RandomInteger[1,5]}

SeedRandom[192673033350623974216981134157370125810];
{RandomReal[{-1,1},5], RandomInteger[1,5]}

Methods

There are now several different methods that you can use for generating the basic sequences for the functions.

Here is a comparison of the speed of some of the generators for generating 10 million random reals.

TableForm[Table[
BlockRandom[
SeedRandom[Method->method];
{method, First[Timing[RandomInteger[2^16 - 1, 10^7];]],First[Timing[RandomReal[1, 10^7];]]}],
{method, {"ExtendedCA", "Rule50025CA", "Rule30CA", "MersenneTwister", "Congruential", "MKL"}}],
TableHeadings->{{},{"Method", "Integers", "Reals"}}]

Many of the methods have parameters that control or parameterize them

Generate the linear congruential sequence (Subscript[x, n+1])= 1+13 (Subscript[x, n]) mod 47 converting the elements directly into reals.

BlockRandom[
SeedRandom[1, Method->{"Congruential", "Increment"->1, "Modulus"->47, "Multiplier"->13, "ConvertToRealsDirectly"->True}];
RandomReal[47, 10]
]

Generate a Niederreiter low discrepancy sequence for dimension 2.

ListPlot[BlockRandom[
SeedRandom[Method->{"MKL", Method->{"Niederreiter", "Dimension"->2}}];RandomReal[1,{10000,2}]], 
AspectRatio->1]

Compare with the default generator

ListPlot[BlockRandom[RandomReal[1,{10000,2}]], AspectRatio->1]

Extensibility

The new framework allows extensibility for random distributions and basic random number generators.  Adding a basic generator is somewhat complicated, so I will not cover it in this talk.  Adding a new distribution is quite easy and is done by adding a rule for  Random`DistributionVector.

Defining Distributions

It is a good idea to associate the rules for a new distribution with the new distribution name.  This can be done by using UpSet

Define a distribution that generates values based on selecting k coins from a pile of cash with p pennies, etc, and adding up their value (in cents).

coins /: Random`DistributionVector[coins[k_, {p_, n_, d_, q_}], m_, ___] := Module[{cash = Flatten[MapThread[Table[#1, {#2}]&,{ {1,5,10,25}, {p,n,d,q}}]]},
If[k > Length[cash], 
Message[coins::greed];
Throw[$Failed],
Table[Total[RandomSample[cash, k]], {m}]]]

coins::greed = "Don't be so greedy";

A single number from this distribution.

RandomInteger[coins[10, {100, 20, 10, 4}]]

A 10x10 matrix

MatrixForm[RandomInteger[coins[10, {100, 20, 10, 4}], {10, 10}]]

Too greedy.

RandomInteger[coins[1000,  {100, 20, 10, 4}]]

Acceptance-rejection method for uniformly distributed values in the circle of radius r in the complex plane.

ucircle /: Random`DistributionVector[ucircle[r_], n_, prec_] := 
Module[{left = n, c, s},
Reap[While[left > 0,
s= RandomComplex[r {-1 -I, 1 + I}, left, WorkingPrecision->prec];
Table[If[Abs[v] <= r, left--;Sow[v]],{v, s}];
]][[2,1]]]

Generate a ten thousand sample points and plot them.

c = RandomComplex[ucircle[1], 10^4];
ListPlot[Transpose[{Re[c], Im[c]}], AspectRatio->1]

The Dirichlet distribution

Dirichlet Distribution

DirichletDistribution /: Random`DistributionVector[DirichletDistribution[α_?(VectorQ[#, Positive]&)], n_, prec_] := 
Module[{k, theta, normalization},
k = Length[α];
theta = Map[RandomReal[GammaDistribution[#,1], n, WorkingPrecision->prec]&, α];
normalization = Total[theta];
Transpose[theta]/normalization
]

RandomReal[DirichletDistribution[{1,2}], 10]