MCMC简介

MCMC，第一个MC是马尔可夫(Markov Chain)，指的是从分布中抽样，部分，第二个MC是蒙特卡罗(Monte Carlo)，指的是用蒙特卡罗方法计算。Metropolis-Hastings算法是MCMC的代表算法。

MCMC的理论依据是如下的遍历定理：

定理：设有马尔科夫链X=\{X _ 0,X _ 1,X _ 2,\dots,X _ t,\dots\}，有离散状态空间为S（连续状态空间表述完全类似，不再重复），若马尔科夫链X是不可约的、非周期的、正常返的，则：X有唯一平稳分布\pi=\{\pi _ 1,\pi _ 2,\dots\}^\mathsf T，且转移概率的极限分布是这个平稳分布，即
\mathbb P(X _ t=i\mid X _ 0=j)\to\pi _ i, \quad\forall i,j\in\mathbb N^*.

（我们发现，和初始分布无关）而且，若f是定义在状态空间上的函数、\mathbb E|f|<\infty，则：以概率1成立
\frac1t\sum _ {s=1}^t f(x _ s)\to\mathbb E _ {x\sim\pi}f(x).
对于要随机抽样或者求数学期望的情形，假设要抽样的分布是p(\boldsymbol x)，MCMC的基本想法是：在随机变量\boldsymbol x的状态空间S上定义一个满足遍历定理的马尔可夫链，使其平稳分布就是抽样的目标分布p(\boldsymbol x)。然后在这个马尔可夫链上进行随机游走，每个时刻得到一个样本。根据遍历定理，当时间趋于无穷时，样本的分布趋于平稳分布，样本的函数均值趋近函数的数学期望。所以，当时间足够长时（时刻t\geqslant m），在之后的时间（时刻t\leqslant n）里随机游走得到的样本\{\boldsymbol x _ {m+1},\boldsymbol x _ {m+2},\dots,\boldsymbol x _ n\}就是目标分布的抽样结果，得到的函数均值就是需要的值：
\widehat {\mathbb E}f=\frac1{n-m}\sum _ {s=m+1}^nf(\boldsymbol x _ s).

连续状态马尔可夫链

如果马尔可夫链定义在连续状态空间上，那么转移概率分布由转移核表示。

若对任意的\boldsymbol x，定义在S上\sigma-代数的集合函数P(\boldsymbol x,\cdot)是概率测度，那么就称之为转移核。转移核P(\boldsymbol x,A)表示从\boldsymbol x转移到属于A的状态的概率。如果P(\boldsymbol x,\cdot)有密度p(\boldsymbol x,\cdot)，那么也可将p(\boldsymbol x,\cdot)称为转移核。

分布\pi(\boldsymbol x)称为平稳分布，如果
\pi(y)=\int _ S p(\boldsymbol x,\boldsymbol y)\pi(\boldsymbol x)\,\mathrm d\boldsymbol x,\quad\forall \boldsymbol y\in S,

简写为\pi=P\pi（和离散时的定义统一了）。

Metropolis-Hastings算法

如何构建具体的马尔可夫链成为这个方法的关键。Metropolis-Hastings算法可描述如下：任选初值后，假设已有i-1个样本，设状态\boldsymbol x _ {i-1}=\boldsymbol x，按照提议分布q(\boldsymbol x,\boldsymbol x^\prime)（提议分布是另一个马尔可夫链的转移核，便于采样）采样，得到一个候选状态\boldsymbol x^\prime。然后计算接受概率
\alpha(\boldsymbol x,\boldsymbol x^\prime )=\min\left\{1,\frac{p(\boldsymbol x^\prime )q(\boldsymbol x^\prime ,\boldsymbol x)}{p(\boldsymbol x)q(\boldsymbol x,\boldsymbol x^\prime )}\right\}.

再以\alpha(\boldsymbol x,\boldsymbol x^\prime )的概率将状态转移到\boldsymbol x^\prime，操作是，从区间(0,1)中随机抽取一个数u，若u\leqslant \alpha(\boldsymbol x,\boldsymbol x^\prime )，则状态\boldsymbol x _ i=\boldsymbol x^\prime；否则，状态\boldsymbol x _ i=\boldsymbol x。重复上述步骤，得到\{\boldsymbol x _ {m+1},\boldsymbol x _ {m+2},\dots,\boldsymbol x _ n\}后，按上面的式子计算\widehat {\mathbb E}f。

可以证明，这个算法对应于转移核为p(\boldsymbol x,\boldsymbol x^\prime )=q(\boldsymbol x,\boldsymbol x^\prime )\alpha(\boldsymbol x,\boldsymbol x^\prime )+\mathbb I(\boldsymbol x^\prime =\boldsymbol x)[\int(1-\alpha(\boldsymbol x,\boldsymbol x^\prime ))p(\boldsymbol x,\boldsymbol x^\prime )\,\mathrm d\boldsymbol x^\prime ]的马尔可夫链，这个马尔可夫链是可逆马尔可夫链，因而满足遍历定理；且其平稳分布就是p(\boldsymbol x)，即要抽样的分布，这就具体实现了MCMC。

证明

证明分两部分，第一部分我们证明转移核为p(\boldsymbol x,\boldsymbol x^\prime )的马尔可夫链是可逆的，即我们有细致平衡方程成立：
p(\boldsymbol x)p(\boldsymbol x,\boldsymbol x^\prime )=p(\boldsymbol x^\prime )p(\boldsymbol x^\prime ,\boldsymbol x),

并且p(\boldsymbol x)是该马尔可夫链的平稳分布。如果能证明，那么由上面的定理，这个平稳分布正是该马尔可夫链的唯一平稳分布，就可以应用遍历定理。

当\boldsymbol x=\boldsymbol x^\prime，上式显然成立。若不然，
\begin{aligned}p(\boldsymbol x) p\left(\boldsymbol x, \boldsymbol x^{\prime}\right) &=p(\boldsymbol x) q\left(\boldsymbol x,\boldsymbol x^{\prime}\right) \min \left\{1, \frac{p\left(\boldsymbol x^{\prime}\right) q\left(\boldsymbol x^{\prime}, \boldsymbol x\right)}{p(\boldsymbol x) q\left(\boldsymbol x, \boldsymbol x^{\prime}\right)}\right\} \\&=\min \left\{p(\boldsymbol x) q\left(\boldsymbol x, \boldsymbol x^{\prime}\right), p\left(\boldsymbol x^{\prime}\right) q\left(\boldsymbol x^{\prime}, \boldsymbol x\right)\right\} \\&=p\left(\boldsymbol x^{\prime}\right) q\left(\boldsymbol x^{\prime}, \boldsymbol x\right) \min \left\{\frac{p(\boldsymbol x) q\left(\boldsymbol x, \boldsymbol x^{\prime}\right)}{p\left(\boldsymbol x^{\prime}\right) q\left(\boldsymbol x^{\prime}, \boldsymbol x\right)}, 1\right\} \\&=p\left(\boldsymbol x^{\prime}\right) p\left(\boldsymbol x^{\prime}, \boldsymbol x\right).\end{aligned}

仍然成立。

验证是平稳分布：
\begin{aligned}\int _ S p(\boldsymbol x) p\left(\boldsymbol x, \boldsymbol x^{\prime}\right) \mathrm{d} x &=\int _ S p\left(\boldsymbol x^{\prime}\right) p\left(\boldsymbol x^{\prime}, \boldsymbol x\right) \mathrm{d} x \\&=p\left(\boldsymbol x^{\prime}\right) \int _ S p\left(\boldsymbol x^{\prime}, \boldsymbol x\right) \mathrm{d} x \\&=p\left(\boldsymbol x^{\prime}\right).\end{aligned}

对照前面的定义式即知确实是平稳分布。

下面证明第二部分，按照Metropolis-Hastings算法得到的马尔可夫链的转移核就是
p(\boldsymbol x,\boldsymbol x^\prime )=q(\boldsymbol x,\boldsymbol x^\prime )\alpha(\boldsymbol x,\boldsymbol x^\prime )+\mathbb I(\boldsymbol x^\prime =\boldsymbol x)\Big[\int _ S(1-\alpha(\boldsymbol x,\boldsymbol x^\prime ))q(\boldsymbol x,\boldsymbol x^\prime )\,\mathrm d\boldsymbol x^\prime \Big].

就是说，对A\subseteq S，有
\mathbb P(\boldsymbol x\to A)=\int _ A p(\boldsymbol x,\boldsymbol x^\prime )\,\mathrm d\boldsymbol x^\prime .
考虑两种基本情况，第一种\boldsymbol x\notin A，第二种\{\boldsymbol x\}=A。若\boldsymbol x\notin A，要将状态从\boldsymbol x转移到\boldsymbol x^\prime，首先采样\boldsymbol y\sim q(\boldsymbol x,\cdot)，然后以概率\alpha(\boldsymbol x,\boldsymbol y)转移状态：\boldsymbol x^\prime =\boldsymbol y（现在的叙述又引入了\boldsymbol y，避免混淆），否则\boldsymbol x^\prime =\boldsymbol x。我们发现这正类似一次拒绝采样法。类似拒绝采样中的证明，
\begin{aligned}\mathbb P(\boldsymbol x^\prime \in A)&=\mathbb P(\boldsymbol y\neq \boldsymbol x,u\leqslant\alpha(\boldsymbol x,\boldsymbol y),\boldsymbol y\in A)\\&=\int _ A\int _ 0^{\alpha(\boldsymbol x,\boldsymbol y)}q(\boldsymbol x,\boldsymbol y)\, \mathrm du\mathrm d\boldsymbol y\\&=\int _ Aq(\boldsymbol x,\boldsymbol y)\alpha(\boldsymbol x,\boldsymbol y)\, \mathrm d\boldsymbol y.\end{aligned}
对第二种情况，
\begin{aligned}\mathbb P(\boldsymbol x\to\boldsymbol x)&=\mathbb P(\boldsymbol y=\boldsymbol x)+\mathbb P(\boldsymbol x^\prime =\boldsymbol x,\, \boldsymbol y\neq\boldsymbol x)\\&=q(\boldsymbol x,\boldsymbol x)+\int _ {S\backslash\{\boldsymbol x\}}\int _ {\alpha(\boldsymbol x,\boldsymbol y)}^1q(\boldsymbol x,\boldsymbol y)\, \mathrm du\mathrm d\boldsymbol y\\&=q(\boldsymbol x,\boldsymbol x)+\int _ S(1-\alpha(\boldsymbol x,\boldsymbol y))q(\boldsymbol x,\boldsymbol y)\,\mathrm d\boldsymbol y\\&=q(\boldsymbol x,\boldsymbol x)\alpha(\boldsymbol x,\boldsymbol x)+\mathbb I(\boldsymbol x=\boldsymbol x)\Big[\int _ S(1-\alpha(\boldsymbol x,\boldsymbol x^\prime ))q(\boldsymbol x,\boldsymbol x^\prime )\,\mathrm d\boldsymbol x^\prime \Big].\end{aligned}
一般情况则可以由两种情况组合得到。这就完成了第二部分的证明。

另证

前面提到，和拒绝采样不同的是不会像拒绝采样法那样丢弃样本，而是将采样结果改为\boldsymbol x，所以可以得知条件分布p(\boldsymbol x^\prime \mid\boldsymbol x)在\boldsymbol x点处会有概率质量。但对于非\boldsymbol x点处，还是有类似拒绝采样中的结论，即
p(\boldsymbol x,\boldsymbol x^\prime )\propto q(\boldsymbol x,\boldsymbol x^\prime )\alpha(\boldsymbol x,\boldsymbol x^\prime ),\quad\forall\boldsymbol x^\prime\neq\boldsymbol x.

但是我们需要的是等于号，要成立等号这时候只需要再考虑证明\boldsymbol x处
p(\boldsymbol x,\boldsymbol x)=1-\int _ {S\backslash\{\boldsymbol x\}}q(\boldsymbol x,\boldsymbol x^\prime )\alpha(\boldsymbol x,\boldsymbol x^\prime )\, \mathrm d\boldsymbol x^\prime .