# Metropolis-Hastings and underlying Markov process

I tried to understand the workings of the Metropolis-Hasting algorithm. As far as I can understand, it allows to draw samples from an unknown distribution $$T(x)$$ as long as a function proportional to it, $$F(x)$$ is known. In practice, it seems a great advantage is not to have to calculate the partition function $$Z$$.

The algorithm constructs a Markov process with transition probabilities $$P_{ij}$$ which is then made evolving to its stationary distribution $$\pi$$, which uniquely exists under some conditions.

I would like to explicitly compute the transition matrix $$P$$ and find the stationary distribution solving a linear system. The transition probability for Metropolis-Hastings is split in two steps,$$P_{ij} = A_{ij} g_{ij}$$, getting a candidate via distribution $$g$$ and accepting/rejecting it according to $$A$$. Getting a candidate is done via a distribution, to be chosen, $$g(j | i)$$, giving the conditional probability of picking element $$j$$ if currently at $$i$$. The acceptance $$A_{ij}$$ is defined as $$A_{ij} = \min \Big( 1, \frac{F_j}{F_i)} \frac{g(i | j)}{g(j | i)}\Big)$$ As said $$P_{ij} = A_{ij} g_{ij}$$ and this is the transition matrix I would like to use, for example to calculate the stationary distribution $$\pi$$ using the equation $$\pi P = \pi$$. I tried to apply this but I am not getting it to work. I present the simplest example. I have a particle that can be in one of four positions $$1,2,3,4$$,

The energy level of each position equals $$0,1,2,-1$$ respectively. In the canonical setting, the stationary distribution should be, taking the inverse temperature $$\beta = 1$$ for simplicity $$F_i = \frac{\exp\big(- \beta E_i\big)}{Z}$$ and $$Z = \exp\big(0 \big) + \exp\big(- 1 \big) + \exp\big(-2 \big) + \exp\big(1 \big)$$ As far as $$g$$ is concerned, I assume that only nearest-neighbour (imagining the $$1,2,3,4$$ configurations are located on a ring) jumps as well as jumps to the current position are possible, and they are all equiprobable. Then $$g$$ could be represented as a matrix $$\begin{equation*} g_{ij} = \begin{pmatrix} 1/3 & 1/3 & 0 & 1/3 \\ 1/3 & 1/3 & 1/3 & 0 \\ 0 & 1/3 & 1/3 & 1/3 \\ 1/3 & 0 & 1/3 & 1/3 \end{pmatrix} \end{equation*}$$

It seems I am already going wrong here, as looking at the formula for $$A$$ one needs $$g$$ never to $$0$$, but surely there must some subtlety I am missing here, as I have seen Metropis-Hasting applications using nearest-neighbour jumps only, at least I think so...

Coming to the $$A_ij$$, it seems to me that it is either $$1$$, if $$E_j, else equal to $$\exp(E_j-E_i)$$.

Putting it all together, the transition probability matrix should be

$$\begin{equation*} P = \begin{pmatrix} 1-\exp(-1)/3- \exp(1)/3& \exp(-1)/3 & 0 & \exp(1)/3 \\ 1/3 & 2/3 - \exp(-1)/3 & \exp(-1)/3 & 0 \\ 0 & 1/3 & 2/3 - \exp(-1)/3 & \exp(-1)/3 \\ \exp(1)/3 & 0 & \exp(-3)/3 & 1-\exp(1)/3-\exp(-3)/3 \end{pmatrix} \end{equation*}$$

But it does not satisfy the equation $$\pi P = \pi$$ that is should if $$\pi$$ is to be the stationary distribution.

What is that I could be doing wrong?

What is also bugging me as mad is that the stationary distribution according to the exact canonical calculation should equal

[0.23688282, 0.08714432, 0.0320586 , 0.64391426]

while computing $$\pi P$$ I get

[0.10120567, 0.08714432, 0.03881357, 0.77283644]

the second element being identical!

I would be glad to share a simple Python snippet that performs the calculation, if anybody were interested at all, thanks

• There is some imprecision in notation, as you denote states sometimes by $i,j$ and sometimes by $x_i,x_j$. The point of MH is that the desired distribution is the equilibrium distribution of the Markov chain, which has therefore to be made reversible (it cannot get stuck in a state) - this is usually imposing the detailed balance condition, $W_{ij}p_i=W_{ji}p_j$, which is implicit in the acceptance probability. I am posting this as a comment, since it is hard to follow through the formulae in the question. It could be better, if you focus on conceptual points. Commented May 23, 2022 at 14:16
• @RogerVadim, thanks for this, well I though $I,j$ represent states, and $x_i, x)j$ represent their coordinate that is relevant for the energy, I will improve it. For the second part, I am not sure I see where you are tending. Detailed Balance as you say is implicit in the acceptance probability, that is properly how Metropolis is built. I want to build a Markov chain with the same stationary distribution. MH is sampling a Markov process, I want the proper transition matrix, I cannot figure out why it is not working (apart from the dubious choice of $g$ maybe, zero somewhere). Commented May 23, 2022 at 14:25
• @RogerVadim, the fact is I thought I had a good grasp of the concepts, until I put it to the test with an example. Would you please tell which formulae are not clear, for me to improve their clarity? I must be missing something but I think working out a trivial concrete example must be the best way to get rid of the misunderstanding. From the abstract point of view, I went numerous times through the derivation I feel I got it fully, but it is clearly not so if I cannot work out a 4x4 stochastic matrix. If you let me know where to improve the question, I would be delighted to do so. Commented May 23, 2022 at 14:44
• There seems to be quite a few typos in your equations: the energies in the partition function are different from those cited in the sentence above, and $g_{ij}$ has zeros in different places in the last equation and the one before the last... or is the last one for $A_{ij}$? Anyhow, the first $g_{ij}$ seems to have has wrong positions of zeros for detailed balance. Commented May 23, 2022 at 14:52
• @Smerdjakov Classic! Hope that everything goes smoothly now :)
– Javi
Commented May 23, 2022 at 20:37

I will answer as I found the mistake, my apologies for having been careless in posting I can assure you I spent much time doing hand calculations but sometimes things just evade your attention.

The entry $$P_{14}$$ in the probability matrix is wrong, should equal $$1/3$$ as the energy level $$E_1$$ is large than $$E_4$$. The method works, which is of course not surprising as it is based on MCMC definitions.

For completeness I attach the code that I put together and a bug of which started the whole question. It is not corrected and one can check in the last line that indeed the stationary distribution calculated via the canonical distribution, $$p_0$$ equals the stationary distribution $$\pi$$ of the Markov chain, satisfying $$\pi P = \pi$$

import numpy as np
#from scipy import linalg
energyLevels = np.array([0.,1.,2.,-1.],dtype=np.double)
BoltzFactors = np.exp(-energyLevels,dtype=np.double)
Z = np.sum(BoltzFactors,dtype=np.double)
# Canonical Distribution
p0 = np.exp(-energyLevels) / Z

# this dummy case assumes there are 3 transitions possible per vertex
ProbTrans = np.zeros((4,4))
indeces = np.argwhere(aux!=0)
for i in range(indeces.shape[0]):
np.minimum (  1, np.exp(-(energyLevels[indeces[i][1]] - energyLevels[indeces[i][0]]))))
ProbTrans = ProbTrans
ProbTrans[i,i] = 1 - np.sum(ProbTrans[i,:])
return ProbTrans
#Check if Markov stationary distribution equals p0
np.allclose(np.matmul(np.transpose(PP), np.transpose(p0) ),p0)
(*** TRUE ****)


`

Jumps between energy levels are not equally probabilistic, your transition matrix is just too arbitrary.
First according to your definition $$g_{ij}$$ means that the particle was previously at energy level i and then transitioned to energy level j.Then the equilibrium distribution is set to $$F_i = e^{-E_i}$$, which does not need to be normalized.
At this time, the transformation matrix needs to meet two conditions: The first is $$\sum_{j=1}^4 g_{ij} = 1$$ It means that the particle's previous position is deterministic, whether it stays in place or transitions to another energy level next.
The second is: $$\sum_{j=1}^4 g_{ij} e^{-E_j} = e^{-E_i}$$ It means that the delicate balance condition has been met at this time. Your transition matrix needs to be solved under these two conditions, I totally don't understand how you can arbitrarily write a transition matrix ignoring the second condition, which doesn't even satisfy this condition at all.

• Thanks for your answer. Both those conditions are enforced, the first on $g$ by my closest-neighbour jump assumption, and the second one implicitly by the Metropolis Hasting algorithm. I calculate the elements of the transitions matrix using $P_{ij} = A_{ij} g{ij}$ and the $A_{ij}$ ar calculated using the Metropis Hasting equation, first equation at the beginning of the post. the methods works as it must, I was messing it up with a silly mistake, sorry for the confusion. Commented May 25, 2022 at 13:01