Hastings - Metropolis Algorithm implementation For Problem 8.5 part(a)

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 problem 8.5 part(b)

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

SeedRandom[121212];

Define the data given in the problem

v={1,2,3,4};
edges={{1,2},{1,3},{2,3},{2,4},{3,4}};
f={1,8,6,4};

Define the functions π(x) and q(x,y) to use in the implementation

pi[x_,f_]:=((f[[x]])/((∑_(i=1)^(Length[f]))f[[i]]))

q[x_,y_]:=Module[{r},
r=Count[edges,{x,any_}]+Count[edges,{any_,x}];
If[(Count[edges,{x,y}]>0|| Count[edges,{y,x}]>0),1/r,0]
]

Find the stationary distribution

w=Table[pi[x,f],{x,1,v[[-1]]}]

Generate the original q matrix

(originalMatrix=Table[q[x,y],{x,1,v[[-1]]},{y,1,v[[-1]]}])//MatrixForm

Bring up the user interface

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

Define a function for cummulative sum

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];
(*Print["X=",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[maxN_]:=Module[{i,j,nStates,n,s,y,α,u,x,sample,pts,sum,k,found},
nStates=Length[originalMatrix];
n=1;
x=Table[0,{maxN}];
x[[n]]=1;  (*pick any state to start from*)
i=1;
While[i< maxN,
{
y=sampleFromQConditional[originalMatrix,x[[n]]];
α=beta[x[[n]],y,w,originalMatrix];
u=RandomReal[];
n++;
If[u≤ α,x[[n]]=y,x[[n]]=x[[n-1]]]; (*acceptance step*)
i++;
}
];

generatePMatrixFromStatePath[nStates,x]
]