从随机过程到马尔科夫链蒙特卡洛方法（MCMC）

从随机过程到马尔科夫链蒙特卡洛方法

1. Introduction

第一次接触到 Markov Chain Monte Carlo (MCMC) 是在 theano 的 deep learning tutorial 里面讲解到的 RBM 用到了 Gibbs sampling，当时因为要赶着做项目，虽然一头雾水，但是也没没有时间仔细看。趁目前比较清闲，把 machine learning 里面的 sampling methods 理一理，发现内容还真不少，有些知识本人也是一知半解，所以这篇博客不可能面面俱到详细讲解所有的 sampling methods，而是着重讲一下这个号称二十世纪 top 10 之一的算法—— Markov chain Monte Carlo。在介绍 MCMC 之前，我们首先了解一下 MCMC 的 Motivation 和在它之前用到的方法。本人也是初学者，错误在所难免，欢迎一起交流。

这篇博客从零开始，应该都可以看懂，主要内容包括：

随机采样
拒绝采样
重要性采样
Metropolis-Hastings Algorithm
Gibbs Sampling

2. Sampling

我们知道，计算机本身是无法产生真正的随机数的，但是可以根据一定的算法产生伪随机数（pseudo-random numbers)。最古老最简单的莫过于 Linear congruential generator：

式子中的 a 和 c 是一些数学知识推导出的合适的常数。但是我们看到，这种算法产生的下一个随机数完全依赖现在的随机数的大小，而且当你的随机数序列足够大的时候，随机数将出现重复子序列的情况。当然，理论发展到今天，有很多更加先进的随机数产生算法出现，比如 python 数值运算库 numpy 用的是 Mersenne Twister 等。但是不管算法如何发展，这些都不是本质上的随机数，用冯诺依曼的一句话说就是：

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

要检查一个序列是否是真正的随机序列，可以计算这个序列的 entropy 或者用压缩算法计算该序列的冗余。

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

假如在 y 轴上产生（0,1）之间的均匀分布的随机数，水平向右投影到高斯累计分布函数上，然后垂直向下投影到 x 轴，得到的就是高斯分布。可见高斯分布的随机数实际就是均匀分布随机数在高斯分布的 CDF 函数下的逆映射。当然，在实际操作中，更有效的计算方法有 Box–Muller_transform (an efficient polar form)，Ziggurat algorithm 等，这些方法 tricky and faster，没有深入了解，这里也不多说了。

3. Motivation

MCMC 可解决高维空间里的积分和优化问题：

上面一个例子简单讲了利用高斯分布的 CDF 可以产生高斯随机数，但是有时候我们遇到一些分布的 CDF 计算不出来（无法用公式表示），随机数如何产生？
遇到某些无法直接求积分的函数，如 e^{x^2}，在计算机里面如何求积分？
如何对一个分布进行高效快速的模拟，以便于抽样？
如何在可行域很大(or large number of possible configurations)时有效找到最优解——RBM 优化目标函数中的问题。
如何在众多模型中快速找到更好的模型——MDL, BIC, AIC 模型选择问题。

3.1 The Monte Carlo principle

实际上，Monte Carlo 抽样基于这样的思想：假设玩一局牌的赢的概率只取决于你抽到的牌，如果用穷举的方法则有 52! 种情况，计算复杂度太大。而现实中的做法是先玩几局试试，统计赢的概率，如果你不太确信这个概率，你可以尽可能多玩几局，当你玩的次数很大的时候，得到的概率就非常接近真实概率了。

上述方法可以估算随机事件的概率，而用 Monte Carlo 抽样计算随即变量的期望值是接下来内容的重点：X 表示随即变量，服从概率分布 p(x), 那么要计算 f(x) 的期望，只需要我们不停从 p(x) 中抽样

当抽样次数足够的时候，就非常接近真实值了：

Monte Carlo 抽样的方法还有一个重要的好处是：估计值的精度与 x 的维度无关（虽然维度越高，但是每次抽样获得的信息也越多），而是与抽样次数有关。在实际问题里面抽样二十次左右就能达到比较好的精度。

但是，当我们实际动手的时候，就会发现一个问题——如何从分布 p(x) 中抽取随机样本。之前我们说过，计算可以产生均匀分布的伪随机数。显然，第二小节产生高斯随机数的抽样方法只对少数特定的问题管用，对于一般情况呢？

3.2 Rejection Sampling

Reject Sampling 实际采用的是一种迂回( proposal distribution q(x) )的策略。既然 p(x) 太复杂在程序中没法直接采样，那么我设定一个程序可抽样的分布 q(x) 比如高斯分布，然后按照一定的方法拒绝某些样本，达到接近 p(x) 分布的目的：

具体操作如下，设定一个方便抽样的函数 q(x)，以及一个常量 k，使得 p(x) 总在 kq(x) 的下方。（参考上图）

x 轴方向：从 q(x) 分布抽样得到 a。(如果是高斯，就用之前说过的 tricky and faster 的算法更快）
y 轴方向：从均匀分布（0, kq(a)) 中抽样得到 u。
如果刚好落到灰色区域： u > p(a), 拒绝，否则接受这次抽样
重复以上过程

用均匀分布拒绝抽样来近似两个高斯混合分布的代码如下：

rejectionsampling.py

# -*- coding=utf8 -*-

# Code from Chapter 14 of Machine Learning: An Algorithmic Perspective
# The basic rejection sampling algorithm

from pylab import *
from numpy import *

def qsample():
    return random.rand()*4.

def p(x):
    return 0.3*exp(-(x-0.3)**2) + 0.7* exp(-(x-2.)**2/0.3) 

def rejection(nsamples):
    
    M = 0.72#0.8
    samples = zeros(nsamples,dtype=float)
    count = 0
    for i in range(nsamples):
        accept = False
        while not accept:
            x = qsample()
            u = random.rand()*M
            if u<p(x):
                accept = True
                samples[i] = x
            else: 
                count += 1
    print count   
    return samples

x = arange(0,4,0.01)
x2 = arange(-0.5,4.5,0.1)
realdata = 0.3*exp(-(x-0.3)**2) + 0.7* exp(-(x-2.)**2/0.3) 
box = ones(len(x2))*0.75#0.8
box[:5] = 0
box[-5:] = 0
plot(x,realdata,'k',lw=6)
plot(x2,box,'k--',lw=6)

import time
t0=time.time()
samples = rejection(10000)
t1=time.time()
print "Time ",t1-t0

hist(samples,15,normed=1,fc='k')
xlabel('x',fontsize=24)
ylabel('p(x)',fontsize=24)
axis([-0.5,4.5,0,1])
show()

在高维的情况下，Rejection Sampling 会出现两个问题，第一是合适的 q 分布比较难以找到，第二是很难确定一个合理的 k 值。这两个问题会导致拒绝率很高，无用计算增加。

3.3 Importance Sampling

Importance Sampling 也是借助了容易抽样的分布 q (proposal distribution)来解决这个问题，直接从公式出发：

其中，p(z) / q(z) 可以看做 importance weight。我们来考察一下上面的式子，p 和 f 是确定的，我们要确定的是 q。要确定一个什么样的分布才会让采样的效果比较好呢？直观的感觉是，样本的方差越小期望收敛速率越快。比如一次采样是 0, 一次采样是 1000, 平均值是 500,这样采样效果很差，如果一次采样是 499, 一次采样是 501, 你说期望是 500,可信度还比较高。在上式中，我们目标是 p×f/q 方差越小越好，所以 |p×f| 大的地方，proposal distribution q(z) 也应该大。举个稍微极端的例子：

第一个图表示 p 分布，第二个图的阴影区域 f = 1，非阴影区域 f = 0, 那么一个良好的 q 分布应该在左边箭头所指的区域有很高的分布概率，因为在其他区域的采样计算实际上都是无效的。这表明 Importance Sampling 有可能比用原来的 p 分布抽样更加有效。

但是可惜的是，在高维空间里找到一个这样合适的 q 非常难。即使有 Adaptive importance sampling 和 Sampling-Importance-Resampling(SIR) 的出现，要找到一个同时满足 easy to sample 并且 good approximations 的 proposal distribution, it is often impossible！

4. MCMC Algorithm

上面说了这么多采样方法，其实最终要突出的就是 MCMC 的过人之处。MCMC 的绝妙之处在于：通过稳态的 Markov Chain 进行转移计算，等效于从 P(x) 分布采样。但是在了解 MCMC 具体算法之前，我们还要熟悉 Markov Chain 是怎么一回事。

4.1 Markov Chain

Markov Chain 体现的是状态空间的转换关系，下一个状态只决定与当前的状态(可以联想网页爬虫原理，根据当前页面的超链接访问下一个网页)。如下图：

这个状态图的转换关系可以用一个转换矩阵 T 来表示：

举一个例子，如果当前状态为 u(x) = (0.5, 0.2, 0.3), 那么下一个矩阵的状态就是 u(x)T = (0.18, 0.64, 0.18), 依照这个转换矩阵一直转换下去，最后的系统就趋近于一个稳定状态 (0.22, 0.41, 0.37) (此处只保留了两位有效数字)。而事实证明无论你从那个点出发，经过很长的 Markov Chain 之后都会汇集到这一点。

考虑一般的情况，满足什么条件下经过很长的 Markov Chain 迭代后系统分布会趋近一个稳定分布，即最后的 u(x) 等效于从目标分布 p(x) 采样。大概的条件如下（自己随便总结的，可能有遗漏和错误）：

Irreducibility. 即图是联通的，T 矩阵不能被切豆腐一样划分成小方块，举个例子，比如爬虫爬不到内部局域网的网页
Aperiodicity. 即图中遍历不会陷入到一个死圈里，有些网站为了防机器人，会专门设置这种陷阱
Detailed Balance，这是保证系统有稳态的一个重要条件，详细说明见下面。

假设 p(x) 是最后的稳态，那么 detailed balance 可以用公式表示为：

什么意思呢？假设上面状态图 x1 有 0.22 元， x2 有 0.41 元，x3 有 0.37 元，那么 0.22×1 表示 x1 需要给 x2 钱，以此类推，手动计算，可以发现下一个状态每个人手中的钱都没有变。值得说明的是，这里体现了一个很重要的特性，那就是从一个高概率状态 xi 向一个低概率状态 x(i-1) 转移的概率等于从这个低概率状态向高概率状态转移的概率（reversible，至于要不要转移又是另外一回事)。当然，在上面一个例子中，情况比较特殊，等号两边其实都是同一个东西。马氏链的收敛性质主要由转移矩阵决定, 所以基于马氏链做采样的关键问题是如何构造转移矩阵,使得平稳分布恰好是我们要的分布p(x)。但是考虑一维的情况，假设 p(x) 是一维高斯分布，x 是根据 markov chain 得到的一个样本，依照上面的等式，那么我们可以根据转移矩阵 T左和 T右（这里实际是 proposal distribution）来得到 p(xi) 和 p(x(i-1)) 的比率，进而按照一定的概率对这两个样本进行选择。通过大量这样的处理，得到样本就符合原始的 p(x) 分布了。这就是 MH 算法的基本原理。

4.2 Metropolis-Hastings Algorithm

Metropolis Hasting（下面简称MH）是蒙特卡罗马尔科夫链中一种重要的抽样方法。本文简要介绍MH算法，并给出一个实例。

MH算法在参数空间随机取值，作为起始点。按照参数的概率分布生成随机的参数，按照这一系列参数的组合，计算当前点的概率密度。依据当前点和起始点概率密度比值是否大于(0,1)之间的随机数来判断是否保留当前点。若当前点的概率密度大于该随机数，就称这种状态为接受状态，此时，在满足参数概率分布的前提下，继续随机抽取参数的组合，作为下一点，计算下一点的概率密度，并计算下一点概率密度和概率密度的比值，并继续循环。若当前点不能被接受，则继续在满足参数概率分布的前提下，继续生成随机数，作为新的参数组合，直到参数组合能够被接受为止。

举个例子，我们要用 MH 算法对标准高斯分布进行采样，转移函数(对称)是方差为 0.05 的高斯矩阵，上述算法过程如下：

选取一个随机点 x0，作为一个采样点
以 x0 为中心，以转移函数为分布采取随机点 x1
以算法中的 A 概率接受 x1, 否则接受 x0
重复第二步第三步

注意到高斯分布是一个径向基函数，上面算法画波浪线的部分相等。

matlab 代码如下：

n = 250000;
x = zeros(n, 1);
x(1) = 0.5;
for i = 1: n-1
    x_c = normrnd(x(1), 0.05);
    if rand < min(1, normpdf(x_c)/normpdf(x(i)))
        x(i+1) = x_c;
    else
        x(i+1) = x(i);
    end
end

MH 算法中的 proposal distribution q(x) 也是需要小心确定的，详细知识可以查阅这篇介绍论文 (An introduction to MCMC for machine learning, Andrieu, Christophe). 可以看到，这个算法和模拟退火算法的思想是非常相似的，但是在模拟退火算法过程中，随着时间的增加，接受值大的区域的概率越来越高，直到找到最高点。

4.3 Gibbs Sampling

Gibbs Sampling 实际上是 MH 算法的一个变种。具体思路如下：假设在一定温度下一定量的分子在容器里做无规则的热运动，如何统计系统的能量呢？同样，我们用 Monte Carlo 的思想进行统计计算。我们假设所有的分子静止在某一个时刻，这是初识状态。固定其他的分子，根据分子间的作用力对其中一个分子进行移动，也就是说在该分子以一定的概率移动到领域的某一个地方，移动完了之后再静止。然后基于移动后的状态对下一个分子进行同样的移动操作...直到所有的分子移动完毕，那么现在的状态就是 Monte Carlo 采样的第二个样本。依照这样的顺序采样下去，我们对于这个系统就能计算一个统计意义上的能量了。从条件分布的角度来看，算法过程如下：

总体来讲，Gibbs Sampling 就是从条件概率中选择一个变量（分子），然后对该变量（分子）进行采样。当所有变量采样完毕之后，就得到了后面的一个状态，从而完成了对系统配置的采样。在 deep learning 的 RBM 中，gibbs 采样是已知权重参数和一个 v 变量，通过采样得到 h。通过 h 采样又可以得到另一个 v ，如此交替采样，就可以逐渐收敛于联合分布了。

Gibbs.py (对高斯分布进行 Gibbs 采样)

# -*- coding=utf8 -*-

# Code from Chapter 14 of Machine Learning: An Algorithmic Perspective
# A simple Gibbs sampler

from pylab import *
from numpy import *

def pXgivenY(y,m1,m2,s1,s2):
    return random.normal(m1 + (y-m2)/s2,s1)

def pYgivenX(x,m1,m2,s1,s2):
    return random.normal(m2 + (x-m1)/s1,s2)

def gibbs(N=5000):
    k=20
    x0 = zeros(N,dtype=float)
    m1 = 10
    m2 = 20
    s1 = 2
    s2 = 3
    for i in range(N):
        y = random.rand(1)
        # 每次采样需要迭代 k 次
        for j in range(k):
            x = pXgivenY(y,m1,m2,s1,s2)
            y = pYgivenX(x,m1,m2,s1,s2)
        x0[i] = x
    
    return x0

def f(x):
    return exp(-(x-10)**2/10)

# 画图
N=10000
s=gibbs(N)
x1 = arange(0,17,1)
hist(s,bins=x1,fc='k')
x1 = arange(0,17, 0.1)
px1 = zeros(len(x1))
for i in range(len(x1)):
    px1[i] = f(x1[i])
plot(x1, px1*N*10/sum(px1), color='k',linewidth=3)

show()

ps:

modified in 2013.11.1: 偶然发现统计之都有一篇类似的博客，gibbs采样解释得更加详细更加恰当，^_^，请点击这里

参考文献：

[1] PRML, chapter 11

[2] An introduction to MCMC for machine learning, Andrieu, Christophe

[3] 随机模拟的基本思想和常用采样方法（sampling）

[4] youtube 上的讲解 MCMC 的视频

作者：daniel-D 出处：http://www.cnblogs.com/daniel-D/ 欢迎转载或分享，但请务必声明文章出处。

posted @ 2016-12-08 13:49 stardsd 阅读(16859) 评论(2) 编辑收藏举报

刷新页面返回顶部

赏月斋

慎终如始宁静致远

从随机过程到马尔科夫链蒙特卡洛方法（MCMC）