Markov Chain Monte Carlo on Finite State Spaces

TL;DR

This paper provides a clear, mathematically rigorous explanation of MCMC methods on finite state spaces, including convergence proofs and fundamental algorithms like Gibbs and Metropolis-Hastings.

Contribution

It offers a simple, coherent presentation of MCMC theory with proofs of key theorems and discussion of core sampling methods, accessible with basic mathematical knowledge.

Findings

Proves a convergence theorem for finite Markov chains

Introduces a minimal Perron-Frobenius theorem version

Discusses fundamental MCMC algorithms: Gibbs and Metropolis-Hastings

Abstract

We elaborate the idea behind Markov chain Monte Carlo (MCMC) methods in a mathematically coherent, yet simple and understandable way. To this end, we proof a pivotal convergence theorem for finite Markov chains and a minimal version of the Perron-Frobenius theorem. Subsequently, we briefly discuss two fundamental MCMC methods, the Gibbs and Metropolis-Hastings sampler. Only very basic knowledge about matrices, convergence of real sequences and probability theory is required.

Equations

$$\mathbb{P}\left(nN<\mathbf{T}\right)=\mathbb{P}\left({\mathbf{X}}_i\ne {\mathbf{Z}}_i\text{ for all }i\le nN\right)\le \mathbb{P}\left({\mathbf{X}}_{i\cdot N}\ne {\mathbf{Z}}_{i\cdot N}\text{ for all }i\le n\right)$$

$$=\sum _{z_0\ne s_0}{\pi }_{z_0}\sum _{s_1\in \mathcal{S}}{\kappa }_{s_0{s}_1}^N\sum _{z_1\ne s_1}{\kappa }_{z_0{z}_1}^N...\sum _{s_n\in \mathcal{S}}{\kappa }_{s_{n-1}{s}_n}^N\underset{\le 1-ϵ}{\underbrace{\sum _{z_n\ne s_n}{\kappa }_{z_{n-1}{z}_n}^N}}\le (1-ϵ{)}^n\stackrel{n\to \infty }{⟶}0$$

$$\mathbb{P}\left({\mathbf{Y}}_{0:n}=s_{0:n}\right)=\mathbb{P}\left({\mathbf{Y}}_{0:n}=s_{0:n},n<\mathbf{T}\right)+\sum _{\ell =0}^n\mathbb{P}\left({\mathbf{Y}}_{0:n}=s_{0:n},\mathbf{T}=\ell \right)$$

$$=\text{mathantt}{p}_{s_{1:n}}\mathbb{P}\left({\mathbf{Z}}_{0:n}\ne s_{0:n}\right)+\sum _{\ell =0}^n\text{mathantt}{p}_{s_{1:\ell }}\mathbb{P}\left({\mathbf{Z}}_{\ell }=s_{\ell },{\mathbf{Z}}_{0:\ell -1}\ne s_{0:\ell -1}\right)\text{mathantt}{p}_{s_{\ell +1:n}}=\text{mathantt}{p}_{s_{1:n}}$$

$$({\kappa }^n{)}_{s_{0}s}$$

$${\pi }_s$$

$$\sum _{z\in N}{x}_z{\kappa }_{zs}+\sum _{z\in \bar{N}}{x}_z{\kappa }_{zs}=x_s\Rightarrow \sum _{z\in N}{x}_z\sum _{s\in \bar{N}}{\kappa }_{zs}+\sum _{z\in \bar{N}}{x}_z\sum _{s\in \bar{N}}{\kappa }_{zs}=\sum _{s\in \bar{N}}{x}_s$$

$$\displaystyle\Leftrightarrow\ \$$

$${\kappa }_{zs}^j=\frac{\text{scalebox}1.2\text{mathds}1}{\text{scalebox}1.2\text{mathds}1}$$

$$\text{mathantt}{a}_{zs}=\min \left\{1,\frac{{\pi }_s{\rho }_{sz}}{{\pi }_z{\rho }_{zs}}\right\}$$

$$\displaystyle\pi_{z}\rho_{zs}\mathantt{a}_{zs}=\min\Big{\{}\pi_{z}\rho_{zs},\pi_{s}\rho_{sz}\Big{\}}=\pi_{s}\rho_{sz}\mathantt{a}_{sz}$$

Full text

Keywords: Metropolis-Hastings, Gibbs sampling, Convergence of Markov chains, Perron-Frobenius theorem

1 Introduction

MCMC techniques aim at drawing samples from prespecified distributions. They do so in an indirect, approximate fashion through Markov chains. This is important since the distributions deployed in practice are often too complex to be dealt with directly or even unavailable in closed form.

There exists a tremendous number of scientific articles and books about MCMC. See for example Bishop and Mitchell (2014) for a vivid and more comprehensive introduction without mathematical proofs.

Let $\pi$ be a distribution over a finite state space $\mathcal{S}$ and $\pi_{s}\in[0,1]$ the probability of state $s\in\mathcal{S}$ under $\pi$, whereby $\sum_{s\in\mathcal{S}}\pi_{s}=1$. By virtue of the strong law of large numbers, independent samples from $\pi$ (so-called i.i.d. samples) can be used universally to approximate expectations w.r.t. $\pi$. Thus, for a set of such samples, say $x_{1},...,x_{m}\in\mathcal{S}$, and an arbitrary function $\mathantt{f}:\mathcal{S}\rightarrow\mathbb{R}^{\ell}$, we get that $\frac{1}{m}\sum_{i=1}^{m}\mathantt{f}(x_{i})\approx\sum_{s\in\mathcal{S}}\pi_{s}\mathantt{f}(s)$. This simple recipe poses one of the most powerful tools in statistics.

An example for $\mathantt{f}$ is the indicator function \raisebox{-0.85358pt}{\scalebox{1.2}{\mathds{1}}}\{s\in A\} for events $A\subseteq\mathcal{S}$. It is one if the condition in the brackets is true and zero otherwise. Its expectation yields the probability of $A$.

A Markov chain over $\mathcal{S}$ is defined through an (arbitrary) initial state $s_{0}\in\mathcal{S}$ and a transition kernel $\kappa$. $\kappa$ is a non-negative function over $\mathcal{S}\times\mathcal{S}$ such that $\sum_{s\in\mathcal{S}}\kappa_{zs}=1$ for all $z\in\mathcal{S}$. $\kappa_{z}$ can be interpreted as a conditional distribution.

The Markov chain starts in state $s_{0}$ and evolves according to $\kappa$ in an iterative fashion: The distribution of the first link in the chain is $\kappa_{s_{0}}$ and given the first link, say $s_{1}$, the distribution of the second link is $\kappa_{s_{1}}$ and so forth. This results in a potentially infinite sequence of random variables $\mathbf{X}=(\mathbf{X}_{0},\mathbf{X}_{1},\mathbf{X}_{2},...)$, whereby $\mathbf{X}_{n}$ represents the $n$-th link in the chain. Consequently, we get that $\mathbb{P}\left(\mathbf{X}_{n}=s\mid\mathbf{X}_{n-1}=z\right)=\kappa_{zs}$ for $n>0$ and $\mathbb{P}\left(\mathbf{X}_{0}=s_{0}\right)=1$.

Later on, we will deal with the (unconditional) distributions of the $n$-th links. To this end, w.l.o.g. we assume that $\mathcal{S}=\{1,...,k\}$. Therewith, we describe $\kappa$ as a matrix in $[0,1]^{k\times k}$ and $\pi$ as a column vector in $[0,1]^{k}$. Any quadratic matrix with non-negative entries and rows that sum to one is called stochastic matrix. Therewith, $\mathbb{P}\left(\mathbf{X}_{n}=s\right)=(\kappa^{n})_{s_{0}s}$, i.e. the $n$-fold matrix product of $\kappa$ evaluated at $s_{0}$ and $s$. A further generalization is to set the law of $\mathbf{X}_{0}$ to an arbitrary distribution $\pi$. This yields $\mathbb{P}\left(\mathbf{X}_{n}=s\right)=(\pi^{t}\kappa^{n})_{s}$, whereby the $t$ indicates transposition.

We say that a distribution $\pi$ is an invariant distribution of $\kappa$ if $\pi^{t}\kappa=\pi^{t}$. Thus, transitioning according to $\kappa$ doesn’t affect $\pi$. Once a link in a Markov chain with transition kernel $\kappa$ follows the law $\pi$, all subsequent links will do likewise. In this case, the chain is considered to be in equilibrium. Equilibrium can be enforced by setting the distribution of $\mathbf{X}_{0}$ to $\pi$, but it may also be reached (approximately) in the long run through convergence.

The foundation of MCMC sampling is that under some circumstances the distributions of the $n$-th links of a Markov chain converge towards an invariant distribution regardless of the initial state. Thus, by simulating such a chain until equilibrium is reached sufficiently, we may obtain an approximate sample of this invariant distribution. On these grounds, MCMC methods provide schemes to build Markov chains with a predefined unique invariant distribution.

The Gibbs sampler (Geman and Geman, 1984) builds a Markov chain with invariant distribution $\pi$ by decomposing $\pi$ into simpler conditional versions. This facilitates sampling of complex joint distributions, but is somewhat restricted in its ability to explore the state space.

The well-known Metropolis-Hastings algorithm (Metropolis et al., 1953; Hastings, 1970) is capable of incorporating user defined proposal distributions, which enables the exploration of the state space in any desired fashion. It further facilitates the application of complex statistical models to observed data.

1.1 Convergence towards and existence of invariant distributions

At first, we consider the convergence of the distributions of the $n$-th links of certain Markov chains. This convergence forms the very basis of MCMC sampling. At second, we provide a version of the Perron-Frobenius theorem, which gives further insights into the existence of invariant distributions. The following two fundamental properties of $\kappa$ are imposed by these theorems.

$\kappa$ is called irreducible if for every $z,s\in\mathcal{S}$ there exists an $n>0$ such that $(\kappa^{n})_{zs}>0$. Thus, regardless of the state the Markov chain starts in, every state can eventually be reached with positive probability.

$\kappa$ is called aperiodic if there exists an $N$ such that $(\kappa^{n})_{zs}>0$ for all $n>N$ and all $z,s\in\mathcal{S}$. Thus, regardless of the state the Markov chain starts in, in the long run, it can always reach any state immediately with positive probability.

The name irreducible suggests that the Markov chain does not divide $\mathcal{S}$ into separate, mutually inaccessible classes. In turn, aperiodicity excludes the case where parts of $\mathcal{S}$ can only be reached in a periodic fashion, for example only through either an even or odd number of transitions. Aperiodicity implies irreducibility, but not the other way around.

It is obvious that a periodic behavior may impede convergence of the distributions of the $n$-th links. For a general convergence theorem, aperiodicity is thus a necessary condition. We will now see that it is also sufficient. The following theorem is a simplified version of a convergence theorem for Markov chains over countable state spaces provided inter alia by Koenig (2005).

Theorem 1**.**

For an aperiodic stochastic matrix $\kappa$ with invariant distribution $\pi$, we get that $\lim\limits_{n\rightarrow\infty}(\kappa^{n})_{zs}=\pi_{s}$ for all $z,s\in\mathcal{S}$.

Proof.

Assume that $\mathbf{X}=(\mathbf{X}_{0},\mathbf{X}_{1},...)$ is a Markov chain with transition kernel $\kappa$ that starts with an arbitrary but fixed state $s_{0}$. Furthermore, consider the Markov chain $\mathbf{Z}=(\mathbf{Z}_{0},\mathbf{Z}_{1},...)$ with transition kernel $\kappa$ and initial distribution $\pi$, i.e. $\mathbb{P}\left(\mathbf{Z}_{0}=s\right)=\pi_{s}$ for all $s\in\mathcal{S}$. $\mathbf{X}$ and $\mathbf{Z}$ are supposed to be independent of each other.

Let $\mathbf{T}$ be the random variable that represents the first $n$ where $\mathbf{X}$ and $\mathbf{Z}$ equal, i.e. $\mathbf{T}=\min\{n\mid\mathbf{X}_{n}=\mathbf{Z}_{n}\}$. We want to show that $\mathbf{T}$ is finite with probability one. Due to the aperiodicity of $\kappa$, we may choose an $N$ such that $\kappa^{N}$ has solely positive entries. Let $\epsilon>0$ be the smallest entry of $\kappa^{N}$ and consider that

[TABLE]

Define further \mathbf{Y}_{n}=\raisebox{-0.85358pt}{\scalebox{1.2}{\mathds{1}}}\{n\leq\mathbf{T}\}\mathbf{X}_{n}+\raisebox{-0.85358pt}{\scalebox{1.2}{\mathds{1}}}\{n>\mathbf{T}\}\mathbf{Z}_{n}. The Markov chain $\mathbf{Y}=(\mathbf{Y}_{0},\mathbf{Y}_{1},...)$ starts with copying $\mathbf{X}$ and switches to $\mathbf{Z}$ as soon as both equal the first time. We are interested in the distribution of $\mathbf{Y}$. To this end, consider an arbitrary path $s_{1:n}\in\mathcal{S}^{n}$, define $\mathantt{p}_{s_{j:i}}=\prod_{\ell=j}^{i}\kappa_{s_{\ell-1}s_{\ell}}$ for $j,i=1,...,n$ and observe that

[TABLE]

whereby we used that $\mathbb{P}\left(\mathbf{Z}_{0:n}\neq s_{0:n}\right)+\sum_{\ell=0}^{n}\mathbb{P}\left(\mathbf{Z}_{\ell}=s_{\ell},\mathbf{Z}_{0:\ell-1}\neq s_{0:\ell-1}\right)=1$. This shows that $\mathbf{Y}$ is a Markov chain with transition kernel $\kappa$ and initial state $s_{0}$.

Altogether, we may state that

[TABLE]

and $|(\kappa^{n})_{s_{0}s}-\pi_{s}|=|\mathbb{P}\left(\mathbf{Y}_{n}=s,n\leq\mathbf{T}\right)-\mathbb{P}\left(\mathbf{Z}_{n}=s,n\leq\mathbf{T}\right)|\ \overset{n\rightarrow\infty}{\longrightarrow}\ 0$ for all $s_{0}\in\mathcal{S}$. ∎

Given a distribution $\pi$, MCMC methods seek an aperiodic transition kernel $\kappa$ with invariant distribution $\pi$. Thus, it is possible to sample approximately from $\pi$ by simulating a Markov chain with transition kernel $\kappa$ until equilibrium is reached to a sufficient extend. The last sample within this chain is then taken as a single approximate sample from $\pi$. In particular, this procedure is independent of the state the Markov chain has started in. The pace with which equilibrium is approached is referred to as the mixing time.

Now, we consider a version of the well-known Perron-Frobenius theorem (Frobenius, 1912). It is usually stated in a more general context and corresponding proofs can be fairly complicated. In turn, we provide our own convenient proof based on simple arithmetics and matrix algebra.

Theorem 2** (Perron-Frobenius Theorem).**

An irreducible transition kernel $\kappa$ has a unique invariant distribution $\pi$.

Proof.

Since any stochastic matrix has a right eigenvector with corresponding eigenvalue 1, it also has such a left eigenvector. In particular, any such left 1-eigenvector exhibits non-zero elements. Let $x$ be a left 1-eigenvector of $\kappa$. If $x$ has only non-negative or non-positive entries, we can immediately derive an invariant distribution $\pi$ of $\kappa$ through normalizing $x$, i.e. \pi=x\big{/}\sum_{s\in\mathcal{S}}x_{s}.

Assume now that $x$ exhibits positive entries for $s\in\bar{N}$ and negative entries for $s\in N$. The following applies

[TABLE]

Hence, the l.h.s and r.h.s. of (1) have to be zero, which implies that $\kappa_{zs}=0$ for all $z\in\bar{N}$ and $s\not\in\bar{N}$. Consequently, $(\kappa^{n})_{zs}=0$ for all $z\in\bar{N}$, $s\not\in\bar{N}$ and $n\in\mathbb{N}$.

Since the existence of positive and negative entries implies reducibility, we conclude that irreducibility implies that any left 1-eigenvector has either solely non-positive or non-negative entries. Thus, an irreducible transition kernel $\kappa$ exhibits at least one invariant distribution $\pi$.

Finally, assume that there is a second invariant distribution $\pi^{\prime}\neq\pi$. In order to represent a distribution, not all components of $\pi$ can be either larger or smaller than the components of $\pi^{\prime}$. Thus, $\pi-\pi^{\prime}$ must have positive and negative entries. However, $\pi-\pi^{\prime}$ is a left 1-eigenvector of $\kappa$ and thus, $\kappa$ can’t be irreducible, which contradicts the existence of $\pi^{\prime}$. ∎

The Perron-Frobenius theorem shows that invariant distributions can certainly be found for an abundance of stochastic matrices, especially for aperiodic ones. In the context of MCMC, it is, however, only a nice-to-have result and not utterly necessary. In fact, there is great freedom in choosing aperiodic transition kernels that exhibit a prespecified invariant distribution and each MCMC method provides its very own approach to do so.

1.2 The Gibbs Sampler

The Gibbs sampler (Geman and Geman, 1984) is a primal MCMC sampling algorithm which is based on a decomposition of the objective distribution $\pi$ into conditional versions. It is mainly used to sample from the joint distribution of a set of random variables. Thereby each step may involve sampling from a single random variable given the remaining random variables conditioned on the last sample.

Assume that $\mathcal{S}=\mathcal{V}^{m}$, with a finite space $\mathcal{V}$ and let $\pi$ be a distribution over $\mathcal{S}$. We proceed in a step wise manner. In step $i$, given the sample $z\in\mathcal{S}$ from the last step $i-1$, we choose $j=(i\mod m)+1$ and sample from the transition kernel $\kappa^{j}_{z}$ which is defined through

[TABLE]

$\kappa^{j}_{zs}$ is the conditional version of $\pi_{s}$ conditioned on $s_{\ell}=z_{\ell}$ for all $\ell$ except the $j$-th one.

We get that $\pi$ is an invariant distribution of $\kappa^{j}$ since $\sum_{z\in\mathcal{S}}\pi_{z}\kappa^{j}_{zs}=\sum_{z\in\mathcal{S}}\pi_{s}\kappa^{j}_{sz}=\pi_{s}$. However, a single Gibbs step is generally not aperiodic since it manipulates only one single coordinate. This may be fixed by considering $m$ subsequent Gibbs steps as one transition within the Markov chain of the ensuing MCMC algorithm.

Unfortunately, the possible transitions are completely determined by $\pi$. For example, assume that $\mathcal{S}=\{0,1\}^{2}$ and $\pi_{(0,0)}=\pi_{(1,1)}=0$. In this case, we can never get from $(1,0)$ to $(0,1)$ and thus, the chain is not irreducible. This poses one of several reasons why it can be important to be able to traverse the state space in a more arbitrary fashion.

A famous application of the Gibbs sampler is the Ising model (Ising, 1925). There, $\mathcal{S}$ consists of the positive or negative states of the grid points of a finite grid, whereby independence is induced by spatial separation. This yields very simple sampling steps, each conducted on a single grid point given all the other, but essentially only its neighboring grid points. Higdon (1998) provides very vivid and more sophisticated treatments of the Ising model through Gibbs sampling.

1.3 The Metropolis-Hastings algorithm

The well-known Metropolis-Hastings algorithm is an MCMC sampler that traverses through the state space by means of a user defined proposal distribution. Inherent for this sampler is that each proposed state undergoes an accept-reject step which decides whether the proposed state or the previous sample is chosen to be the next sample. This step alone secures the required invariance and thus, gives the user great freedom in exploring the state space. A primal version was first published in Metropolis et al. (1953) and then extended in Hastings (1970).

The user provides a transition kernel $\rho$ over $\mathcal{S}$ which is referred to as the proposal. In step $i$, given the previous sample $z$, a new state $s$ is proposed according to $\rho_{z}$ and accepted with probability

[TABLE]

whereby we agree that dividing by zero yields $\infty$. The new sample is than either $s$ if accepted or $z$ otherwise.

The invariance can be shown by means of the following condition. A transition kernel $\kappa$ preserves the detailed balance w.r.t. $\pi$ if $\pi_{s}\kappa_{sz}=\pi_{z}\kappa_{zs}$ for all $s,z\in\mathcal{S}$. Detailed balance implies invariance since $\sum_{z\in\mathcal{S}}\pi_{z}\kappa_{zs}=\pi_{s}\sum_{z\in\mathcal{S}}\kappa_{sz}=\pi_{s}$, but the opposite implication does not hold in general.

A Markov chain that preserves the detailed balance w.r.t. another distribution $\pi$ is said to be time-reversible in equilibrium or simply called reversible. That is, if $\mathbf{X}_{n}$ follows the law $\pi$, we get $\mathbb{P}\left(\mathbf{X}_{n+1}=s\mid\mathbf{X}_{n}=z\right)=\mathbb{P}\left(\mathbf{X}_{n}=z\mid\mathbf{X}_{n+1}=s\right)$.

The Metropolis-Hastings algorithm describes a reversible Markov chain. This is trivial for transitions from $z$ to $s$ with $s=z$. For $s\neq z$, consider that

[TABLE]

The aperiodicity of the resulting Markov chain has to be met jointly by the acceptance probability and the proposal. Equation (3) shows that a non-trivial transition from $z$ to $s$ can only take place if the corresponding backward transition is also feasible, i.e. if $\pi_{s}\rho_{sz}>0$. In the most extreme case, this could mean that we apply an irreducible and aperiodic proposal, but the resulting Metropolis-Hastings Markov chain will never move away from the initial state and is thus not irreducible.

Apart from the ability to traverse the state space in a user defined fashion, a very important advantage of this sampling algorithm is that normalization constants w.r.t. $\pi$ cancel out in (3). Thus, we may sample from conditional versions of $\pi$ by using $\pi$ directly in the acceptance probability. This enables convenient data-driven inference approaches by means of carefully designed and complex statistical models, without the need to compute cumbersome marginalizations over unobserved states.