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<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
 A: 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 ****)
```

 

A: 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.
