5
$\begingroup$

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$, toy canonical Markov Chain

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<E_i$, 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

$\endgroup$
11
  • $\begingroup$ 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. $\endgroup$
    – Roger V.
    Commented May 23, 2022 at 14:16
  • $\begingroup$ @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). $\endgroup$
    – Smerdjakov
    Commented May 23, 2022 at 14:25
  • $\begingroup$ @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. $\endgroup$
    – Smerdjakov
    Commented May 23, 2022 at 14:44
  • $\begingroup$ 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. $\endgroup$
    – Roger V.
    Commented May 23, 2022 at 14:52
  • 1
    $\begingroup$ @Smerdjakov Classic! Hope that everything goes smoothly now :) $\endgroup$
    – Javi
    Commented May 23, 2022 at 20:37

2 Answers 2

3
$\begingroup$

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
AdjMat = np.array([[1,1,0,1],[1,1,1,0],[0,1,1,1],[1,0,1,1]])

# Assembly Transition Probability Matrix given adjacency matrix Adjmat, 
# this dummy case assumes there are 3 transitions possible per vertex
def MarkovTrans(AdjMat, energyLevels):
    ProbTrans = np.zeros((4,4))
    aux = AdjMat - np.diag(np.ones(4))
    indeces = np.argwhere(aux!=0)
    for i in range(indeces.shape[0]):
        ProbTrans[indeces[i][0],indeces[i][1]] = ( (1/3)*AdjMat[indeces[i][0],indeces[i][1]]*
            np.minimum (  1, np.exp(-(energyLevels[indeces[i][1]] - energyLevels[indeces[i][0]]))))
        ProbTrans = ProbTrans 
    for i in range(AdjMat.shape[0]):
        ProbTrans[i,i] = 1 - np.sum(ProbTrans[i,:])
    return ProbTrans
PP = MarkovTrans(AdjMat, energyLevels)
#Check if Markov stationary distribution equals p0
np.allclose(np.matmul(np.transpose(PP), np.transpose(p0) ),p0)
(*** TRUE ****)
```

 
$\endgroup$
0
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – Smerdjakov
    Commented May 25, 2022 at 13:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.