马尔科夫链蒙特卡罗法（Markov Chain Monte Carlo，MCMC）

$$\begin{align} \end{align}$$

1 蒙特卡洛法

蒙特卡洛法要解决的问题是，假设概率分布的定义已知，通过抽样获得概率分布的随机样本，并通过得到的随机样本对概率分布的特征进行分析。

蒙特卡洛的法有直接抽样法、接受-拒绝抽样法、重要性抽样法等，其核心是随机抽样。

接受-拒绝算法：
输入：抽样的目标概率分布的概率密度函数$p\left(x\right)$；
输出：概率分布的随机样本$x_1,x_2,\cdots,x_n$。

1. 选择概率密度函数为$q\left(x\right)$的概率分布作为建议分布，使其对任一$x$满足$cq\left(x\right)\geq \left(x\right)$，其中$c>0$。
2. 按照建议分布$q\left(x\right)$随机抽样得到样本$x’$，再按照均匀分布在$\left(0,1\right)$范围内抽样得到$u$。
3. 如果$u\leq\frac{p\left(x’\right)}{cq\left(x’\right)}$，则将$x’$作为抽样结果；否则，回到步骤2.。
4. 直至得到$n$个随机样本，结束。

2 马尔科夫链

假设在时刻$0$的随机变量$X_0$遵循概率分布$P\left(X_0\right)=\pi_0$，称为初始状态分布。在某个时刻$t\geq 1$的随机变量$X_t$与前一个时刻的随机变量$X_{t-1}$之间有条件分布$P\left(X_t|X_{t-1}\right)$，如果$X_t$只依赖于$X_{t-1}$，这一性质称为马尔科夫性，即

$$P\left(X_t|X_0,X_1,\cdots,X_{t-1}\right)=P\left(X_t|X_{t-1}\right), t=1,2,\cdots$$

具有马尔科夫性的随机序列$X=\{X_0,X_1,\cdots,X_t,\cdots\}$称为马尔科夫链或马尔科夫过程。每个随机变量$X_t\left(t=0,1,2,\cdots\right)$的取值集合相同，称为状态空间$\mathcal{S}$。

离散状态马尔科夫链$X=\{X_0,X_1,\cdots,X_t,\cdots\}$，随机变量$X_t\left(t=0,1,2,\cdots\right)$定义在离散空间$\mathcal{S}$，转移概率分布可以由矩阵表示。

若马尔科夫链咋时刻$t-1$处于状态$j$，在时刻$t$移动到状态$i$，将转移概率记作

$$p_{ij}=\left(X_t=i|X_{t-1}=j\right),i=1,2,\cdots;j=1,2,\cdots$$

满足

$$p_{ij}\geq 0, \sum_i p_{ij}=1$$

马尔科夫链的转移矩阵可表示为

$$P=\left[p_{ij}\right]$$

马尔科夫链$X=\{X_0,X_1,\cdots,X_t,\cdots\}$在时刻$t,t=0,1,2,\cdots$的概率分布，称为时刻$t$的状态分布，记作

$$\pi\left(t\right)=\left[\begin{matrix} \pi_1\left(t\right) \\ \pi_2\left(t\right) \\ \vdots \\ \end{matrix}\right]$$

其中$\pi_i\left(t\right)$表示时刻$t$状态为$i$的概率$P\left(X_t=i\right)$

$$\pi_i\left(t\right)=P\left(X_t=i\right),i=1,2,\cdots$$

特别地，马尔科夫链的初始状态分布可表示为

$$\pi\left(0\right)=\left[\begin{matrix} \pi_1\left(0\right) \\ \pi_2\left(0\right) \\ \vdots \\ \end{matrix}\right]$$

其中$\pi_i\left(0\right)$表示时刻$0$状态为$i$的概率$P\left(X_t=i\right)$

$$\pi_i\left(0\right)=P\left(X_t=i\right)$$

通常初始分布$\pi\left(0\right)$的向量只有一个分量是$1$，其余分量是$0$，表示马尔科夫链从一个具体状态开始。

马尔科夫链$X$在时刻$t$的状态分布，可以由在时刻$t-1$的状态分布以及转移概率分布决定

$$\pi\left(t\right)=P\pi\left(t-1\right)$$

设有马尔科夫链$X=\{X_0,X_1,\cdots,X_t,\cdots\}$，其状态空间为$\mathcal{S}$，转移概率矩阵为$P=\left[p_{ij}\right]$，如果状态空间$\mathcal{S}$上的一个分布

$$\pi=\left[\begin{matrix} \pi_1 \\ \pi_2 \\ \vdots \\ \end{matrix}\right]$$

使得

$$\pi=P\pi$$

则称$\pi$为马尔科夫链$X=\{X_0,X_1,\cdots,X_t,\cdots\}$的平稳分布。

连续状态马尔科夫链$X=\{X_0,X_1,\cdots,X_t,\cdots\}$，随机变量$X_t\left(t=0,1,2,\cdots\right)$定义在连续状态空间$\mathcal{S}$，转移状态分布由概率转移核或转移核表示。

设$\mathcal{S}$是连续状态空间，对任意的$x\in\mathcal{S},A\subset\mathcal{S}$，转移核定义为

$$P\left(x,A\right)=\int_A p\left(x,y\right)\mathrm{d}y$$

其中，$p\left(x,\centerdot\right)$是概率密度函数，满足$p\left(x,\centerdot\right)\geq0,P\left(x,\mathcal{S}\right)=\int_{\mathcal{S}}p\left(x,y\right)\mathrm{d}y=1$。转移核$P\left(x,A\right)$表示从$x\thicksim A$的转移概率

$$P\left(X_t=A|X_{t-1}=x\right)=P\left(x,A\right)$$

有时也将概率密度函数$p\left(x,\centerdot\right)$称为转移核。

若马尔科夫链的转态空间$\mathcal{S}$上的概率分布$\pi\left(x\right)$满足条件

$$\pi\left(y\right)=\int p\left(x,y\right)\pi\left(x\right)\mathrm{d}x,\quad\forall y\in\mathcal{S}$$

则称分布$\pi\left(x\right)$为该马尔科夫链的平稳分布。等价地，

$$\pi\left(A\right)=\int P\left(x,A\right)\pi\left(x\right)\mathrm{d}x,\quad\forall A\subset\mathcal{S}$$

或简写成

$$\pi=P\pi$$

3 马尔科夫链蒙特卡洛法——Metropolis-Hastings算法

假设要抽样的概率分布为$p\left(x\right)$。MH算法采用转移核为$p\left(x,x^{‘}\right)$的马尔科夫链：

$$p\left(x,x^{'}\right)=q\left(x,x^{'}\right)\alpha\left(x,x^{'}\right)$$

其中，$q\left(x,x^{‘}\right)$为建议分布，$\alpha\left(x,x^{‘}\right)$为接受分布。

建议分布$q\left(x,x^{‘}\right)$是另一个马尔科夫链的转移核，并且其概率值恒不为$0$，同时是一个容易抽样的分布。

接受分布$\alpha\left(x,x^{‘}\right)$是

$$\alpha\left(x,x^{'}\right)=\min\left\{1,\frac{p\left(x^{'}\right)q\left(x^{'},x\right)}{p\left(x\right)q\left(x,x^{'}\right)}\right\}$$

转移核可表示为

$$p\left(x,x^{'}\right)= \left\{ \begin{array} & q\left(x,x^{'}\right) & p\left(x^{'}\right)q\left(x^{'},x\right)\geq p\left(x\right)q\left(x,x^{'}\right)\\ q\left(x,x^{'}\right) \frac{p\left(x^{'}\right)}{p\left(x\right)} & p\left(x^{'}\right)q\left(x^{'},x\right)< p\left(x\right)q\left(x,x^{'}\right) \end{array} \right.$$

转移核为$p\left(x,x^{‘}\right)$的马尔科夫链上的随机游走以以下方式进行：如果在时刻$t-1$处于状态$x$，即$X_{t-1}=x$，则先按建议分布$q\left(x,x^{‘}\right)$抽样一个候选状态$x^{‘}$，然后按照接受分布$\alpha\left(x,x^{‘}\right)$抽样决定是否接受状态$x^{‘}$。以概率$\alpha\left(x,x^{‘}\right)$接受$x^{‘}$，$X_t=x^{‘}$；以概率$1-\alpha\left(x,x^{‘}\right)$拒绝$x^{‘}$，$X_t=x$。

具体地，在区间$\left(0,1\right)$上的均匀分布中抽取一个随机数$u$，决定时刻$t$的状态

$$X_t= \left\{ \begin{array} & x^{'}, & u\leq\alpha\left(x,x^{'}\right) \\ x, & u>\alpha\left(x,x^{'}\right) \end{array} \right.$$

可以证明，转移核为$p\left(x,x^{‘}\right)$的马尔科夫链的平稳状态就是$p\left(x\right)$，即要抽样的目标分布。

建议分布的对称形式，即对任意的$x$和$x^{‘}$有

$$q\left(x,x^{'}\right)=q\left(x^{'},x\right)$$

这样的建议分布称为Metropolis选择。此时，接受分布$\alpha\left(x,x^{‘}\right)$简化为

$$\alpha\left(x,x^{'}\right)=\min\left\{1,\frac{p\left(x^{'}\right)}{p\left(x\right)}\right\}$$

建议分布的独立抽样形式，即$q\left(x,x^{‘}\right)$与当前状态$x$无关，$q\left(x,x^{'}\right)=q\left(x^{'}\right)$
此时，接受分布$\alpha\left(x,x^{‘}\right)$可写成

$$\alpha\left(x,x^{'}\right)=\min\left\{1,\frac{w\left(x^{'}\right)}{w\left(x\right)}\right\}$$

其中，$w\left(x^{‘}\right)=p\left(x^{‘}\right)/q\left(x^{‘}\right)$，$w\left(x\right)=p\left(x\right)/q\left(x\right)$。

Metropolis-Hastings算法：
输入：待抽样的目标分布的密度函数$p\left(x\right)$，收敛步骤$m$，迭代步骤$n$
输出：$p\left(x\right)$的随机样本$x_{m+1},x_{m+2},\cdots,x_n$

1. 任意选择一个初始状态值$X_0$
2. 对$i=1,2,\cdots,n$ 循环执行
   (a)设状态 $X_{i-1}=x$，按照建议分布$q\left(x,x^{‘}\right)$随机抽样一个候选状态$x^{‘}$
   (b)计算接受概率 $$\alpha\left(x,x^{'}\right)=\min\left\{1,\frac{p\left(x^{'}\right)q\left(x^{'},x\right)}{p\left(x\right)q\left(x,x^{'}\right)}\right\}$$ (c)从区间$\left(0,1\right)$中按均匀分布随机抽取数$u$
   若$u\leq\alpha\left(x,x^{‘}\right)$，则状态$X_i=x^{‘}$；否则，状态$X_i=x$。
3. 得到样本集合$\{x_{m+1},x_{m+2},\cdots,x_n\}$

4 单分量Metropolis-Hastings算法

假设马尔科夫链的状态由$k$维随机变量表示

$$x=\left(x_1,x_2,\cdots,x_k\right)^\top$$

其中，$x_j$表示随机变量$x$的第$j$个分量，$j=1,2,\cdots,k$

马尔科夫链在时刻$i$的状态

$$x^{\left(i\right)}=\left(x_1^{\left(i\right)},x_2^{\left(i\right)},\cdots,x_k^{\left(i\right)}\right)$$

其中，$x_j^{\left(i\right)}$是随机变量$x^{\left(i\right)}$的第$j$个分量，$j=1,2,\cdots,k$。

单分量Metropolis-Hastings算法由下面的$k$步迭代实现Metropolis-Hastings算法的一次迭代。

设在第$i-1$次迭代结束时分量$x_j$的取值为$x_j^{\left(i-1\right)}$，在第$i$次迭代的第$j$步，对分量$x_j$根据Metropolis-Hastings算法更新，得到其新的取值$x_j^{\left(i\right)}$。

首先，由建议分布$q\left(x_j^{\left(i-1\right)},x_j|x_{-j}^{\left(i\right)}\right)$抽样产生分量$x_j$的候选值$x_j^{‘\left(i\right)}$。$x_{-j}^{\left(i\right)}$表示在第$i$次迭代的第$j-1$步后的$x^{\left(i\right)}$除去$x_j^{\left(i-1\right)}$的所有值，即

$$x_{-j}^{\left(i\right)}=\left(x_1^{\left(i\right)},\cdots,x_{j-1}^{\left(i\right)},x_{j+1}^{\left(i-1\right)},\cdots,x_{k}^{\left(i-1\right)}\right)^\top$$

其中分量$x_1^{\left(i\right)},\cdots,x_{j-1}^{\left(i\right)}$已经更新。

然后，按照接受概率

$$\alpha\left(x_j^{\left(i-1\right)},x_j^{'\left(i\right)}|x_{-j}^{\left(i\right)}\right)=\min\left\{1,\frac{p\left(x_j^{'\left(i\right)}|x_{-j}^{\left(i\right)}\right)q\left(x_j^{'\left(i\right)},x_j^{\left(i-1\right)}|x_{-j}^{\left(i\right)}\right)}{p\left(x_j^{\left(i-1\right)}|x_{-j}^{\left(i\right)}\right)q\left(x_j^{\left(i-1\right)},x_j^{'\left(i\right)}|x_{-j}^{\left(i\right)}\right)}\right\}$$

抽样决定是否接受候选值$x_j^{‘\left(i\right)}$。如果接受，则令$x_j^{\left(i\right)}=x_j^{‘\left(i\right)}$；否则，令$x_j^{\left(i\right)}=x_j^{\left(i-1\right)}$。其余分量在第$j$步不改变。

马尔科夫链的转移概率为

$$p\left(x_j^{\left(i-1\right)},x_j^{'\left(i\right)}|x_{-j}^{\left(i\right)}\right)=q\left(x_j^{\left(i-1\right)},x_j^{'\left(i\right)}|x_{-j}^{\left(i\right)}\right)\alpha\left(x_j^{\left(i-1\right)},x_j^{'\left(i\right)}|x_{-j}^{\left(i\right)}\right)$$

由于建议分布可能不被接受，Metropolis-Hastings算法可能在一些相邻的时刻不产生移动。

5 吉布斯抽样

吉布斯抽样用于多元变量联合分布的抽样和评估。

吉布斯抽样是单分量Metropolis-Hastings算法的特殊情况。定义建议分布是当前变量$x_j\left(j=1,2,\cdots,k\right)$的满概率分布

$$q\left(x,x^{'}\right)=p\left(x_j^{'}|x_{-j}\right)$$

则接受概率

$$\begin{align} \alpha\left(x,x^{'}\right) &= \min\left\{1,\frac{p\left(x^{'}\right)q\left(x^{'},x\right)}{p\left(x\right)q\left(x,x^{'}\right)} \right\} \\ &=\min\left\{1, \frac{p\left(x_j^{'}|x_{-j}^{'}\right)p\left(x_{-j}^{'}\right)p\left(x_j|x_{-j}^{'}\right)}{p\left(x_j|x_{-j}\right)p\left(x_{-j}\right)p\left(x_j^{'}|x_{-j}\right)}\right\} \end{align}$$

由于在对第$j$个分量的采样过程中其余分量不变，即$x_{-j}=x_{-j}^{‘}$，所以可得$\alpha\left(x,x^{‘}\right)=1$，即转移概率为$1$。

由以上，可得吉布斯抽样的转移核

$$\begin{align} p\left(x,x^{'}\right) &= q\left(x,x^{'}\right)\alpha\left(x,x^{'}\right) \\ &= p\left(x_j^{'}|x_{-j}\right) \end{align}$$

即一次按照单变量的满条件概率分布$p\left(x_j^{‘}|x_{-j}\right)$进行随机抽样，就能实现但分量Metropolis-Hastings算法。

吉布斯抽样对每次抽样的结果都接受，没有拒绝。抽样会在样本点之间持续移动。

吉布斯抽样算法：
输入：目标概率分布的密度函数$p\left(x\right)$，迭代步数$n$，收敛步数$m$；
输出：$p\left(x\right)$的随机样本$x_{m+1},x_{m+2},\cdots,x_n$；

1. 给出初始样本$x^{\left(0\right)}=\left(x_1^{\left(0\right)},x_2^{\left(0\right)},\cdots,x_k^{\left(0\right)}\right)^\top$
2. 对$i\left(i=1,2,\cdots,n\right)$循环执行
   设第$i-1$轮迭代结束时的样本为$x^{\left(i-1\right)}=\left(x_1^{\left(i-1\right)},x_2^{\left(i-1\right)},\cdots,x_k^{\left(i-1\right)}\right)^\top$
   在第$i$轮迭代中，对$j\left(j=1,2,\cdots,k\right)$循环执行
   $\quad$由满条件分布$p\left(x_j|x_1^{\left(i-1\right)},\cdots,x_{j-1}^{\left(i-1\right)},x_{j+1}^{\left(i-1\right)}\cdots,x_k^{\left(i-1\right)}\right)$抽取$x_j^{\left(i\right)}$
   得到第$i$轮迭代值$x^{\left(i\right)}=\left(x_1^{\left(i\right)},x_2^{\left(i\right)},\cdots,x_k^{\left(i\right)}\right)^\top$
3. 得到样本集合
   $\left\{x^{\left(m+1\right)},x^{\left(m+2\right)},\cdots,x^{\left(n\right)}\right\}$