马尔可夫链蒙特卡罗MCMC

马尔可夫链蒙特卡罗

马尔可夫链蒙特卡罗(Markov Chain Monte Carlo, MCMC)是一类利用马尔可夫链从特定分布中随机采样的算法

具体来说，MCMC定义了一个满足遍历定理的马尔可夫链，使其平稳分布就是目标分布

这样的马链在时间趋于无穷时，样本分布趋于平稳分布，即假设已长时间运行到时刻$m$，则之后的时间里得到的样本集合${X{m+1},X{m+2},…X_{m+n}}$就是目标概率分布的抽样结果

那么MCMC的关键问题就是如何定义满足条件的马尔可夫链

常用的MCMC方法有Metropolis-Hastings算法(M-H算法)和Gibbs抽样

Metropolis-Hastings

首先令$Q$是一个不可约马尔可夫链，满足任意$q(i,j)>0$，称为建议分布(proposal distribution)

再定义$\alpha(i,j)$，称为接受分布(acceptance distribution)，则M-H算法构造的马尔可夫链$P$转移概率为

$$\begin{aligned} &P_{ij}=q(i,j)\alpha(i,j)&&,i\neq j\\ &P_{ii}=q(i,i)+\sum_{k\neq i}q(i,k)(1-\alpha(i,k)) \end{aligned}$$

其直观解释是，假设$Xn=i$，生成一个随机变量$Y$满足$P(Y=j)=q(i,j)$，若$Y=j$，则$X{n+1}$以概率$\alpha(i,j)$等于$j$，以概率$1-\alpha(i,j)$等于$i$

令$\pi$是马尔可夫链$P$的平稳分布，则接受分布$\alpha(i,j)$表达式为

$$\alpha(i,j)=\min\{1, \frac{\pi_jq(j,i)}{\pi_iq(i,j)}\}$$

则转移概率$P_{ij}$化为

$$P_{ij}=\begin{cases} q(i,j) && ,\pi_jq(j,i)\geq \pi_iq(i,j)\\ q(j,i)\frac{\pi_j}{\pi_i} && ,\pi_jq(j,i)< \pi_iq(i,j) \end{cases}$$

容易验证这个马尔可夫链是时间可逆的，因为

$$\begin{aligned} \pi_iP_{ij} &=\pi_iq(i,j)\min\{1, \frac{\pi_jq(j,i)}{\pi_iq(i,j)}\} \\ &=\min\{\pi_iq(i,j), \pi_jq(j,i)\}\\ &=\pi_jq(j,i)\min\{\frac{\pi_iq(i,j)}{\pi_jq(j,i)}, 1\} \\ &=\pi_jq(j,i) \end{aligned}$$

通俗的解释，若$\alpha(i,j)=\frac{\pijq(j,i)}{\pi_iq(i,j)}$，则必有$\alpha(j,i)=1$
因此$\pi_iP{ij}=\piiq(i,j)\frac{\pi_jq(j,i)}{\pi_iq(i,j)}=\pi_jq(j,i)$和$\pi_jP{ji}=\pi_jq(j,i)$，反之亦然

对于可逆马尔可夫链，从运行时间足够长的某时刻开始，可通过随机游走获得平稳分布的样本

综上所述，令马尔可夫链$P$的平稳分布$\pi$等于目标分布$p(x)$，则有如下M-H算法

1. 选择任意初值$x_0$

2. 对$t=1,2,…,n$循环执行

   (1) 设$x_{t=1}=i$，则按照建议分布$q(i,j)$随机抽取一个状态$j$

   (2) 计算$\alpha(i,j)=\min{1, \frac{p(j)q(j,i)}{p(i)q(i,j)}}$

   (3) 按照$u\sim U(0,1)$ 随机抽取一个数$u$，若$u\leq \alpha(i,j)$，则$x_t=j$，否则$x_t=i$

3. 运行时间足够长后，通过随机游走得到样本集合${X{m+1},X{m+2},…X_{m+n}}$

其中建议分布$q(i,j)$的一般选择标准是采样效率高、且形状与目标分布尽可能相似

此外为计算接受分布$\alpha(i,j)$，建议分布必须可逆，常用的建议分布有高斯分布、均匀分布等等

Gibbs抽样

Gibbs抽样用于多元变量联合分布的抽样，假设目标多元变量的概率分布为$p(x)=P(X_1,X_2,…,X_n)$，算法首先指定初始样本$X^{(0)}=(x_1^{(0)}, x_2^{(0)},…, x_n^{(0)})$

对于第$i$轮迭代，样本$X^{(i)}=(x_1^{(i)}, x_2^{(i)},…, x_n^{(i)})$按照如下条件概率依次抽取$x_1^{(i)}, x_2^{(i)},…, x_n^{(i)}$获得

$$\begin{aligned} &p(x_1^{(i)}|x_2^{(i-1)},...,x_n^{(i-1)})\\ &p(x_j^{(i)}|x_1^{(i)},...,x_{j-1}^{(i)}, x_{j+1}^{(i-1)},...,x_n^{(i-1)}), \quad j=2,...,n-1 \\ &p(x_n^{(i)}|x_1^{(i)},...,x_{n-1}^{(i)}) \end{aligned}$$

例如使用Gibbs抽样从二元正态分布中随机抽样

$$(X,Y)\sim N(0,\Sigma),\quad \Sigma=\begin{bmatrix}1 & \rho \\ \rho & 1\end{bmatrix}$$

条件概率为一元正态分布

$$P(X|Y=y)=N(\rho y,1-\rho^2) \\ P(Y|X=x)=N(\rho x,1-\rho^2) \\$$

假设初始样本为$X^{(0)}=(x_1^{(0)},x_2^{(0)})$

第一轮，按照$x_1^{(1)}\sim N(\rho x_2^{(0)},1-\rho^2)$抽样得$x_1^{(1)}$，按照$x_2^{(1)}\sim N(\rho x_1^{(1)},1-\rho^2)$抽样得$x_2^{(1)}$，获得样本$(x_1^{(1)},x_2^{(1)})$

依此类推，第$i$轮，按照$x_1^{(i)}\sim N(\rho x_2^{(i-1)},1-\rho^2)$抽样得$x_2^{(i)}$，按照$x_2^{(i)}\sim N(\rho x_1^{(i)},1-\rho^2)$抽样得$x_2^{(i)}$，获得样本$(x_1^{(i)},x_2^{(i)})$

最终得到样本集${x^{(m)},x^{(m+1)},…,x^{(m+n)}}$

实际上Gibbs抽样可以看作特殊的M-H算法，即建议分布为条件分布，而接受分布始终为1