My question concerns Markov Chain Monte Carlo methods based on the Metropolis algorithm.

Suppose we have a system with $N$ levels, and at each instant in time $n$, the probability distribution among the levels is given by $\vec{\pi}^{(n)}=\langle\pi^{(n)}_0, \pi^{(n)}_1, ..., \pi^{(n)}_N\rangle$. Let the transition probability matrix be $\mathbb{P}$. We can obtain the probability distribution at the next time step as

$$\vec{\pi}^{(n+1)}=\mathbb{P}\ \vec{\pi}^{(n)}.$$

Let's write this condition in the following equivalent way, noting that we are considering each component now:

$$\begin{aligned} \pi^{(n+1)}(x)&=\pi^{(n)}(x)\mathbb{P}(x\rightarrow x)+\sum_{y\ \neq\ x}\pi^{(n)}(y)\mathbb{P}(y\rightarrow x) \\ &=\pi^{(n)}(x)\left[1-\sum_{y\ \neq\ x} \mathbb{P}(x\rightarrow y)\right]+\sum_{y\ \neq\ x}\pi^{(n)}(y)\mathbb{P}(y\rightarrow x) \\ \end{aligned}$$

$$\tag{1} \Rightarrow \quad \pi^{(n+1)}(x)-\pi^{(n)}(x)= \sum_{y\ \neq\ x}\pi^{(n)}(y)\mathbb{P}(y\rightarrow x)-\pi^{(n)}(x)\mathbb{P}(x\rightarrow y).$$

If we seek the steady state solution, we seek to make the right hand side of Eq (1) equal to zero. One such way is to require

$$\tag{2} \pi^{(n)}(y)\mathbb{P}(y\rightarrow x)=\pi^{(n)}(x)\mathbb{P}(x\rightarrow y),$$

a condition called the Detailed Balance equation. If Eq (2) is satisfied, we are guaranteed to have a steady state probability distribution from our Markov chain process. The Metropolis algorithm

$$\tag{3} \mathbb{P}(x\rightarrow y) = \mathrm{min}\left[1,\frac{\pi(y)}{\pi(x)}\right],$$

will satisfy Eq (2) (here, I am letting the apriori to be 1).

My question

As a physicist, I'll illustrate my confusion through a physical example. Say we have a system with three levels of energy $\{\epsilon_i\}_{i=0,1,2}$. We already know from Boltzmann statistics that the probability for each level, up to the normalization, which doesn't matter here, is

$$\pi_i=\mathrm{e}^{-\beta \epsilon_i}.$$

We use this in the Metropolis algorithm in Eq (3) to advance with our Markov chain.

  1. What probability distribution are we sampling at equilibrium? Is it $\pi_i$? If not, then what?

It would be redundant if it's $\pi_i$, no? We already have a mathematical expression for that PDF and we could just use simple towering or rejection method to sample it.

I'm confused because it seems Markov chains are used to sample the PDF $\pi$, while the framework itself relies on $\pi$ having a mathematical expression for the Metropolis part. But, if we have an expression for the PDF, then why bother with Markov chains for sampling?

  • $\begingroup$ Knowing the distribution up to an unknown scaling constant is different from knowing the distribution. The Metropolis-Hastings sampler is helpful when you know the distribution except for an unknown scaling constant and can't evaluate that constant because it would be impractical to integrate over the typically very high dimensional probability distribution. $\endgroup$ – Brian Borchers Sep 14 '19 at 4:47
  • $\begingroup$ Can't you just do Monte Carlo integration? I learned the MCMC method as a very powerful tool, and I must be missing something here because given my example it doesn't seem to be doing much! So I'm trying to find what I'm missing. $\endgroup$ – Ptheguy Sep 14 '19 at 15:56

Your Answer

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

Browse other questions tagged or ask your own question.