Supplemental Notebook for Monte Carlos

Random Number Generator

The following program generates Nev random numbers, uniformly distributed on [0,1].  We just use the default random number generator in Mathematica (which is the Mersenne Twister).

rand[Nev_]:=Module[{randtab,r},
randtab={};
For[i=1,i≤Nev,i++,
r=RandomReal[{0,1}];
AppendTo[randtab,r];
];
Return[randtab];
];

As an example of its use, here are 10 random numbers:

rand[10]

Sampling Exponential Distribution

The following program samples the exponential distribution with half-life of λ with Nev events.  The program outputs a histogram of Nbins bins that ranges from 0 to maxpoint.  The distribution is sampled by generating a random number on [0,1] and inverting the cumulative distribution.

expdist[λ_,Nev_,maxpoint_,Nbins_]:=Module[{randtab,r,t,binno},
randtab=Table[{maxpoint*i/Nbins,0.0},{i,0,Nbins-1}];
For[i=1,i≤Nev,i++,
r=RandomReal[{0,1}];
t=-λ*Log[1-r]/Log[2];
If[t<maxpoint,binno=Ceiling[Nbins*t/maxpoint];
randtab[[binno,2]]+=Nbins/(maxpoint*Nev);
];
];
Return[randtab];
];

Let’s see how this works.  We’ll set the half-life to be 1s, we’ll generate 10000 events, we’ll plot the histogram out to 7s and it will have 40 bins.  We can also compare our histogram to the exact expression for the probability distribution.

tab=expdist[1,10000,7,40];
histplot=ListLinePlot[tab,InterpolationOrder->0,PlotStyle->Black,Frame->True,FrameStyle->Thick,FrameLabel->{"t","p(t)"}];
trueplot=Plot[Log[2](E^(-(t/1)Log[2.])),{t,0,7},PlotStyle->{Black,Dashed},Frame->True,FrameStyle->Thick,FrameLabel->{"t","p(t)"}];
Show[{histplot,trueplot}]

Pretty good agreement!  The agreement will improve when more events are included in the pseudodata sample.

Sampling a 2D distribution

Here is a program that samples the 2D distribution p(x,y) = x + y, as an illustration of using conditional probability.  The function xygen generates Nev events that are distributed according to the probability distribution.  To do this requires calculating and inverting the cumulative distribution of y and the conditional distribution of x given y.

xygen[Nev_]:=Module[{outtab,r1,r2,x,y},
outtab={};
For[i=1,i≤Nev,i++,
r1=RandomReal[{0,1}];
r2=RandomReal[{0,1}];
y=-1/2+(Sqrt[1/4+2*r1]);
x=-y+(Sqrt[y*y+(2*y+1)*r2]);
AppendTo[outtab,{x,y}];
];
Return[outtab];
];

Sampling a Distribution with the Veto Method

The following function vetowell samples the probability distribution of the position of the particle in the ground state of the infinite square well using the veto method.  Its inputs are the number Nev events to generate and the number of bins Nbins to put in the output histogram.  The function bounds the probability distribution by a constant, and first finds a random number on [0,1] to determine the position being considered.  Then, it picks another random number to determine if the event is kept, according to the value of the probability distribution at that point.

vetowell[Nev_,Nbins_]:=Module[{outtab,r1,r2,binno},
outtab=Table[{i/Nbins,0.0},{i,0,Nbins-1}];For[i=1,i≤Nev,i++,
r1=RandomReal[{0,1}];
r2=RandomReal[{0,1}];
binno=Ceiling[r1*Nbins];
If[r2<((Sin[π*r1])^2),outtab[[binno,2]]+=2.0*Nbins/Nev];
];
Return[outtab];
];

We can test this out and see how well it agrees with the exact probability distribution.  Here, we generate 1000000 events and output a histogram with 100 bins.  The histogram from the Monte Carlo is plotted in solid black while the exact result is dashed.

plottab=vetowell[1000000,100];
plota=ListLinePlot[plottab,InterpolationOrder->0,PlotStyle->Black,Frame->True,FrameStyle->Thick,FrameLabel->{"x","(SuperStar[ψ])ψ"}];
plotb=Plot[2((Sin[π x])^2),{x,0,1},PlotStyle->{Black,Dashed},Frame->True,FrameStyle->Thick,FrameLabel->{"t","p(t)"}];
Show[plota,plotb]

Excellent agreement between our Monte Carlo and the exact result!