Random numbers

A random variable is a variable that takes a value with a certain probability. A random variable is always associated with a probability distribution.

An example is a random variable U[0, 1) that takes a real value between 0 and 1 EQUALLY. Meaning of “equally” is mathematically described as follows using p(x) as the probability density function of X = U[0, 1).

Prob[ x < X < x + dx] = p(x) dx

with sufficiently small interval dx.

A (pseudo) random number can be generated by a random number generator. “Mersenne twister” is a good algorithm to generate U[0, 1) (C programing).

But in this lecture we use a Mathematica built-in command “Random[]” to generate U[0, 1) and exercise that various random numbers can be generated using U[0, 1).

Random numbers uniformly distributed U[a, b) (Mathematica)

Random[]

Table[Random[],{10}]

data=Table[Random[],{100000}];

data

Min[data]
Max[data]

Mean[data]

Variance[data]

Histogram[data,{0,1,0.05}]

Histogram[data,{0,1,0.05},"PDF"]

Random[Real,{1,3}]

Table[Random[Real,{1,3}],{10}]

data=Table[Random[Real,{1,3}],{100000}];

Histogram[data,10,"PDF",PlotRange->{{0,5},All},AxesOrigin->{0,0}]

You will notice the property :

U[0, 1) a + b = U[b, b + a)

where equal “=” means equality in distribution. Once you can generate U[0, 1), you can generate random numbers that are uniformly distributed in arbitrary range.

Random numbers exponentially distributed Exp[r] (Mathematica)

Y = –(log U (0, 1])/g follows the exponential distribution with parameter g > 0. Let’s see this property heuristically.

data=Table[Random[],{100000}];

Take[data,10]

data2=-Log[data]/1.0;

Take[data2,10]

Min[data2]
Max[data2]

Mean[data2]

Variance[data2]

ghist=Histogram[data2,{0,13,0.2},"PDF"]

g=Plot[gamma Exp[- gamma x]/.gamma->1.0,{x,0,20},PlotStyle->RGBColor[1,0,0],PlotRange->All]

Show[ghist,g]

Random numbers with Gaussian (Normal) distribution N(m, sd^2) (Mathematica)

In Mathematica, Random[NormalDistribution[0, 1]] generates a random number that follows N (0, 1).

data=Table[Random[NormalDistribution[0,1]],{10000}];

Min[data]
Max[data]

Mean[data]

Variance[data]

ghist=Histogram[data,50,"PDF"]

gauss=Exp[-(x-m)^2/2/sd^2]/Sqrt[2Pi]/sd/.{m->0, sd->1}
g=Plot[gauss,{x,-4,4},PlotStyle->RGBColor[1,0,0],PlotRange->All]

Show[ghist,g]

Random number of Normal distribution N (0, 1) can be generated by Box Muller method.

seq=Table[ Sqrt[-2Log[Random[]]]Cos[2Pi Random[]],10000]

Min[seq]
Max[seq]

Mean[seq]

Variance[seq]

ghistBM=Histogram[seq,50,"PDF"]

Show[ghistBM,g]

This produces N (3,4)

seq2=seq*2 + 3

Mean[seq2]

Variance[seq2]

ghistBM2=Histogram[seq2,50,"PDF"]

gauss=Exp[-(x-m)^2/2/sd^2]/Sqrt[2Pi]/sd/.{m->3, sd->2}
g2=Plot[gauss,{x,-4,12},PlotStyle->RGBColor[1,0,0],PlotRange->All]

Show[ghistBM2,g2]

Bernoulli distribution

Using the uniformly distributed random number X = U[0, 1), we can implement a Bernoulli distribution X = 0 or 1 with the rule: X = 0 with probability P and X = 1 with 1-P. Flipping a coin corresponds to a Bernoulli process (head or tail) with Bernoulli distribution with P.

myBernoulli[probP_]:=Module[{rvUni, result},
rvUni=Random[];
If[rvUni < probP,result=0,result=1];
Return [result];
]

myBernoulli[0.4]

data=Table[myBernoulli[0.4],{1000}]

Count[data,0]
Count[data,1]

Mean[data]

Variance[data]

Histogram[data,Automatic,"PDF"]

Binomial distribution

Sum of Bernoulli realization results in Binomial distribution. Let’s check this. As an example, assume a distorted coin; we have head 0 with probability P = 0.4 and tail 1 with P =0.6. We toss the coin 10 times and count the number of having head 0.

data=Table[myBernoulli[0.4],{10}]

Count[data,0]

If we repeat this process for a large number of realization, say 1000 times, we have the following result:

dataBinomial={};
Do[
data=Table[myBernoulli[0.4],{10}];
AppendTo[dataBinomial,Count[data,0]],{1000}
]

dataBinomial

Mean[dataBinomial]

Variance[dataBinomial]

gSim=Histogram[dataBinomial,Automatic,"PDF"]

gPDF=DiscretePlot[PDF[BinomialDistribution[10,0.4],n],{n,0,10},PlotRange->All,PlotMarkers->Automatic]

Show[gSim,gPDF]

Poisson process and Poisson distribution

Events that occur with the rate R (R > 0) occurs with time interval (inter-event time or waiting time) T that follows the exponential distribution with parameter R.

rateR=2

seq=Table[-Log[Random[]]/rateR,{10}]

The number of events that occurred within the UNIT time follows the Poisson distribution with parameter R. The rate R corresponds to the average number of events that occurs within the unit time. 

The number of events that actually occurred, N, can be calculated as follows using the sequence of interevent times. Naturally, N is non-negative integer.

seq

index = 1; (* List index starts from 1 in Mathematica *)
sumT = 0.0;
sumT+=seq[[index]];
count=0;
While[sumT<=1,
index++;
sumT+=seq[[index]];
count++;
]
count

This process can be defined as a function of Mathematica.

myPoisson[rateR_]:=Module[{seqIntereventTime={},index=1,sumT=0.0,count=0},
seqIntereventTime=Table[-Log[Random[]]/rateR,{100}];
sumT+=seqIntereventTime[[index]];
count=0;
While[sumT<=1,
index++;
sumT+=seqIntereventTime[[index]];
count++;
];
Return[count];
]

myPoisson[3]

data=Table[myPoisson[3],{1000}]

Mean[data]

Variance[data]

gSim=Histogram[data,Automatic,"PDF"]

gPDF=DiscretePlot[PDF[PoissonDistribution[3],n],{n,0,15},PlotRange->All,PlotMarkers->Automatic]

Show[gSim,gPDF]