MCMC

SetDirectory[NotebookDirectory[]]

Resource Functions

The Wolfram Function Repository is a curated cloud repository of functions set up to be instantly usable in the Wolfram Language . The Function Repository includes a growing number of functions designed for a wide range of application areas . The Wolfram Language supports creation of private resource functions, which can then be submitted for inclusion in the Wolfram Function Repository .

https://resources.wolframcloud.com/FunctionRepository

Metropolis Hastings

(* install the RF permamently *)
ResourceFunction["PersistResourceFunction"]["MetropolisHastingsSequence"]

ResourceFunction["MetropolisHastingsSequence"]

(* generate a MCMC sampling of a likelihood; set initial value, plus burn-in "0" (no burn-in) and thinning "1" (no thinning)*)

mc=[Exp[-1./2.(x^2.+y^2.)],{{x,0},{y,0}},{50000,0,1}];

mc[[1;;10]]

Animate[ListPlot[mc[[1;;i]],Frame->True,PlotRange->{{-4,4},{-4,4}}],{i,1,300,1},AnimationRunning->False]

ListPlot[mc,Frame->True]

Histogram3D[mc]

SmoothHistogram3D[mc]

(* convert into a prob distribution *)

dist=SmoothKernelDistribution[mc]

Plot3D[PDF[dist,{x,y}],{x,-3,3},{y,-3,3},PlotRange->All,PlotPoints->50]

ContourPlot[PDF[dist,{x,y}],{x,-3,3},{y,-3,3},PlotRange->All]

Covariance[dist]//MatrixForm

(* generate a MCMC sampling of a likelihood *)

cov={{1,.5},{.5,1}};invcov=Inverse[cov];var={x,y};mc=[Exp[-1./2.*var.invcov.var],{{x,0},{y,0}},{50000,0,1}];
(* convert into a prob distribution *)
dist=SmoothKernelDistribution[mc];
ContourPlot[PDF[dist,{x,y}],{x,-3,3},{y,-3,3},PlotRange->All]
Covariance[dist]//MatrixForm

max=NMaximize[PDF[dist,{x,y}],{x,y}]
maxl=max[[1]]

(* integrate over a region high than a given level *)

level=0.15;r=ImplicitRegion[PDF[dist,{x,y}]>level && -2<x<2 && -2<y<2,{x,y}];
NIntegrate[PDF[dist,{x,y}],{x,y}∈r]

(* do that for several levels *)

res={};tablev=Table[level,{level,0,.95*maxl,maxl/5}];
Do[{r=ImplicitRegion[PDF[dist,{x,y}]>tablev[[i]] && -2<x<2 && -2<y<2,{x,y}];
AppendTo[res,{NIntegrate[PDF[dist,{x,y}],{x,y}∈r],tablev[[i]]}]},{i,Length[tablev]}]

res

hei=Interpolation[res]

ContourPlot[{PDF[dist,{x,y}]==hei[.68],PDF[dist,{x,y}]==hei[.9]},{x,-3,3},{y,-3,3},PlotRange->All]