蒙特卡洛方法

常见使用场景

机器学习中经常会遇到对复杂的分布做加和或积分，例如在贝叶斯方法中，往往要对参数做积分，\(P(t|X)=\int p(t|\theta)p(\theta|X)d\theta\)，频率派中EM算法的E步也是一个求期望的过程，\(Q(\theta,\theta_{old})=\int p(Z|X,\theta_{old})log\,p(Z,X|\theta)dZ\)，这些积分或者期望往往都是intractable的。对于这个问题，我们可以使用变分来解决，变分是通过最优化方法寻找一个和\(p(\theta|X)\)或\(p(Z|X,\theta_{old})\)近似且好积分的分布。而在这篇博文中，我们将介绍sampling的方法，举个简单的例子：比如我们遇到这种形式\(E[f]=\int f(z)p(z)\,dz\)，如果我们从p(z)中sampling一个数据集\(z^{(l)}\)，然后再求个平均\(\hat f=\frac{1}{L}\sum_{l=1}^L f(z^{(l)})\)来近似\(f(z)\)的期望，so问题就解决了，关键是如何从\(p(z)\)中做无偏的sampling。现在的问题转换为给定一个概率分布，我们如何在计算机中生成它的样本。一般而言均匀分布的样本是相对容易生成的。 通过线性同余发生器可以生成伪随机数，这些伪随机数序列的各种统计指标和均匀分布的理论计算结果非常接近。这样的伪随机序列就有比较好的统计性质，可以被当成真实的随机数使用。而我们常见的概率分布，无论是连续的还是离散的分布，都可以基于均匀分布的样本生成。

蒙特卡洛方法

1.Inverse transform sampling

我们知道，计算机本身是无法产生真正的随机数的，但是可以根据一定的算法产生伪随机数。最古老最简单的莫过于 线性同余发生器：$$x_{n+1}=(ax_{n}+c);mod;m$$式子中的\(a\)和\(c\)是一些数学知识推导出的合适的常数。但是我们看到，这种算法产生的下一个随机数完全依赖现在的随机数的大小，而且当你的随机数序列足够大的时候，随机数将出现重复子序列的情况。当然，理论发展到今天，有很多更加先进的随机数产生算法出现，比如python数值运算库numpy用的是Mersenne Twister等。但是不管算法如何发展，这些都不是本质上的随机数，用冯诺依曼的一句话说就是：

Anyone who considers arithmetic methods of producing random digits is, of course, in a state of sin.

OK，根据上面的算法现在我们有了均匀分布的随机数，但是如何产生满足其他分布（比如高斯分布）下的随机数呢？一种可选的简单的方法是 Inverse transform sampling。它的原理是利用累积分布函数（CDF，cumulative distribution function）来处理。

• Inverse transform sampling
• 假设：已经一个生成器可以产生(0, 1)区间上的均匀分布随机数
• 问题：从分布p(z)中采样
• 算法：假设有分布\(p(x)\)的累计函数是\(P(x)\)，若\(y\)是(0, 1)上的均匀分布随机变量，那么\(P^{-1}(y)\)的就是服从分布\(p(x)\)的随机变量
• 缺陷：\(P^{-1}(y)\)不容易计算

举个列子，我们想对高斯分布采样，首先在 y 轴上产生（0,1）之间的均匀分布的随机数，水平向右投影到高斯累计分布函数上，然后垂直向下投影到 x 轴，得到的就是高斯分布的一个实例。可见高斯分布的随机数实际就是均匀分布随机数在高斯分布的 CDF 函数下的逆映射。当然，在实际操作中，更有效的计算方法有Box-Muller_transform ，Ziggurat algorithm 等，这些方法 tricky and faster，没有深入了解，这里也不多说了。

2.Rejection Sampling

Rejection Sampling和Importance Sampling都是基于proposal distribution的sampling，这两种方法都是假设直接从分布p(z)采样很困难，但可以找到一个相对更容易采用的分布q(z)，即proposal distribution。
假设分布\(p(z)=\frac{1}{Z}\hat p(z)\)，其中Z是\(p(z)\)中与z无关的一个因子。之所以要分离开来写，是因为有时候\(p(z)\)中与z无关的部分可能比较复杂，难以计算。例如\(p(z|x)=\frac{p(x,z)}{p(x)},p(x)=\int p(x,z)dz\)，其中，\(p(x)\)与z无关且难以计算。因此，在大多数的sampling中，只需要利用与z有关部分即可进行采样的算法，这样就避开了对复杂的Z的计算。

• Rejection Sampling
• 假设：直接从分布p(z)中采样困难
• 问题：从分布p(z)中采样
• 算法：首先为p(z)找一个proposal distribution q(z)，而且q(z)必须是很容易进行sampling的。然后找到一个尽可能小的常数k，使得\(kq(z)\geq p(z)\)，对任意z成立。然后从q(z)中sample出一个数\(z_0\)；然后从均匀分布[0,\(kq(z_0)\)]中sample出另一个数\(u_0\)；这时候，平面上的点（\(z_0\)，\(u_0\)）是\(kq(z)\)下方区域中的均匀分布。如果\(u_0>\hat p(z_0)\)，则拒绝样本\(z_0\)并重复前面步骤，否则接收\(z_0\)为符合分布p(z)的点。
• 缺陷：上述过程看出了k为什么要尽可能小。K越小，才能使\(z_0\)被拒绝的概率尽可能小，从而提高rejection sampling的效率，因此我们需要选择一个合适的k 。维数越高，拒绝率越高，采样效率越低。例如高维的球，可计算其测度主要集中在球的表面；而rejection sampling中，\(u_0>\hat p(z_0)\)的部分正是高维几何体的表层。这就导致在高维情况下，有很高的拒绝率。

z从q(z)分布抽样，而接受率是\(\frac{\hat p(z)}{kq(z)}\)，所以

\begin{split}
p(accept)&=\int {\hat p(z)/kq(z)}q(z)dz \
&= \frac{1}{k} \int \hat p(z)dz
\end{split}

3.Importance Sampling

Importance Sampling和Rejection Sampling类似都是假设对p(z)采样比较困难，不过对于一个给定的z，却可容易的计算其概率值p(z)。但是与Rejection Sampling不同的是，Importance Sampling不是求p(z)样例，而是直接计算函数f(z)在该分布下的期望。因为p(z)本身采样困难，所以我们还是得像Rejection sampling那样，找到一个更容易采样的分布q(z)，并且假设从q(z)采样了L个样本。那么：$$E[f(z)]=\int f(z)p(z)dz=\int \frac{f(z)p(z)}{q(z)}q(z)dz\approx \frac{1}{L}\sum_{l=1}^L\frac{p(z^l)}{q(z^l)}f(z^l) $$
其中\(r_l=\frac{p(z^l)}{q(z^l)}\)被成为importance weights。
再进一步假设：对分布p(z)的认识集中在p(z)的与z相关的部分\(\hat p(z)\)，其normalization constant \(Z_p\)还未知。同时也从q(z)中分离出一个常数\(Z_q\)，那么：

\begin{split}
E[f(z)]&=\int f(z)p(z)dz=\frac{Z_q}{Z_p}\int \frac{f(z)\hat p(z)}{\hat q(z)}q(z)dz\
&\approx \frac{Z_q}{Z_p}\frac{1}{L}\sum_{l=1}^L\frac{\hat p(z^l)}{\hat q(z^l)}f(z^l)
\end{split}

其中\(\hat r_l=\frac{\hat p(z^l)}{\hat q(z^l)}\)，同样地，可以计算：

\[\frac{Z_p}{Z_q}=\frac{1}{Z_q}\int \hat p(z)dz=\int \frac{\hat p(z)}{\hat q(z)}q(z)dz\approx \frac{1}{L} \sum_{l=1}^L\hat r_l \]

于是最终得到：

\[E[f(z)]=\sum_{l=1}^Lw_lf(z^l) \]

其中，\(z^l\)是分布从q(z)采样的L个样本，而

\[w_l=\frac{\hat r^l}{\sum_{m=1}^L\hat r^m} \]

注意，在\(w_l\)的计算中，已经只需要分布p(z)的与z有关部分\(\hat p(z)\)。从而达到了避开Z的目的。

• Importance Sampling
• 假设：直接从分布p(z)中采样困难
• 问题：计算函数f(z)在p(z)分布下的期望
• 算法：首先为p(z)找一个proposal distribution q(z)，而且q(z)必须是很容易进行sampling的并且和目标分布相似。然后从q(z)中sample出L个数\(z_l\)，带入公式\(E[f(z)]=\sum_{l=1}^Lw_lf(z^l)\)求的f(z)在p(z)分布下的期望
• 缺陷：但是可惜的是，在高维空间里找到一个这样合适的q非常难。即使有 Adaptive importance sampling 和 Sampling-Importance-Resampling(SIR) 的出现，要找到一个同时满足容易抽样并且和目标分布相似的proposal distribution，通常是不可能的！

马尔可夫蒙特卡洛

由于Rejection sampling和Importance sampling这两种方法在高维下都会失效。我们的目标还是对于给定的概率分布p(x),我们希望能有便捷的方式生成它对应的样本。由于马氏链能收敛到平稳分布，于是一个很的漂亮想法是：如果我们能构造一个转移矩阵为P的马氏链，使得该马氏链的平稳分布恰好是p(x), 那么我们从任何一个初始状态x0出发沿着马氏链转移, 得到一个转移序列 x0,x1,x2,⋯xn,xn+1⋯,， 如果马氏链在第n步已经收敛了，于是我们就得到了样本xn,xn+1⋯。

这个绝妙的想法在1953年被Metropolis想到了，为了研究粒子系统的平稳性质，Metropolis考虑了物理学中常见的波尔兹曼分布的采样问题，首次提出了基于马氏链的蒙特卡罗方法，即Metropolis算法，并在最早的计算机上编程实现。Metropolis 算法是首个普适的采样方法，并启发了一系列 MCMC方法，所以人们把它视为随机模拟技术腾飞的起点。Metropolis的这篇论文被收录在《统计学中的重大突破》中，Metropolis算法也被遴选为二十世纪的十个最重要的算法之一。

我们接下来介绍的MCMC 算法是Metropolis算法的一个改进变种，即常用的 Metropolis-Hastings 算法。由上一节的例子和定理我们看到了，马氏链的收敛性质主要由转移矩阵P 决定, 所以基于马氏链做采样的关键问题是如何构造转移矩阵P,使得平稳分布恰好是我们要的分布p(x)。如何能做到这一点呢？我们主要使用如下的定理。

如果非周期马氏链的转移矩阵\(T(z^*,z)\)和分布\(\pi(x)\)满足$$\pi(i)T(i,j)=\pi(j)T(j,i)$$
则\(\pi(x)\)是马氏链的平稳分布，上式被称为细致平稳条件(detailed balance condition)。
满足细致平稳条件就能收敛到平稳分布，平稳分布的定义

\begin{split}
&\pi(z)=\sum_{z^} T(z^,z)\
&\pi=\pi P
\end{split}

\begin{split}
&\pi(z)=\sum_{z^} T(z^,z)\pi(z^*)\
&\pi=\pi P
\end{split}

我们将细致平稳条件公式带入到上式右边，即可推出平稳分布定义

\[\sum_{z^*} T(z^*,z)\pi(z^*)=\sum_{z^*} T(z,z^*)\pi(z)=\pi(z)\sum_{z^*} p(z^*|z)=\pi(z) \]

1.Metropolis-Hastings

假设我们已经有一个转移矩阵q(i,j)(表示从状态i转移到状态j的概率，也可以写为q(j|i)), 显然，通常情况下$$p(i)q(i,j)\neq p(j)q(j,i)$$
也就是细致平稳条件不成立，所以p(x)不太可能是这个马氏链的平稳分布。我们可否对马氏链做一个改造，使得细致平稳条件成立呢？譬如，我们引入一个a(i,j) , 我们希望$$p(i)\underbrace{ q(i,j)a(i,j) }{q'(i,j)}=p(j)\underbrace{ q(j,i)a(j,i) }{q'(j,i)}$$
取什么样的 以上等式能成立呢？最简单的，按照对称性，我们可以取$$a(i,j)=p(j)q(j,i)\qquad a(j,i)=p(i)q(i,j)$$
于是细致平稳条件就成立了，我们把原来具有转移矩阵Q的一个很普通的马氏链，改造为了具有转移矩阵Q'的马氏链，而Q'恰好满足细致平稳条件，由此马氏链的平稳分布就是p(x)!

在改造Q的过程中引入的a(i,j)称为接受率，物理意义可以理解为在原来的马氏链上，从状态i以q(i,j)的概率转跳转到状态j的时候，我们以a(i,j)的概率接受这个转移，于是得到新的马氏链Q'的转移概率为q(i,j)a(i,j)。

假设我们已经有一个转移矩阵Q(对应元素为q(i,j))，把以上的过程整理一下，我们就得到了如下的用于采样概率分布的算法

• basic Metropolis
• 假设：直接从分布p(z)中采样困难
• 问题：从分布p(z)中采样
• 算法：

1. 首先初始化马氏链的初始状态\(X_0=x_0\)
2. 第t时刻状态为\(x_t\)，对t+1时刻采样\(y\sim p(x|x_t)\)
3. 从均匀分布采样\(u\sim Uniform[0,1]\)
4. 如果\(u<a(x_t,y)=p(y)q(x_t|y)\)则接受转移\(x_t\sim y\)，即\(X_{t+1}=y\)
5. 否则不接受转移\(X_{t+1}=x_t\)
6. 跳转第2步

• 缺陷：马氏链在转移的过程中的接受率a(i,j)可能偏小，这样采样过程中马氏链容易原地踏步，拒绝大量的跳转，这使得马氏链遍历所有的状态空间要花费太长的时间，收敛到平稳分布p(x)的速度太慢。

有没有办法提升一些接受率呢?
假设a(i,j)=0.1 a(j,i)=0.2 , 此时满足细致平稳条件，于是$$p(i)q(i,j)\times 0.1=p(j)q(j,i)\times 0.2$$
上式两边扩大5倍，我们改写为$$p(i)q(i,j)\times 0.5=p(j)q(j,i)\times 1$$

看，我们提高了接受率，而细致平稳条件并没有打破！这启发我们可以把细致平稳条件式中的a(i,j)和a(j,i)同比例放大，使得两数中最大的一个放大到1，这样我们就提高了采样中的跳转接受率。所以我们可以取

\[a(i,j)=min\{\frac{p(j)q(j,i)}{p(i)q(i,j)},1\} \]

于是，经过对上述 basic Metropolis 采样算法中接受率的微小改造，我们就得到了如下教科书中最常见的 Metropolis-Hastings 算法。

• Metropolis-Hastings
• 假设：直接从分布p(z)中采样困难
• 问题：从分布p(z)中采样
• 算法：

1. 首先初始化马氏链的初始状态\(X_0=x_0\)
2. 第t时刻状态为\(x_t\)，对t+1时刻采样\(y\sim p(x|x_t)\)
3. 从均匀分布采样\(u\sim Uniform[0,1]\)
4. 如果\(u<a(x_t,y)=min\{\frac{p(y)q(x_t|y)}{p(x_t)q(y|x_t)},1\}\)则接受转移\(x_t\sim y\)，即\(X_{t+1}=y\)
5. 否则不接受转移\(X_{t+1}=x_t\)
6. 跳转第2步
   对于分布p(x) ,我们构造转移矩阵Q使其满足细致平稳条件。

2.Gibbs Sampling

对于高维的情形，由于接受率a的存在(通常a<1), 以上 Metropolis-Hastings 算法的效率不够高。能否找到一个转移矩阵Q使得接受率a=1呢？我们先看看二维的情形，假设有一个概率分布p(x,y) ，考察x坐标相同的两个点\(A(x_1,y_1)，A(x_1,y_2)\)我们发现

\begin{split}
& p(x_1,y_1)p(y_2|x_1)=p(x_1)p(y_1|x_1)p(y_2|x_1) \
& p(x_1,y_2)p(y_1|x_1)=p(x_1)p(y_2|x_1)p(y_1|x_1)
\end{split}

所以得到

\begin{split}
& p(x_1,y_1)p(y_2|x_1)=p(x_1,y_2)p(y_1|x_1) \qquad (1) \
& p(A)p(y_2|x_1)=p(B)p(y_1|x_1)
\end{split}

基于以上等式，我们发现，在\(x=x_1\)这条平行于y轴的直线上，如果使用条件分布\(p(y|x_1)\)做为任何两个点之间的转移概率，那么任何两个点之间的转移满足细致平稳条件。同样的，如果我们在\(y=y_1\)这条直线上任意取两个点 ，也有如下等式

\[p(A)p(x2|y1)=p(C)p(x1|y1) \]

于是我们可以如下构造平面上任意两点之间的转移概率矩阵Q

\begin{split}
& Q(A,B)=p(y_B|x_1)\qquad 如果x_A=x_B=x_1\
& Q(A,C)=p(x_C|y_1)\qquad 如果y_A=y_C=y_1\
& Q(A,D)=0\qquad
\end{split}

有了如上的转移矩阵Q, 我们很容易验证对平面上任意两点 , 满足细致平稳条件

\[p(X)Q(X,Y)=p(Y)Q(Y,X) \]

于是这个二维空间上的马氏链将收敛到平稳分布 。以上的过程我们很容易推广到高维的情形，对于(1)式，如果\(x_1\)变为多维情形\(\vec x_1\)，可以看出推导过程不变，所以细致平稳条件同样是成立的

\[p(\vec x_1,y_1)p(y_2|\vec x_1)=p(\vec x_1,y_2)p(y_1|\vec x_1) \]

此时转移矩阵Q由条件分布\(p(y|\vec x_1)\)定义。上式只是说明了一根坐标轴的情形，和二维情形类似，很容易验证对所有坐标轴都有类似的结论。所以n维空间中对于概率分布\(p(x_1,x_2,...,x_n)\)可以如下定义转移矩阵

\begin{split}
& Q(A,B)=p(x_i|x_1,...,x_{i-1},x_{i+1}...,x_n)\qquad 如果沿着x_i这根坐标轴做转移的时候\
& Q(A,D)=0\qquad others
\end{split}

• Gibbs Sampling
• 假设：直接从分布p(z)中采样困难
• 问题：从分布p(z)中采样
• 算法：

1. 首先初始化马氏链的初始状态\(X_0=x_0\)
2. 第t时刻状态为\(x_t\)，对t+1时刻采样
3. \(x_1^{(t+1)}\sim p(x_1|x_2^{(t)},x_3^{(t)},...,x_n^{(t)})\)
4. \(x_2^{(t+1)}\sim p(x_1|x_1^{(t+1)},x_3^{(t)},...,x_n^{(t)})\)
5. ...
6. \(x_j^{(t+1)}\sim p(x_1|x_1^{(t+1)},x_{j-1}^{(t+1)},...,x_{j+1}^{(t)},x_n^{(t)})\)
7. ...
8. \(x_n^{(t+1)}\sim p(x_1|x_1^{(t+1)},x_2^{(t)},...,x_{n-1}^{(t+1)})\)
9. 跳转第3步

以上算法收敛后，得到的就是概率分布\(p(x_1,x_2,..,x_n)\)的样本，当然这些样本并不独立，但是我们此处要求的是采样得到的样本符合给定的概率分布，并不要求独立。同样的，在以上算法中，坐标轴轮换采样不是必须的，可以在坐标轴轮换中引入随机性，这时候转移矩阵 中任何两个点的转移概率中就会包含坐标轴选择的概率，而在通常的Gibbs Sampling算法中，坐标轴轮换是一个确定性的过程，也就是在给定时刻，在一根固定的坐标轴上转移的概率是1。

应用案例

1.LDA

参考资料

1. http://www.52nlp.cn/lda-math-mcmc-和-gibbs-sampling2
2. PRML, chapter 11
3. http://blog.quantitations.com/inference/2012/11/24/rejection-sampling-proof/
4. http://thexbar.me/2014/11/07/reject-sample/
5. Notes on Pattern Recognition and Machine Learning (Jian Xiao)
6. Pattern Recognition And Machine Learning 读书会, chapter 11

Written with StackEdit.