Define a Probability Distribution
Define a Probability Distribution
Lets take a simple nontrivial non-Gaussian Distribution in two dimensions. E.g. a mixture of two Gaussian with different mean and variance:
Lets take a simple nontrivial non-Gaussian Distribution in two dimensions. E.g. a mixture of two Gaussian with different mean and variance:
In[]:=
mixturedist=MixtureDistribution,,{MultinormalDistribution[{-1,-1},{{1,0},{0,1}}],MultinormalDistribution[{2,5},{{2,0},{0,2}}]};
1
2
1
2
Plot the probability distribution, should give us two little hills in the plane.
Plot the probability distribution, should give us two little hills in the plane.
In[]:=
pdfxy[x_,y_]:=PDF[mixturedist,{x,y}]
In[]:=
myniceplot=Plot3D[pdfxy[x,y],{x,-10,10},{y,-10,10},PlotRange->All,PlotPoints->25,PlotTheme->"ZMesh"]
Out[]=
Define Markov Chain Monte Carlo Algorithm
Define Markov Chain Monte Carlo Algorithm
Algorithm
Algorithm
In[]:=
chainfun[whereweare_]:=Module{deltax,ru,ratio,whereweget}, Withrandomdelta:=RandomVariate[NormalDistribution[]], randomuniform:=RandomVariate[UniformDistribution[]], deltax={randomdelta,randomdelta}; ru = randomuniform; ratio = pdfxy@@(whereweare+deltax)pdfxy@@(whereweare); If[ru < ratio,whereweget = whereweare + deltax,whereweget=whereweare]
... that’s all
... that’s all
Let’s compile it for efficiency
Let’s compile it for efficiency
In[]:=
sampleMany=Compile[{{startpoint,_Real,1},{n,_Integer}},NestList[chainfun,startpoint,n]]
Out[]=
CompiledFunction
simple test
simple test
In[]:=
startpoint={0.,0.};
In[]:=
sampleMany[startpoint,4]
Out[]=
{{0.,0.},{0.,0.},{0.,0.},{-0.741135,-0.710525},{-1.39183,-0.950601}}
Sample from the distribution
Sample from the distribution
with a few points...
with a few points...
In[]:=
mypoints=sampleMany[startpoint,100];
In[]:=
mypointplot=ListPointPlot3D[Append[#,0.1]&/@mypoints]
Out[]=
In[]:=
Show[myniceplot,mypointplot]
Out[]=
this doesn’t look too promising, but that’s the “burn-in” period...
this doesn’t look too promising, but that’s the “burn-in” period...
sample a huge number of points, takes a minute...
sample a huge number of points, takes a minute...
In[]:=
manypoints=sampleMany[startpoint,100000];
chop off first 10.000 points, the burn-in period and plot
chop off first 10.000 points, the burn-in period and plot
In[]:=
manyendpoints=Take[manypoints,-90000];
In[]:=
manypointplot=ListPointPlot3D[Append[#,0.1]&/@manyendpoints]
Out[]=
In[]:=
Show[myniceplot,manypointplot]
Out[]=
reconstruct the distribution from the sampled points
reconstruct the distribution from the sampled points
Application
Application