Markov Chain Monte Carlo（1）

最近因为种种原因不得已去钻研一下一些基础的机器学习的知识。当初统计和矩阵论的知识都忘得差不多了。好在我爱自然语言处理这里的文章写得都比较清楚，剩下不懂的自己写写算算和朋友讨论讨论，也有一些收获。大部分都是转载，附加了一些学习过程中自己的笔记（免得将来忘了又得从头思考）

---

随机模拟(或者统计模拟)方法有一个很酷的别名是蒙特卡罗方法(Monte Carlo Simulation)。统计模拟中有一个重要的问题就是给定一个概率分布$p(x)$，我们如何在计算机中生成它的样本。一般而言均匀分布$Uniform(0,1)$的样本是相对容易生成的。 通过线性同余发生器可以生成伪随机数，我们用确定性算法生成$[0,1]$之间的伪随机数序列后，这些序列的各种统计指标和均匀分布$Uniform(0,1)$ 的理论计算结果非常接近。这样的伪随机序列就有比较好的统计性质，可以被当成真实的随机数使用。

而我们常见的概率分布，无论是连续的还是离散的分布，都可以基于$Uniform(0,1)$的样本生成。例如正态分布可以通过著名的 Box-Muller 变换得到。不过我们并不是总是这么幸运的，当$p(x)$的形式很复杂，或者 $p(x)$是个高维的分布的时候，样本的生成就可能很困难了。此时就需要使用一些更加复杂的随机模拟的方法来生成样本。而本节中将要重点介绍的 MCMC(Markov Chain Monte Carlo) 和 Gibbs Sampling算法就是最常用的一种，这两个方法在现代贝叶斯分析中被广泛使用。要了解这两个算法，我们首先要对马氏链的平稳分布的性质有基本的认识。

马氏链及其平稳分布

马氏链的数学定义很简单，也就是状态转移的概率只依赖于前一个状态。
$P(X_{t+1}=x|X_t,X_{t−1},⋯)=P(X_{t+1}=x|X_t)$

假设有初始分布$\pi_0=(\pi_1, \pi_2, ... \pi_m)$和状态转移矩阵$P= \bigl( \begin{smallmatrix} p_{11} & p_{12} & \cdots & p_{1m} \\ p_{21} & p_{22} & \cdots & p_{2m} \\ \vdots & \vdots & \ddots & \vdots \\ p_{m1} & p_{m2} & \cdots & p_{mm} \end{smallmatrix} \bigr)$ ，那么有
马氏链定理： 如果一个非周期马氏链具有转移概率矩阵$P$,且它的任何两个状态是连通的，那么 $\lim_{n \to \infty}P_{ij}^n$存在且与$i$无关（这里$P_{ij}^n$表示矩阵$P$的$n$次幂的$i$行$j$列的元素），则$lim_{n \to \infty}P_{ij}^n=\pi(j)$, 并且我们有：

1. $\pi P=\pi$，其中$\pi=\lim_{n \to \infty}\pi_0 P^n$且$\pi$与$\pi_0$无关；记$\pi=[\pi(0), \pi(1), \cdots \pi(m)]$，$\pi$称为马氏链的平稳分布。
2. $\pi$是方程$\pi P=\pi$的唯一非负解。
3. $\pi(j)=\sum_{i=0}^m\pi(i)P_{ij}$
4. 
   $$\lim_{n\to \infty}P=\left[\begin{matrix} \pi(0) & \pi(1) & \cdots & \pi(m) \\ \pi(0) & \pi(1) & \cdots & \pi(m) \\ \vdots & \vdots & \ddots & \vdots \\ \pi(0) & \pi(1) & \cdots & \pi(m) \\ \end{matrix} \right]$$

这个马氏链的收敛定理非常重要，所有的 MCMC(Markov Chain Monte Carlo) 方法都是以这个定理作为理论基础的。 定理的证明相对复杂，一般的随机过程课本中也不给证明，所以我们就不用纠结它的证明了，直接用这个定理的结论就好了。就定义3做一些说明。
由矩阵乘法，$\pi_n(j)=\sum_{i=0}^n\pi_{n-1}(i)p_{ij}$，对等式两侧n取极限有：
$\lim_{n\to \infty}\pi_n(j)=\sum_{i=0}^n\lim_{n\to \infty}\pi_{n-1}(i) p_{ij}$
由于$\lim_{n\to \infty}\pi_n(j)=\pi(j)$，故：$\pi(j)=\sum_{i=0}^m\pi(i)P_{ij}$