Hastings - Metropolis Algorithm implementation

by Naser Abbasi. Mathematics 504, Spring 2008. CSUF

This below is an implementation of the Hastings - Metropolis algorithm. A simple GUI interface allows the user to specify the number of steps to run the algorithm for. At each step, the current P matrix and the current calculated stationary distribution for this P matrix are shown to help observe the convergence.

The input to this run below is that of example 8.3.1 from Lecture notes of Math 504 by Professor B. Gearhart, CSUF. Below I show the user interface, and few screen shots showing the progress of the convergence to the final P matrix.

Few seed the random number generator and display the q and the π distribution used

SeedRandom[121212];
MatrixForm[q={{0,.5,.5},{.5,0,.5},{.5,.5,0}}]

MatrixForm[pi={1/6,1/3,1/2}]

m=Manipulate[First@{x=hastings[pi,q,maxN]; Grid[{{"stationary distribution w=",N[pi]},{ "current stationary distribution=",MatrixPower[N[x],100][[1,All]]},{,},{"Current P Matrix=",N[MatrixForm[x]]}},Alignment->Left]},
{{maxN,1,"number of steps"},1,50000,25,ContinuousAction->False,Appearance->"Labeled"},
AutorunSequencing->{{1,300}}
]

This a snap shot of the above output for number of iterations

Appending

Algorithm implementation

cumSum[list_]:=Module[{i,sum,s,k},
sum=0;
k=Length[list];
s=Table[0,{k}];
For[i=1,i≤ k,i++,
{
sum=sum+list[[i]];
s[[i]]=sum;
}
];
s
]

Function to calculate β (x, y)

beta[x_,y_,pi_,q_]:=Module[{},
Min[1,((pi[[y]] q[[y,x]] )/(pi[[x]] q[[x,y]]))]
]

Function called at the end of the run to generate P from the path of states travelled

generatePMatrixFromStatePath[nStates_,x_]:=Module[{i,j,p,allPairs,n,m},
n=Length[x];
p=Table[0,{nStates},{nStates}];
allPairs = Partition[x, 2, 1];
For[i=1,i≤nStates,i++,
{
m=Count[allPairs, {i, y_}];
For[j=1,j≤nStates,j++,
If[m≠0,p[[i,j]]=Count[allPairs, {i, j}]/m,p[[i,j]]=0]
]
}
];

p
]

Function to sample from q using uniform distribution

sampleFromQConditional[q_,x_]:=Module[{s,found,j,k,sample,y},
s=Flatten[Position[q[[ x,All]],Except[0],1,Heads->False] ];
sample=q[[ x,s]];
sample=cumSum[sample];
y=RandomReal[];
found=False;
For[j=1,j≤Length[sample],j++,
If[ Not[found],If[y≤sample[[j]],{k=j;found=True}]]
];

y=s[[k]] 
]

The Hastings algorithm main loop

hastings[pi_,q_,maxN_]:=Module[{i,j,nStates,n,s,y,α,u,x,sample,pts,sum,k,found},
nStates=Length[q];
n=1;
x=Table[0,{maxN}];
x[[n]]=1;  (*pick any state to start from*)
i=1;
While[i< maxN,
{
y=sampleFromQConditional[q,x[[n]]];
α=beta[x[[n]],y,pi,q];
u=RandomReal[];
n++;
If[u≤ α,x[[n]]=y,x[[n]]=x[[n-1]]]; (*acceptance step*)
i++;
}
];

generatePMatrixFromStatePath[nStates,x]
]