Problem 8.5 part(b)

by Naser Abbasi. Mathematics 504, Spring 2008. CSUF

This below construct the P matrix using direct computation of the Hasting-Meropolis method

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)  to use in the implementation

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

Define the functions q(x,y)  to use in the implementation. This does something similar to the adjcancy matrix normally used. I used the Count[] function in Mathematica which automatically counts the edges from the edges list above, so there is really no need to construct an adjancy matrix as such.

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

define Beta function

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

define a function to calculate a non - diagonal entry in the P matrix

makeEntryInPMatrix[x_,y_,pi_,q_]:=Module[{},
If[q[[x,y]]==0,0,q[[x,y]] beta[x,y,pi,q]]
]

Construct the P matrix for the off-diagonal elements only

nStates=Length[originalMatrix];
p=Table[0,{nStates},{nStates}];

For[i=1,i≤ nStates,i++,
For[j=1,j≤nStates ,j++,
If[i≠ j,p[[i,j]]=makeEntryInPMatrix[i,j,w,originalMatrix]]
]
];

Now calculate the diagonal elements of the P matrix

For[i=1,i≤ nStates,i++,
For[j=1,j≤nStates ,j++,
If[i==  j,p[[i,j]]=1-Total[p[[i,All]]]]
]
];

Print the P matrix

p//N//MatrixForm

Raise the final p matrix to some large power to verify it is regular

MatrixPower[p,50]//N//MatrixForm