Lesson #7: Stochastic Simulation

Generating Pseudo-Uniforms (Our own)

a=1664525;(*1664525*)
c=1013904223;(*1013904223*)
m=(2^32);
next[x_]:=Mod[a x + c,m]

x0=3;
nn=(10^5);
lst=N[(1/m)NestList[next,x0,nn]];
Mean[lst]
Variance[lst]
Variance[UniformDistribution[{0,1}]]//N
((Covariance[Drop[lst,1],Drop[lst,-1]])/(Variance[lst]))
Manipulate[ListPlot[Take[lst,{n,n+1000}],Ticks->{False,Automatic}],{n,1,nn-1000}]
ListPlot[(1/nn)BinCounts[lst,{0,1,0.01}],Filling->Axis,PlotRange->{0,0.01*1.1}]

Using Mathematica's RandomReal[] and SeedRandom[]

SeedRandom[3];
lst=Table[RandomReal[],{nn}];
Mean[lst]
Variance[lst]
((Covariance[Drop[lst,1],Drop[lst,-1]])/(Variance[lst]))

The Inverse Probability Transform

λ=1/2;

inv[x_]:=InverseCDF[ExponentialDistribution[λ],x]

inv[x]

nn=500;

unifs=Table[RandomReal[{0,1}],{nn}];

<<Histograms`

Histogram[Map[inv,unifs]]

element[un_]:={PointSize[Medium],Point[{0,un}],Line[{{0,un},{inv[un],un}}],Point[{inv[un],un}],Line[{{inv[un],un},{inv[un],0}}]}

grps=Map[element,unifs];

Manipulate[
GraphicsColumn[{Plot[CDF[ExponentialDistribution[1/2],x],{x,0,10},
Epilog->grps[[n]]],
ListPlot[(1/n)BinCounts[Map[inv,Take[unifs,n]],{0,10,0.5}],AxesOrigin->{0,0},Joined->False,PlotRange->{0,0.35}]}],
{n,1,nn,1}]

Specialized Methods

Map[InverseCDF[NormalDistribution[],#]&,{0,0.25,0.5,0.75,1}]

Timing[Table[InverseCDF[NormalDistribution[],RandomReal[{0,1}]],{10000}];]

myNormal[]:=Module[{},
angle=2 Pi RandomReal[];
radius = Sqrt[-Log[RandomReal[]]];
{radius Cos[angle],radius Sin[angle]}
]

Timing[Table[myNormal[][[1]],{10000}];]

Timing[Table[RandomReal[NormalDistribution[]],{10000}];]

ListPlot[Table[myNormal[],{10000}] ,AspectRatio->Automatic,PlotRange->{{-3,3},{-3,3}}]

Scaling Random Variables

Manipulate[ 
SeedRandom[seed];
ListPlot[
Map[μ +σ#&,Table[myNormal[],{100}]]
,AspectRatio->Automatic,PlotRange->{{-3,3},{-3,3}}],{{μ,{0,0}},{-2,-2},{2,2}},
{{σ,{1,1}},{0,0},{2,2}},
{seed,1,5,1},ControlPlacement->Left]

rvs=Table[myNormal[],{1000}];

Manipulate[ 
ListPlot[Map[μ +σ#&,rvs],AspectRatio->Automatic,PlotRange->{{-3,3},{-3,3}},
,{{μ,{0,0}},{-2,-2},{2,2}},
{{σ,{1,1}},{0,0},{2,2}},Button["New Sample",(rvs=Table[myNormal[],{1000}])&],
ControlPlacement->Left
]

Monte Carlo Integration

f[x_]:=30+((((((x^2)+(Sqrt[((Cos[x])/(Sin[x]))]))/(Log[x]))))^7)

Plot[f[x],{x,0.05,0.3},Filling->Axis,AxesOrigin->{0,0}]

Integrate[f[x],{x,0.05,0.3}]

NIntegrate[f[x],{x,0.05,0.3}]

(0.30-0.05)*15

Table[
(Table[f[RandomReal[{0.05,0.30}]],{10000}]//Mean)*(0.30-0.05)
,
{10}]//MatrixForm

f[x_,y_,z_,w_]:=Cos[100x]Cos[100y]Cos[100z]Cos[100w]

Integrate[f[x,y,z,w],{x,0,1},{y,0,1},{z,0,1},{w,0,1}]

%//N

NIntegrate[f[x,y,z,w],{x,0,1},{y,0,1},{z,0,1},{w,0,1}]

Table[Apply[f,RandomReal[{0,1},4]],{10000}] //Mean

Class Exercise: Random Process Simulation

Random Walks

Markov Chains

Markov Jump Processes