I wanted to create a simple example of the Metropolis-Hastings algorithm as a sanity check but the check failed and the results do not line up. I have a vector of energies $\vec E=(E_1,E_2,E_3)$ and the expected result is that the stationary distribution follows the Boltzmann distribution $$p(i)=\frac{e^{-\beta E_i}}{\sum_j e^{-\beta E_j}}$$

Here is my code:

e = {1, 2, 3};
beta = 1.;
P = Table[Min[1, Exp[-beta (e[[j]] - e[[i]])]], {i, 1, 3}, {j, 1, 3}]
Do[P[[i]] /= Total@P[[i]], {i, 1, 3}];
M = DiscreteMarkovProcess[pinit, P];

stationaryDistr = PDF[StationaryDistribution[M], Range[3]]

Exp[-beta e]/(Total@Exp[-beta e])

which gives the output $ \{0.540, 0.313, 0.146 \}$ for stationaryDistr and $\{0.665,0.245,0.090\}$ for the last expression. Does anyone know how I can get the expected result?

  • 1
    $\begingroup$ The problem is that when you divide the rows of P by their sums to make it a valid transition matrix, you also violate the detailed balance condition $p(i \to j)/p(j \to i) = e^{-\beta(E_i - E_j)}$ (since you're dividing by a different number on each row.) gwr's method correctly constructs the transition matrix. $\endgroup$ May 3 at 17:23

1 Answer 1


Disclaimer: I am not a physicist... :)

... but I know how Metropolis-Hastings works. So let's got about this with the following idea:

  1. We will use a simple proposal distribution that selects one of the other two states ($j \neq i$) with a probability of $0.5$.

  2. If $E(j) \leq E(i)$ we will select the proposed state with probability of $1$, i.e., we will switch to a proposed state with lower or equal energy level with certainty.

  3. If $E(j) > E(i)$, i.e., the energy level of the proposed state is higher, we will select the proposed state with probability $e^{-\Delta E}$.

A note with regard to the proposal distribution: In the Metropolis-Hastings algorithm the decision criterion is actually:

\begin{align} \alpha &= \frac{g(\theta')/q(\theta'|\theta_{i-1})}{g(\theta_{i-1})/q(\theta_{i-1})|\theta')} \end{align}

where $g(\cdot)$ is the PDF for a distribution that is proportional to the distribution that you want to sample from (here the Boltzmann distribution), $q(\cdot)$ is the PDF for the proposal distribution, $\theta'$ is a candidate value for sampling (here: state), and $\theta_{i-1}$ is the previous state or sampled value. It reduces to the simple form

\begin{align} \alpha &= \frac{g(\theta')}{g(\theta_{i-1})} = \frac{e^{-\beta E_j}}{e^{-\beta E_i}} \end{align}

if the proposal distribution is symmetric, as Michael Seifert has noted, i.e., if $q(\theta'|\theta_{i-1}) = q(\theta_{i-1}|\theta')$, which is the case for the simple proposal function here where $q(j|i) = q(i|j) = 0.5$.

Doing it by hand

Let's quickly do this by hand. First, we can come up with the transition probabilities for the other states, i.e., the diagonal will be calculated in the second step.

enter image description here

Then we can replace the $e_i$ in the diagonal with the probabilities for not switching to another energy level from the fact that each row needs to add up to one.


The transition matrix in numerical form looks like this:

matP // MatrixForm

numerical matrix

We can verify this nicely for the third row: Since each of the other states $E_1,E2$ has a lower energy level, the probability to switch to these states is simply the probability for this state being proposed (no state by itself cannot be a candidate) and thus $0.5$ in each case. Thus being in the third state means that there will be a switch with certainty and the probability for remaining in that state $P_{3,3}$ will be $0$.

If you go about finding the stationary distribution as before you will get your result:

M = DiscreteMarkovProcess[pinit, matP];
stationaryDistr = PDF[StationaryDistribution[M], Range[3]]
(* {0.665241, 0.244728, 0.0900306} *)

Doing it programmatically

pinit = {0, 1, 0};
e = {1, 2, 3};
beta = 1.;
matP = Table[ 
   Boole[ i != j ] 0.5 Min[1., Exp[ - beta (e[[j]] - e[[i]])]] 
   + Boole[ i == j] 0,
   {i, 3}, {j, 3}
Do[ matP[[i]] = matP[[i]] /.  x_ /; x == 0 -> 1 - Total[ matP[[i]] ], {i, 3} ];
M = DiscreteMarkovProcess[pinit, matP];
stationaryDistr = PDF[StationaryDistribution[M], Range[3]]
(* {0.665241, 0.244728, 0.0900306} *)


We can make this solution a bit more general:

transitionMatrix = Function[ {e, b},
         n = Length @ e,
         matrix, row
      matrix = Table[
         Boole[i != j] 1/(n - 1) Min[1., Exp[-b (e[[j]] - e[[i]])]]
         + Boole[i == j] 0.
         {i, n}, {j, n}
         Function[ row, 
            row /. p_ /; p == 0 -> 1 - Total @ row

pint = PadRight[ {1}, 10 ];
M = DiscreteMarkovProcess[ pint, transitionMatrix[ Range[10], 1 ] ];
PDF[ StationaryDistribution[M], Range[10] ]
   0.00425939,0.00156694,0.000576446,0.000212062,0.0000780134} *)
  • $\begingroup$ This is great! I implemented a proper MC system before and I felt stupid that I couldn't get this simple example to work. $\endgroup$ May 3 at 17:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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