统计模拟中有一个重要的问题就是，对于给定的概率分布 $p(x)$ ，如何生成它对应的样本，常见的概率分布都可以基于 $Uniform(0, 1)$ 的样本生成。但当 $p(x)$ 形式复杂或是一个高维分布的时候，就需要请出MCMC（Markov Chain Monte Carlo）方法了。
0 引理
在请出MCMC方法之前，需要先认识一下MCMC的基石：马氏链收敛定理和细致平稳条件。
马氏链收敛定理 若一个非周期马氏链具有转移概率矩阵$P$，且它的任何两个状态是连通的，那么 $\lim_{n \to \infty}{P_{ij}^n}$ 存在且与 $i$ 无关，记 $\lim_{n \to \infty}{P_{ij}^n} = \pi(j)$ ，我们有：
- \[\lim_{n \to \infty}{P^n} = \left[ \begin{matrix} \pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \\ \pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ \pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots \end{matrix} \right]\]
-
$\pi(j) = \sum_{i=0}^{\infty}{\pi(i)P_{ij}}$
- $\pi$是方程 $\pi{P} = \pi$ 的唯一非负解，其中 $\pi = [\pi(1), \pi(2), \cdots, \pi(j), \cdots], 且\sum_{i=0}^{\infty}{\pi(i)} = 1$ ，这里 $\pi$ 称为马氏链的平稳分布。
细致平稳条件 若非周期马氏链的转移矩阵 $P$ 和分布 $\pi(x)$ 满足
\[\pi(i)P_{ij} = \pi(j)P_{ji} \quad for\:all\;i,j\]则 $\pi(x)$ 是马氏链的平稳分布，上式被称为细致平稳条件。
1 MCMC采样算法
假设我们有一个转移矩阵为 $Q$ 的马氏链（ $q(i,j)$ 为 $Q$ 的元素，代表从状态 $i$ 转移到状态 $j$ 的概率），通常情况下：
\[p(i)q(i,j) \neq p(j)q(j,i)\]也就是细致平稳条件不成立，因而$p(x)$不太可能为马氏链的平稳分布。为了使细致平稳条件成立，我们引入$\alpha(i,j)$，希望
\[p(i)\underbrace{q(i,j)\alpha(i,j)}_{q'(i,j)} = p(j)\underbrace{q(j,i)\alpha(j,i)}_{q'(j,i)} \tag{1}\]我们不妨取
\[\alpha(i,j) = p(j)q(j,i),\quad \alpha(j,i) = p(i)q(i,j)\]于是(1)式成立了。这样，把原来具有转移矩阵 $Q$ 的普通马氏链改造成了具有转移矩阵 $Q’$ （ $q’(i,j)$ 为 $Q’$ 的元素，代表从状态 $i$ 转移到状态 $j$ 的概率）的马氏链，而 $Q’$ 恰好满足细致平稳条件，因而 $Q’$ 对应的分布 $p(x)$ 就是平稳分布。
在改造 $Q$ 过程中引入的 $\alpha(i,j)$ 称为接受率，可以理解为在原来的马氏链上，从状态 $i$ 以 $q(i,j)$ 的概率转化为状态 $j$ 时，我们以 $\alpha(i,j)$ 的概率接受这个转移，于是得到了新的马氏链 $Q’$ 的转移概率为 $q(i,j)\alpha(i,j)$ 。如此，我们就有了用于采样概率分布 $p(x)$ 的算法。
该过程中讨论的是离散的情况，事实上即使分布式连续的，MCMC算法仍然有效，只不过 $q(x\vert y)$ 就是连续二元概率分布对应的条件分布。
2 M-H采样算法
上一节的算法虽然已经能够采样任何分布了，但存在一个小问题：效率不高。原因是采样过程中接受率 $\alpha(i,j)$ 可能偏小，使得马氏链拒绝大量跳转，因而马氏链遍历所有状态空间所花费的时间可能很长，导致收敛到平稳分布的速度太慢。如何提升一下接受率呢？
我们看一个例子，假设 $\alpha(i,j) = 0.1, :\alpha(j,i) = 0.2$ 时满足细致平稳条件，于是
\[p(i)q(i,j)\times 0.1 = p(j)q(j,i)\times 0.2\]为了提高接受率，我们给式子两边同乘0.5，得到
\[p(i)q(i,j)\times 0.5 = p(j)q(j,i)\times 1\]此时，细致平稳条件依旧成立！因而我们可以将两边同时放大，使得最大的一个放大到1，即取
\[\alpha(i,j) = \min{\lbrace\frac{p(i)q(i,j)}{p(j)q(j,i)}, 1\rbrace}\]这样就得到了最常见的Metroplis-Hastings算法，简称MH算法。
此处，$x$并不要求是一维的，对于高维空间的 $p(x)$ ，若满足细致平稳条件
\[p(x)q'(x \to y) = p(y)q'(y \to x)\]则M-H算法依旧适用。
3 Gibbs采样算法
对于高维的情形，由于接受率 $\alpha$ 的存在（通常 $\alpha < 1$ ），MH算法的效率不够高，能否找到一个转移矩阵 $Q$ 使得接受率 $\alpha = 1$ 呢？我们先看看二维的情形，假设有一个概率分布 $p(x.y)$ ，考察$x$坐标相等的两个点 $A(x_1, y_1), B(x_1, y_2)$ ，会发现
\[\begin{equation} p(x_1, y_1)p(y_2\vert x_1) = p(x_1)p(y_1 \vert x_1)p(y_2 \vert x_1) \\ p(x_1, y_2)p(y_1 \vert x_1) = p(x_1)p(y_2 \vert x_1)p(y_1 \vert x_1) \end{equation}\]因而得到
\[p(x_1, y_1)p(y_2 \vert x_1) = p(x_1, y_2)p(y_1 \vert x_1)\]即
\[p(A)q(y_2 \vert x_1) = p(B)p(y_1 \vert x_1)\]从上面的式子中我们可以发现，在 $x = x_1$ 这条平行于 $y$ 轴的直线上，如果使用条件分布 $p(y \vert x_1)$ 作为任何两个点之间的转移概率，那么任何两个点之间的转移满足细致平稳条件。同样的，如果我们在 $y = y_1$ 这条直线上任意取两点 $A(x_1, y_1), C(x_2, y_1)$ ，也有如下等式
\[p(A)p(x_2 \vert y_1) = p(C)p(x_1 \vert y_1)\]这样我们就可以构造平面上两点之间的转移矩阵概率 $Q$
\[\begin{cases} Q(A \to B) = p(y_B \vert x_1), & x_A = x_B = x_1 \\ Q(A \to C) = p(x_C \vert y_1), & y_A = y_C = y_1 \\ Q(A \to D) = 0, & 其它 \end{cases}\]很容易验证该转移矩阵 $Q$ 对于平面上的任意两点 $X, Y$ 满足细致平稳条件，即
\[p(X)Q(X \to Y) = p(Y)Q(Y \to X)\]于是这个二维空间上的马氏链将收敛到平稳分布 $p(x, y)$ ，而这个算法就成为Gibbs采样算法。
以上采样过程如图所示：
马氏链的转移只是沿着坐标轴轮换采样，随之产生了样本 $(x_0, y_0), (x_0, y_1), (x_1, y_1), (x_1, y_2), \cdots$ 马氏链收敛后，最终得到的就是 $p(x, y)$ 的样本，而收敛之前的阶段称为burn-in period。其实，坐标轴轮换只是为了方便，在 $t$ 时刻，也可以从各坐标轴之间随机选一个，然后按照条件概率进行转移，最终亦可收敛。
很容易将其推广到多维，得到高铝分布 $p(x_1, x_2, \cdots, x_n)$ 的样本，值得注意的是这些样本并不独立，但此处要求的仅是采样得到符合给定的概率分布的样本，并不要求独立。
4 SMC采样算法
传统的MCMC算法如MH算法很难表征多模式的准后验分布，如神经网络的参数，当遇到这种分布时，由于没有足够的信息建立模式之间的连接，马氏链很可能困在局部模式中，相比之下，SMC算法并行地从先验中抽样，这些原始样本配合着准则信息在调和分布上传播，最后得到准后验。
在算法过程中，随着新的信息被添加进来，算法会舍弃一些质量相对低的抽样，复制一些质量相对高的，接着在这些样本的基础上通过MCMC方法产生新样本。
算法如下：
在实践中，可以使用 $\phi_j = (\frac{j-1}{J-1})^\lambda$ 来初始化 $\phi$，通常选取 $\lambda = 2$ 。
可以看到，在筛选和变异的时候，我们能够并行地处理运算，如此一来，相较于传统的MCMC方法，其余步骤花费的时间就不值一提了。
