马尔科夫链 - 蒙特卡洛方法及其一个例子 Metropolis-Hastings 抽样

感情也存在着可以隐瞒的倾斜的程度。

如果有一个随机变量$X$，服从概率密度函数为$\pi(x)$ 的分布。

我们想对这个随机变量抽一组样本出来时，若$\pi(x)$ 比较简单（如均匀分布，正态分布），则显然比较好模拟。但如果$\pi(x)$ 比较复杂，不好模拟怎么办。

这时候，需要构造一个马尔可夫过程，即构造一个马尔可夫转移矩阵$P$（$P_{yx}$ 表示从 x 转移到 y 的概率），使得：

$$\forall y,\pi(y)=\sum_xP_{yx}*\pi(x)$$

矩阵化地，即有$\pi_{n\times 1}=P_{n\times n}*\pi_{n\times 1}$。一个遍历的马尔科夫链有一个很好的性质：对于任意一个分布列$p\in\mathbb{R}^{n\times 1}$ 满足$||p||=1,p_i\geq 0$，都有$\lim_{n\rightarrow\infty}P^np=\pi$。即常说的，“以任意一个状态开始，不断进行马尔可夫转移，最后都会趋于稳定态”。

需要注意的是，马尔科夫链我认为实际上是随机变量间的转移。譬如一系列随机变量为$X_1,X_2,...,X_n$，它们对应的样本值是$x_1,x_2,...,x_n$。

• 若$X_1\sim U[0,1]$，我们很容易地随机取了个$x_1$ 出来。
• 然后再根据转移矩阵，以概率分布$P_{1x_1},P_{2x_1},P_{3x_1},...,P_{nx_1}$ 随机抽取得$x_2$，则显然$x_2$ 对应的随机变量$X_2$ 的概率分布就不是$U[0,1]$ 了。
• 而重复这个过程，当趋于无穷时，我们就会发现$\lim_{n\rightarrow\infty}X_n\sim\pi$。

那么，其实$X_1,X_2,...,X_n$ 就服从分布$U[0,1],P*U[0,1],...,\pi,\pi,\pi,...$，因此我们考虑足够大的一些样本值，它们都是同分布于$\pi$ 的，即做到了抽取一组服从分布$\pi$ 的样本值。

---

然而如何得到$P$ 就是一个很大的问题。我发现很多地方都写的$P$ 是依赖于$\pi$ 的，那既然都可以模拟$P$ 了为啥不直接模拟$\pi$ 呢？

但其实，其简化模拟的核心在于，抽样$\pi$ 需要知道$\pi$ 的连续信息，而马尔科夫链的转移只需要知道$\pi$ 的点信息。

考虑 Metropolis-Hastings 算法：

• 令$Q_{yx}$ 是一个好模拟的，简单的分布。譬如$Q_{yx}=U[0,1]$ 与 x 无关，或$Q_{yx}=G(x,0.01)$ 即以 x 为均值的正态分布。
• 令$\alpha(y,x)=\pi(x)*Q_{xy}$，则可以构造$P_{yx}=Q_{yx}*\alpha(y,x)$。

这样做有什么意义呢？注意到在已知样本$x_1$，想要得到下一个样本$x_2$ 时，我们的操作并不是把$P_{1x_1},P_{2x_2},...$ 都算出来，然后按照分布抽取。而应该先按照$Q_{yx_1}$ 抽取出来一个$y$，再计算$\alpha(y,x_1)$，然后再以$\alpha(y,x_1)$ 的概率转移到$x_2=y$。以$1-\alpha(y,x)$ 的概率停留在$x_2=x_1$。这就是一个简单分布 + 两点分布解决了转移的过程。简单分布负责给出候选值，再计算$\alpha$，再根据两点分布决定是否转移。而最终转移的概率就是$P_{yx}$。

于是这样做，显然就只需要知道$\pi$ 的点信息。

然而这么做为什么是对的？

首先，$x_2$ 停留在$x_1$ 的概率显然是$\sum_yQ_{yx}(1-\alpha(y,x))=1-\sum_{y}Q_{yx}\alpha(y,x)=1-\sum_y P_{yx}$ 这很合理。

其次，这样构造的$P$ 满足了$P_{yx}*\pi(x)\equiv P_{xy}*\pi(y)$。而一个重要的结论就是，满足这个恒等式（detailed balance）的马尔科夫链是遍历的，即其存在一个稳定的输出$\pi(x)$，即我们需要抽样的随机分布。

---

所以说，我认为 MCMC 之所以可以模拟抽样非常复杂的分布，一个重要内容就是把复杂分布的连续信息，变成了单点 check 的信息。我通过一个简单的分布给出一个候选值，再根据$\pi(x)$ 来决定是否采用，而不是根据$\pi(x)$ 来得到候选值。

# 一个例子 —— 指数分布

譬如我想抽取满足概率密度函数$f(x)$ 的随机变量$X$ 的一组样本:

$$f(x)=\begin{cases}e^{-x}&x\geq 0\\ 0&else\end{cases}$$

很显然它是连续的，无穷的，想进行抽样非常麻烦。

但我们可以令$x_1=10$，然后在$U[10-3,10+3]$ 中随机抽取一个数，譬如$10.3$。

那么计算$\alpha(10.3,10)=f(10)*\frac{1}{6}=\frac{e^{-10}}{6}$。然后以这个概率进行二项分布抽样，譬如抽出来选择不接受。

则$x_2=x_1=10$。

简单分析一下就会发现，这个过程如果重复足够多次，$x$ 是会大概率不断减小的。这就是指数分布在$x$ 较小时概率较大的原因。由于我们初始值$x_1=10$ 这个取得不是很好，因此如果想要逼近稳定输出$f(x)$，需要重复很多很多次。

譬如重复了 10000 次后，我们再取出$x_{10000},...,x_{10010}$ 这 11 个样本值，就可以当作指数分布的一组样本了。