3. Markov Chain & Markov Sampling

#math/chapman-Kolmogorov-equation#math/gibbs-sampling#math/markov-chain#math/mcmc-sampling#math/metropolis-hastings-algorithm

1. Markov chains

(1) State transition probability

Markov property The next state depends only on the current state, not on the past states.

Let $(X_t{)}_{t\ge 0}$ be a discrete-time Markov chain on a state space $\mathcal{S}=\{a_1,a_2,\dots \}$. It satisfies the Markov property:

$$\begin{array}{}\text{(1.1.1)}& Pr(X_{t+1}=j\mid X_t=i,X_{t-1},\dots ,X_{0)}=Pr(X_{t+1}=j\mid X_t=i).\end{array}$$

Thus, conditional on the current state, the next state is independent of the earlier history. For a time-homogeneous chain, this conditional probability does not depend on $t$. Define the one-step transition probabilities and transition matrix by

$$\begin{array}{}\text{(1.1.2)}& P_{ij}=Pr(X_{t+1}=a_j\mid X_t=a_i),P=(P_{ij}).\end{array}$$

This note uses the row-vector convention: row $i$ gives the probabilities of moving from $a_i$ to the possible next states. Therefore,

$$\begin{array}{}\text{(1.1.3)}& P_{ij}\ge 0,\sum _j{P}_{ij}=1.\end{array}$$

For example, a finite-state transition diagram may be represented by the matrix $P$, whose entry $P_{ij}$ is the label on the directed edge $a_i\to a_j$.

(2) Multi-step transition probabilities

We define ==the $n$-step transition probability== as :

$$\begin{array}{}\text{(1.2.1)}& p_{ij}^{(n)}=Pr(X_{t+n}=a_j\mid X_t=a_i).\end{array}$$

Let $P^{(n)}=(p_{ij}^{(n)})$. The Chapman-Kolmogorov equation (C-K equation) is defined as the transition probability relation on state transition i-> k -> j summed on every possible intermediate states $k$ :

$$\begin{array}{}\text{(1.2.2)}& \boxed{{\displaystyle p_{ij}^{(m+n)}=\sum _{k\in \mathcal{S}}{p}_{ik}^{(m)}{p}_{kj}^{(n)}.}}\end{array}$$

$p_{ij}^{(m+n)}=Pr(X_{t+m+n}=a_j\mid X_t=a_i)$ is the probability that the chain reaches state $a_j$ after $m+n$ steps, given that it starts at $a_i$.

In particular, for a time-homogeneous chain, it satisfies the Markov condition, the relation of :

$$\begin{array}{}\text{(1.2.3)}& P^{(n)}=P^n,p_{ij}^{(n)}=\sum _{k\in \mathcal{S}}{P}_{ik}{p}_{kj}^{(n-1)}.\end{array}$$

Equivalently, expanding over all intermediate states gives :

$$\begin{array}{}\text{(1.2.4)}& p_{ij}^{(n)}=\sum _{k_1,\dots ,k_{n-1}\in \mathcal{S}}{P}_{i{k}_1}{P}_{k_1{k}_2}\cdots P_{k_{n-1}j}.\end{array}$$

If the initial distribution is the row vector ${\mu }_0$, then the distribution after $n$ steps is

$$\begin{array}{}\text{(1.2.5)}& {\mu }_n={\mu }_0{P}^n.\end{array}$$

(3) Limiting distribution

A distribution $\pi$ is a limiting distribution of state transformmatrix P when

$$\begin{array}{}\text{(1.3.1)}& \underset{n\to \mathrm{\infty }}{lim}{\mu }_0{P}^n=\pi \end{array}$$

for every initial distribution ${\mu }_0$. For a finite Markov chain, irreducibility and aperiodicity are sufficient for a unique limiting distribution. More generally, suitable recurrence conditions are also required on an infinite state space.

2. Stationary distributions and detailed balance

(1) Stationary distribution

A probability distribution $\pi =({\pi }_1,{\pi }_2,\dots )$ is stationary for $P$ if :

$$\begin{array}{}\text{(2.1.1)}& \boxed{{\displaystyle \pi P=\pi ,\sum _i{\pi }_{i=1,}{\pi }_i\ge 0.}}\end{array}$$

Component wise, the stationarity condition is :

$$\begin{array}{}\text{(2.1.2)}& {\pi }_j=\sum _i{\pi }_i{P}_{ij}.\end{array}$$

If a chain is initialized from $\pi$, it remains distributed as $\pi$ after every transition. In MCMC, $\pi$ is the target distribution; after convergence, the chain states can be used as dependent samples from $\pi$.

(2) Detailed balance

To make sure a distribution is a stationary distribution, A convenient sufficient condition for stationarity is detailed balance (or reversibility) :

$$\begin{array}{}\text{(2.2.1)}& {\pi }_i{P}_{ij}={\pi }_j{P}_{ji}\text{for all }i,j.\end{array}$$

Indeed, summing over $i$ and using $\sum _i{P}_{ji}=1$ gives

$$\begin{array}{}\text{(2.2.2)}& \sum _i{\pi }_i{P}_{ij}=\sum _i{\pi }_j{P}_{ji}={\pi }_j.\end{array}$$

Detailed balance is sufficient, but not necessary, for $\pi$ to be stationary.

3. Markov Chain Monte Carlo (MCMC)

MCMC constructs a Markov chain whose stationary distribution is a chosen target density or mass function $\pi (x)$. It is especially useful when direct independent sampling from $\pi$ is difficult, while evaluating $\pi (x)$ up to a normalizing constant is feasible.

Starting from an arbitrary initial state $X_0=x_0$, run a transition kernel that leaves $\pi$ invariant. After a suitable burn-in period, retain states $X_t$ as samples. Because consecutive states are correlated, effective sample size is generally smaller than the number of retained iterations.

Target density need not be normalized In Metropolis-Hastings, only ratios such as $\pi (y)/\pi (x)$ are used. Hence one may use an unnormalized target $\tilde{\pi }(x)=c\pi (x)$ when the normalizing constant $c$ is unknown.

4. Metropolis-Hastings sampling

(1) Introduction

Let $q(y\mid x)$ be a proposal distribution, if we want to sample from this distribution, we can use the

given the current state $x$, it generates a candidate $y$. The Metropolis-Hastings (MH) algorithm is to accept/reject a sample by the probability :

1. Set $X_0=x_0$.
2. Given $X_t=x$, propose $Y\sim q(\cdot \mid x)$.
3. Compute the acceptance probability :

$$\begin{array}{}\text{(4.1.1)}& \alpha (x,y)=min\{1,\frac{\pi (y)q(x\mid y)}{\pi (x)q(y\mid x)}\}.\end{array}$$

Note here $q$ is the state transforming matrix 4. Draw $U\sim \mathrm{Uniform}(0,1)$. if the probability satisfies, accept $Y$ as the next state.

$$\begin{array}{}\text{(4.1.2)}& X_{t+1}=\{\begin{array}{ll}Y,& U\le \alpha (x,Y),\\ x,& U>\alpha (x,Y).\end{array}\end{array}$$

The resulting transition kernel includes both accepted moves and the probability of staying at the current state. It satisfies detailed balance with respect to $\pi$, so $\pi$ is stationary.

For a symmetric proposal, $q(y\mid x)=q(x\mid y)$, the acceptance probability simplifies to:

$$\begin{array}{}\text{(4.1.3)}& \alpha (x,y)=min\{1,\frac{\pi (y)}{\pi (x)}\}.\end{array}$$

Thus, moves toward higher-density regions are always accepted, while moves toward lower-density regions may still be accepted. The latter feature prevents the chain from becoming trapped near a single mode.

Proposal tuning A proposal with very small steps tends to have high acceptance but strong autocorrelation. A proposal with very large steps tends to have low acceptance. Good sampling requires a balance between movement and acceptance.

(2) An Example : Drug Effect problem

We know the drug design parameter $x$ from the sample, the outcome from patients $y$ have following distribution, $y=1$ means the drug takes effect :

$$\begin{array}{}\text{(4.2.1)}& p(\theta )=p(y_i=1|\theta )=\frac{1}{1+e^{-({\beta }_1+{\beta }_{2}x)}}\theta =({\beta }_1,{\beta }_2)D=y\end{array}$$

Here ${\beta }_1$ and ${\beta }_2$ are unknown drug effect parameters, we know $(x_i,y_i)$, then the likelihood of one observation is :

$$\begin{array}{}\text{(4.2.2)}& \boxed{{\displaystyle P(y_i|\theta ,x_i)=p_i(\theta {)}^{y_i}[1-p_i(\theta ){]}^{1-y_i}}}\end{array}$$

We want to get the posterior probability :

$$\begin{array}{}\text{(4.2.3)}& P(\theta |x,y)=P({\beta }_1,{\beta }_2|x_i,y_i)\end{array}$$

using the Bayes formula, we have :

$$\begin{array}{}\text{(4.2.4)}& P({\beta }_1,{\beta }_2|y,x)=\frac{P(y|{\beta }_1,{\beta }_2,x)P({\beta }_1)P({\beta }_2)}{P(y)}=\frac{[\prod _{i=1}^{n}P(y_i|\theta ,x_i)]\times P({\beta }_1)P({\beta }_2)}{\iint \dots d{\beta }_{1}d{\beta }_2}\end{array}$$

where the first term in nominator comes from (4.2.2). The denominator holds because the ${\beta }_1$ and ${\beta }_2$ covers the space. So, in the M-H sampling, we ==don't compute the denominator since it's intractable==. Instead, we set the proposal distribution that prop to the posterior :

$$\begin{array}{}\text{(4.2.5)}& \tilde{\pi }(\theta )=P(y|x,\theta )P(\theta )=\prod _{i=1}^n[p_i(\theta {)}^{y_i}[1-p_i(\theta ){]}^{1-y_i}]\times P({\beta }_1)P({\beta }_2)\propto P({\beta }_1,{\beta }_2|y)\end{array}$$

where $\tilde{\pi }$ is called unnormalized posterior.

Then use following acceptance computation to sample from posterior distribution :

$$\begin{array}{}\text{(4.2.6)}& \alpha (y,x)=min\{1,\frac{\tilde{\pi }({\theta }_{new})}{\tilde{\pi }({\theta }_{old})}\}\end{array}$$

5. Gibbs sampling

Gibbs sampling is an MCMC method for a multivariate target distribution $\pi (x_1,\dots ,x_d)$ when the full conditional distributions can be sampled.

For a two-dimensional target $\pi (x_1,x_2)$, a systematic Gibbs iteration is sample each according to its own distribution $\pi$:

$$\begin{array}{}\text{(5.1)}& \begin{array}{rl}{X}_1^{(t+1)}& \sim \pi (x_1\mid X_2^{(t)}),\\ X_2^{(t+1)}& \sim \pi (x_2\mid X_1^{(t+1)}).\end{array}\end{array}$$

More generally, update one coordinate at a time from its full conditional distribution:

$$\begin{array}{}\text{(5.2)}& X_k^{(t+1)}\sim \pi (x_k\mid X_1^{(t+1)},\dots ,X_{k-1}^{(t+1)},X_{k+1}^{(t)},\dots ,X_d^{(t)}).\end{array}$$

For example, when $x_1$ is held fixed, moving from $(x_1,x_2)$ to $(x_1,x_2^{\prime })$ has transition density $\pi (x_2^{\prime }\mid x_1)$. The corresponding coordinate-update kernel satisfies

$$\begin{array}{}\text{(5.3)}& \pi (x_1,x_2)\pi (x_2^{\prime }\mid x_1)=\pi (x_1,x_2^{\prime })\pi (x_2\mid x_1),\end{array}$$

because $\pi (x_1,x_2)=\pi (x_1)\pi (x_2\mid x_1)$. Consequently, each coordinate-update kernel leaves $\pi$ invariant, and no accept/reject step is required: every Gibbs proposal is accepted.

Systematic versus random scan Each individual coordinate update is reversible. A systematic sequence of coordinate updates still leaves $\pi$ invariant, although the combined full-sweep kernel is not necessarily reversible. Randomly selecting a coordinate at each iteration yields a reversible mixture kernel.

Gibbs sampling can mix slowly when components are strongly correlated, because each move changes only one coordinate. Blocking correlated variables or using more advanced samplers can improve mixing.

References

1. C. P. Robert and G. Casella, Monte Carlo Statistical Methods, 2nd ed., Springer, 2004.
2. W. R. Gilks, S. Richardson, and D. J. Spiegelhalter (eds.), Markov Chain Monte Carlo in Practice, Chapman & Hall, 1996.
3. Metropolis-Hastings algorithm (Wikipedia)
4. Gibbs sampling (Wikipedia)