4.12.2 Problem 8.5

Part(a)

The following is the Hastings-Metrpolois algorithm implementation.

This algorithm generates a time-reversible M.C. (referred to as \(p\) in the lecture notes) given an irreducible M.C. (called \(q\) or the original chain) and given a stationary distribution \(\pi \) for that chain.

Input: \(f\left ( x\right ) \) defined over the states \(x,\) and \(edge\left ( x\right ) \) which represents the number of edges connected to \(x\)

  1. For each state \(x\) calculate \(\pi \left ( x\right ) =\frac {f\left ( x\right ) }{{\displaystyle \sum \limits _{v\in V}} f\left ( v\right ) }\) and for each state \(x\) calculate \(edge(x)\)
  2. compute \(q\left ( x,y\right ) =\frac {1}{edge\left ( x\right ) }\) whenever \(edge\left ( x\right ) \neq 0\) else set \(q\left ( x,y\right ) =0\)
  3. Select a state \(x\) by random to start from.
  4. Let \(n=1\) and let \(X_{1}=x\)
  5. Let \(S\) be the set of all states that can be reached in one step from \(x\). These will be the states \(y\) in which \(q\left ( x,y\right ) \neq 0\)
  6. Select a state \(y\) from \(S\) by random (using a uniform \(U\left [ 0,1\right ] \) random number generator)
  7. Calculate \(\beta \left ( x,y\right ) =\min \left \{ 1,\frac {\pi \left ( y\right ) q\left ( y,x\right ) }{\pi \left ( x\right ) q\left ( x,y\right ) }\right \} \)
  8. Generate a random number \(u\) from \(U\left [ 0,1\right ] \)
  9. Let \(n=n+1\)
  10. Compare \(u\) to \(\beta \left ( x,y\right ) \).
  11. IF \(u<\beta \left ( x,y\right ) \) THEN \(X_{n}=y\) (select the new state) ELSE \(X_{n}=X_{n-1}\) (stay in same state) ENDIF
  12. Let \(x=X_{n}\)
  13. If \(n>\) some Max number of iterations or if we reached some convergence limit Then go to 15
  14. GOTO 5
  15. Algorithm is complete. Now generate the time reversible MC as follows

    1. Scan the state path generate \(X_{n}\,\ \)and count how many times state \(x\) switches to state \(y\) in one step
    2. Do the above for all the states \(x\)
    3. Divide the above number by the total number of steps made to generate \(p\left ( x,y\right ) \)

Since the problem now asks to implement Hastings-Metropolis, then I used the data given at the end of the problem and implemented the above simulation using that data5. Please see appendix for code and final \(P\) matrix generated.

Part (a1)

This is similar the problem 8.4 part(I). To show that the \(p\) (final M.C.) is irreducible, we need to show that there exist no closed proper subsets. Since the graph \(G\) is connected, then we just need to show whenever there is an edge between vertex \(x\) and \(y\) then there corresponds in the chain representation of the final \(p\) matrix a non-zero \(p\left ( x,y\right ) \) and also a non-zero \(p\left ( y,x\right ) \). This will insure that the each state can transition to each other state, just as each vertex can be visited from each other vertex (since it is a connected graph).

Let us consider any 2 vertices say \(x,y\) with a direct edge between them (this is the only case we need to consider due to the argument above). We need to show the resulting \(p\left ( x,y\right ) \) and \(p\left ( y,x\right ) \) are non-zero

Consider \(p\left ( x,y\right ) \) first. Since

\begin{align*} p\left ( x,y\right ) & =q\left ( x,y\right ) \beta \left ( x,y\right ) \\ & =\frac {1}{edge\left ( x\right ) }\min \left \{ 1,\frac {\pi \left ( y\right ) q\left ( y,x\right ) }{\pi \left ( x\right ) q\left ( x,y\right ) }\right \} \\ & =\frac {1}{edge\left ( x\right ) }\min \left \{ 1,\frac {\frac {f\left ( y\right ) }{{\displaystyle \sum \limits _{v\in V}} f\left ( v\right ) }\frac {1}{edge\left ( y\right ) }}{\frac {f\left ( x\right ) }{{\displaystyle \sum \limits _{v\in V}} f\left ( v\right ) }\frac {1}{edge\left ( x\right ) }}\right \} \end{align*}

Hence

\begin{equation} \fbox {$p\left ( x,y\right ) =\frac {1}{edge\left ( x\right ) }\min \left \{ 1,\frac {f\left ( y\right ) edge\left ( x\right ) }{f\left ( x\right ) edge\left ( y\right ) }\right \} $} \tag {1}\end{equation}

Then it is clear that whenever there is an edge between \(x,y\) then \(p\left ( x,y\right ) \neq 0\) since both \(f\left ( x\right ) \) and \(f\left ( y\right ) \) are positive (not zero) and also edge(x) and edge(y) are non-zero as well. Hence we see that \(p\left ( x,y\right ) \neq 0\). Similar argument shows that \(p\left ( y,x\right ) \neq 0\).

This shows that M.C. represented by \(P\) is irreducible.

Part (a2)

The condition for regular chain \(P\) is that there exist at least one state \(x\) such that \(p\left ( x,x\right ) >0.\)From (1) above we can decide under what conditions this will occur.

Consider a vertex \(x\) with \(edge\left ( x\right ) \) edges from it connected to vertices \(y_{1},y_{2},\cdots ,y_{r}\). Then from (1) we see that

\begin{align*} p\left ( x,y_{i}\right ) & =\frac {1}{edge\left ( x\right ) }\min \left \{ 1,\frac {f\left ( y_{i}\right ) edge\left ( x\right ) }{f\left ( x\right ) edge\left ( y_{i}\right ) }\right \} \\ & =\frac {1}{edge\left ( x\right ) }\min \left \{ 1,\frac {\frac {f\left ( y_{i}\right ) }{edge\left ( y_{i}\right ) }}{\frac {f\left ( x\right ) }{edge\left ( x\right ) }}\right \} \end{align*}

The condition for having \(p\left ( x,x\right ) >0\) is that \(\min \left \{ 1,\frac {\frac {f\left ( y_{i}\right ) }{edge\left ( y_{i}\right ) }}{\frac {f\left ( x\right ) }{edge\left ( x\right ) }}\right \} <1\), since this will cause \(p\left ( x,y_{i}\right ) \) to be some quantity less than \(\frac {1}{r}\) and so when summing over all \(r\) there will be a deficit in the sum and we have to compensate for it to make it \(1\) by adding \(p\left ( x,x\right ) \). But for \(\min \left \{ 1,\frac {\frac {f\left ( y_{i}\right ) }{edge\left ( y_{i}\right ) }}{\frac {f\left ( x\right ) }{edge\left ( x\right ) }}\right \} \) to be less than ONE means that \(\frac {f\left ( y_{i}\right ) }{edge\left ( y_{i}\right ) }<\frac {f\left ( x\right ) }{edge\left ( x\right ) }\)

Hence the condition for finding an Aperiodic state is finding a vertex \(x\) such that the above holds for one of the vertices \(y_{i}\) this vertex is directly connected to. For example, if \(y_{i}\) had the same number of edges from it as does \(x\), then the condition will be that \(f\left ( y_{i}\right ) <f\left ( x\right ) \). And if \(y_{i}\) has less or more edges from it than \(x\) has, then we need the ratio \(\frac {f\left ( y_{i}\right ) }{edge\left ( y_{i}\right ) }\) to be less than \(\frac {f\left ( x\right ) }{edge\left ( x\right ) }\).

The above is the same as saying \(\frac {f\left ( x\right ) }{edge\left (x\right )}\) must be constant for the \(p\) not to be regular.

Part(A3)

Since the chain is irreducible, then there is a reverse Markov chain (proof is on page 8.1 and 8.2 of lecture notes). Hence for an irreducible chain the balance equations hold

\begin{equation} r\left ( x,y\right ) =\frac {\pi \left ( y\right ) p\left ( y,x\right ) }{\pi \left ( x\right ) } \tag {2}\end{equation}

Now if the chain the time reversible as well, then \(r\left ( x,y\right ) =p\left ( x,y\right ) \), Then the balance equation (1) becomes

\begin{align} \pi \left ( x\right ) p\left ( x,y\right ) & =\pi \left ( y\right ) p\left ( y,x\right ) \nonumber \\ \frac {f\left ( x\right ) }{{\displaystyle \sum \limits _{v\in V}} f\left ( v\right ) }q\left ( x,y\right ) \beta \left ( x,y\right ) & =\frac {f\left ( y\right ) }{{\displaystyle \sum \limits _{v\in V}} f\left ( v\right ) }q\left ( y,x\right ) \beta \left ( y,x\right ) \nonumber \\ \frac {f\left ( x\right ) }{{\displaystyle \sum \limits _{v\in V}} f\left ( v\right ) }\frac {1}{edge\left ( x\right ) }\left ( \min \left \{ 1,\frac {f\left ( y\right ) edge\left ( x\right ) }{f\left ( x\right ) edge\left ( y\right ) }\right \} \right ) & =\frac {f\left ( y\right ) }{{\displaystyle \sum \limits _{v\in V}} f\left ( v\right ) }\frac {1}{edge\left ( y\right ) }\left ( \min \left \{ 1,\frac {f\left ( x\right ) edge\left ( y\right ) }{f\left ( y\right ) edge\left ( x\right ) }\right \} \right ) \nonumber \\ \frac {f\left ( x\right ) }{{\displaystyle \sum \limits _{v\in V}} f\left ( v\right ) }\frac {1}{edge\left ( x\right ) }\left ( \min \left \{ 1,\frac {\frac {f\left ( y\right ) }{edge\left ( y\right ) }}{\frac {f\left ( x\right ) }{edge\left ( x\right ) }}\right \} \right ) & =\frac {f\left ( y\right ) }{{\displaystyle \sum \limits _{v\in V}} f\left ( v\right ) }\frac {1}{edge\left ( y\right ) }\left ( \min \left \{ 1,\frac {\frac {f\left ( x\right ) }{edge\left ( x\right ) }}{\frac {f\left ( y\right ) }{edge\left ( y\right ) }}\right \} \right ) \tag {3}\end{align}

Hence we need to show that the equation (3) above holds to show the chain is time reversible.

There are 3 cases to consider:

  1. \(\frac {f\left ( y\right ) }{edge\left ( y\right ) }=\frac {f\left ( x\right ) }{edge\left ( x\right ) }\)
  2. \(\frac {f\left ( y\right ) }{edge\left ( y\right ) }<\frac {f\left ( x\right ) }{edge\left ( x\right ) }\)
  3. \(\frac {f\left ( y\right ) }{edge\left ( y\right ) }>\frac {f\left ( x\right ) }{edge\left ( x\right ) }\)

For case (1), LHS of equation (3) simplifies to \(\frac {f\left ( x\right ) }{{\displaystyle \sum \limits _{v\in V}} f\left ( v\right ) }\frac {1}{edge\left ( x\right ) }\) and the RHS of (3) simplifies to \(\frac {f\left ( y\right ) }{{\displaystyle \sum \limits _{v\in V}} f\left ( v\right ) }\frac {1}{edge\left ( y\right ) }\), but since \(\frac {f\left ( y\right ) }{edge\left ( y\right ) }=\frac {f\left ( x\right ) }{edge\left ( x\right ) }\), then LHS=RHS.

Hence balance equation (3) is satisfied for case (1).

For case(2), LHS of (3) simplifies \(\frac {f\left ( x\right ) }{{\displaystyle \sum \limits _{v\in V}} f\left ( v\right ) }\frac {1}{edge\left ( x\right ) }\left ( \frac {\frac {f\left ( y\right ) }{edge\left ( y\right ) }}{\frac {f\left ( x\right ) }{edge\left ( x\right ) }}\right ) =\frac {f\left ( y\right ) }{{\displaystyle \sum \limits _{v\in V}} f\left ( v\right ) }\frac {1}{edge\left ( y\right ) }\) and RHS of (3) simplifies to \(\frac {f\left ( y\right ) }{{\displaystyle \sum \limits _{v\in V}} f\left ( v\right ) }\frac {1}{edge\left ( y\right ) },\)then LHS=RHS.

Hence balance equation (3) is satisfied for case (2)

.

For case (3), LHS of (3) simplifies \(\frac {f\left ( x\right ) }{{\displaystyle \sum \limits _{v\in V}} f\left ( v\right ) }\frac {1}{edge\left ( x\right ) }\) and RHS of (3) simplifies to \(\frac {f\left ( y\right ) }{{\displaystyle \sum \limits _{v\in V}} f\left ( v\right ) }\frac {1}{edge\left ( y\right ) }\left ( \frac {\frac {f\left ( x\right ) }{edge\left ( x\right ) }}{\frac {f\left ( y\right ) }{edge\left ( y\right ) }}\right ) =\frac {f\left ( x\right ) }{{\displaystyle \sum \limits _{v\in V}} f\left ( v\right ) }\frac {1}{edge\left ( x\right ) },\)then LHS=RHS.

Hence balance equation (3) is satisfied for case (3)

.

Hence in all 3 cases we showed the balance equation is satisfied.

Hence M.C. is time reversible.

Part(b)

A small program written to construct the \(P\) matrix directly following instructions on page 8.4 of lecture notes. The following is the resulting \(P\) matrix

Now to check that the final chain \(P\) is regular, it was raised to some high power to check that all entries in the \(P^{m}>0\). This is the result

The above verifies that the final matrix \(p\) is regular.

Using the Hastings-Metropolis simulation algorithm, the convergence to the above matrix was slow. Had to make 2 million observation to be within 3 decimal points from the above. Here is the \(P\) matrix generated from Hastings algorithm for \(N=2,000,000\)