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.
- 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?
