首发于柴柴的科研分享

什么是马尔可夫链蒙特卡洛（MCMC）？

放放的闲聊社

清华大学电气工程及其自动化本科

许多人在学习贝叶斯统计、机器学习的过程中都听说过马尔可夫链蒙特卡洛（Markov Chain Monte Carlo ），但是其概念却很难理解，希望通过这边整理的文章给大家一个直观理解。

蒙特卡洛马尔可夫链是又由两个概念组成的，分别是蒙特卡洛（Monte Carlo）和马尔可夫链（Markov Chains），下面就将分别介绍一下两个概念。

什么是蒙特卡洛？

蒙特卡洛方法来自于摩纳哥的蒙特卡洛赌场，许多纸牌类游戏需要计算其胜利的概率。我们可以将蒙特卡洛理解为简单的模拟，通过模拟的情景来计算其发生的概率。

应用场景：通常情况下，许多概率的计算非常复杂甚至是无法计算的，但是我们可以通过计算机对期待发生的场景进行大量的模拟，从而计算出其发生的概率。简单的公式表达为：

P(发生的概率)=\frac{模拟情形下事件发生的次数}{模拟的次数}\\

学习过概率的同学应当知道这就是典型的频率学派的观点。

对于蒙特卡洛方法，经典的例子就是计算 \pi 值。在（-1，1）之间随机取两个数，如果在单位圆内，则记一次，在圆外则不计入次数。

P=\frac{点在圆内的次数}{总模拟次数}=\frac{\pi}{4}\\

import random
import numpy as np
from math import sqrt, pi

num1 = 1000 # simulation times 
freq1 = 0
x1 = []
y1 = []
for i in range(1,num1+1):
        x_i, y_i = random.uniform(-1, 1), random.uniform(-1, 1)
        square_distance = x_i**2+y_i**2
        x1.append(x_i)
        y1.append(y_i)
        if square_distance <= 1: # 圆的方程x^2+y^2=1
                freq1 += 1

simulation_pi1 = 4*freq1/num1;

num2 = 10000 # simulation times 
freq2 = 0
x2 = []
y2 = []
for i in range(num2):
        x_i, y_i = random.uniform(-1, 1), random.uniform(-1, 1)
        square_distance = x_i**2+y_i**2
        x2.append(x_i)
        y2.append(y_i)
        if square_distance <= 1: # 圆的方程x^2+y^2=1
                freq2 += 1

simulation_pi2 = 4*freq2/num2;

% matplotlib inline
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(8,8))

theta = np.linspace(0,2*pi,500)
x,y = np.cos(theta), np.sin(theta)
ax.plot(x, y, color='red', linewidth=1.0)

ax.scatter(x1,y1,color='blue',alpha=0.75)
ax.set_xlim(-1,1);
ax.set_ylim(-1,1);

结果为：

真实的pi值： 3.141592653589793
蒙特卡洛方法下得出的pi值（模拟1000次）： 3.188
相对误差（模拟1000次）： 1.4771917153924754 %
蒙特卡洛方法下得出的pi值（模拟10000次）： 3.1292
相对误差（模拟10000次）： 0.39447041536821975 %

小结：蒙特卡洛就是一种模拟事情发生的方法。

什么是马尔可夫链？

在学习马尔可夫链之前，我们首先了解一下什么是马尔可夫性质（或者叫马尔可夫特性，Markov Property）。

假设你有一个系统，这个系统具有 M 种可能的状态，并且在状态之间正在移动。真实世界中，我们常常见到这种例子，例如天气从炎热到温和，或者是股票市场从熊市到牛市再到停滞的状态（stagnant states）。

马尔可夫性质是指在给定的一个随机过程（stochastic process）中，在某一个时间点的状态为x_n ，其下一个状态 x_{n+1}=k 的概率（ k 代表 M 种状态里的任何一种），仅仅取决于现在这个它在给定时刻处于哪一种状态，而与它是如何达到现在的状态的无关。换句话说，该随机过程没有记忆性。

Wiki: In probability theory and statistics, the term Markov property refers to the memoryless property of a stochastic process.

用数学公式表达，即：

P(X_{n+1}=k|X_n=k_n,X_{n-1}=k_{n-1},...,X_1=k_1)=P(X_{n+1}=k|X_n=k_n)\\

如果一个过程具有马尔可夫性质，它也被成为马尔可夫过程（Markov Process）。

那么，为什么马尔可夫链这么重要呢？因为它的静态分布（或叫平稳性分布，stationary distribution）性质。

我们用一下Wiki上给出的股票市场的例子（见下图），这是一个具有静态空间的马尔可夫链（a continuous-time Markov chain with state-space），包含了牛市、熊市、停滞市场三种状态（Bull market, Bear market, Stagnant market）。

它各状态之间由转移概率组成的转移矩阵（ transition rate matrix）如下，记为 Q ：

以 Q_{11}=0.9 为例，其实就代表了上图中从bull market 转移到bull market的概率是0.9。

定义我们的状态为[bull, bear, stagnant]，假设我们从熊市（bear market）开始，则初始状态的矩阵 s_{0}=[0,1,0] 。那么我们可以利用转移矩阵Q来计算下一个状态的概率，即用初始状态 s_0*状态矩阵 Q ，表达为：

s_{1}=s_{0}Q=[0.15,0.8,0.05]\\

看得出来，下一个状态 s_1 整个概率之和仍然为1。对于第2个状态，计算公式应该为：

s_{2}=s_{1}Q=s_{0}Q^2\\

最终，你将会得到一个静止状态 sQ=s, s(Q-I)=0 ，求出的静态分布 s 为

s=[0.625,0.3125,0.0625]\\

利用Python计算一下，代码如下：

import numpy as np
Q = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]])
init_s = np.matrix([[0, 1 , 0]])
epsilon =1
while epsilon>10e-9:
    next_s = np.dot(init_s,Q)
    epsilon = np.sqrt(np.sum(np.square(next_s - init_s)))
    init_s = next_s
print(init_s)

结果为：

matrix([[0.62499998, 0.31250002, 0.0625    ]])

你可以尝试从其他状态出发，但最终的静态分布一定是相同的。那么得出这个静态分布有什么用处呢？

静止状态分布的重要性在于它可以让你在随机的时间里，对于一个系统的任意一个状态确定其概率。（The stationary state distribution is important because it lets you define the probability for every state of a system at a random time.）

对于上述例子，你可以说整个系统，62.5%的概率系统将处于牛市，31.25%的概率将处于熊市，而6.25%的概率将处于停滞态。

直观上说，你可以想想成在一个状态链上进行随机游走。你在一个状态出，然后你决定你前往下一个状态的方法是在给定当前状态下观察下一个状态的概率分布。

小结：马尔科夫链是一个数学对象，包含一系列状态以及状态之间的转移概率，如果每个状态转移到其他状态的概率只与当前状态相关，那么这个状态链就称为马尔科夫链。有这样一个马尔科夫链之后，我们可以任取一个初始点，然后根据状态转移概率进行随机游走。假设能够找到这样一个马尔科夫链，其状态转移概率正比于我们想要采样的分布（如贝叶斯分析中的后验分布），采样过程就变成了简单地在该状态链上移动的过程。

什么是蒙特卡洛马尔可夫链（MCMC）？

在了解什么是MCMC的时候，我们需要考虑一个情况。举例，我们已经知道一个分布（例如beta分布）的概率密度函数PDF，那么我们怎么样从这个分布中提取样本呢？

MCMC给我们提供了一种方式来从概率分布中进行采样，而这种方法在贝叶斯统计中的后验概率进行采样十分方便。

贝叶斯统计的基本公式（贝叶斯定理，Bayes Theorem）：

p(H|D)=\frac{p(H)p(D|H)}{p(D)}\\

1、 p(H|D) ：后验概率（posterior probability）。

2、 p(H) ：先验概率（prior probability）。先验代表着你在看见实际现象前相信事件发生的概率或概率分布，即在看到任何数据之前我们对 H 的理解。

3、 p(D|H) ：如果你的假设是正确的，看到现象发生时的可能性（也叫似然，likelihood）。换句话说，即我们认为的数据分布形式。

4、 p(D) ：正态化的常数，也叫做证据（evidence)。通常是为了满足概率分布或概率和为1的要求。（for example, the likelihood of that evidence under any circumstances）

我们可以对所有可能的参数求积分之后得到这个变量：

p(D)=\int_{H\in \mathcal{H}}p(H,D)dH\\

离散形式下，无非是：

p(D)=\sum_{H\in \mathcal{H}}p(H,D)\\

举例，如果H只取三个值，那么：

p(D)=p(H=H_1)p(D|H=H_1)+p(H=H_2)p(D|H=H_2)+p(H=H_3)p(D|H=H_3)\\

p(D) 的值将很容易进行计算。那如果H的值是一个连续值呢？那么 p(D) 将变成一个积分，很难直接求出。尽管这个公式看起来很简洁，但即使是对于一些很简单的模型，你也很难得到一个封闭形式的后验。

我们需要从后验分布中进行采样，但同时将 p(D) 视为一个常数。

给出MCMC的定义， Wikipedia的解释如下:

Markov Chain Monte Carlo (MCMC) methods are a class of algorithms for sampling from a probability distribution based on constructing a Markov chain that has the desired distribution as its stationary distribution. The state of the chain after a number of steps is then used as a sample of the desired distribution. The quality of the sample improves as a function of the number of steps.
马尔可夫链蒙特卡罗方法是一类以期望分布为平稳分布的马尔可夫链为基础，对概率分布进行抽样的算法。经过一系列步骤之后，链的状态将用作期望分布的样本。样品的质量随着步骤数量的增加而提高。

举例：假设我们想从一个Beta分布中提取样本，beta分布的概率密度函数为

f(x)=Cx^{\alpha-1}(1-x)^{\beta-1}, C=\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\\

不用被C的表达形式吓到，你可以简单理解为一个常数。事实上，C的存在是为了使概率密度函数的积分为1，其值是一个与 \alpha 、 \beta 有关的函数。那么根据这个概率密度函数，取不同的 \alpha 和 \beta 的值，分布也同样不同。

MCMC方法提供了可以创建一个以Beta分布为其平稳分布（stationary distribution）的马尔科夫链的算法，从而使取样变得简单，因为我们可以从一个均匀分布中取样(这相对容易)。

如果我们基于一些算法，重复地从一个随机的状态开始，遍历到下一个状态，我们将会创建一个马尔可夫链，一个以beta分布作为其平稳分布的马尔可夫链。

这类算法，一个典型的代表就是Metropolis-Hastings Algorithm算法。

Metropolis-Hastings Algorithm算法

为了更形象地理解这个算法，我们用下面这个例子来类比。假设我们想知道某个湖的水容量以及这个湖中最深的点，湖水很浑浊以至于没法通过肉眼来估计深度，而且这个湖相当大，网格近似法显然不是个好办法。为了找到一个采样策略，我们请来了两个好朋友小马和小萌。经过讨论之后想出了如下办法，我们需要一个船（当然，也可以是竹筏）和一个很长的棍子，这比声呐可便宜多了，而且我们已经把有限的钱都花在了船上。

如何决定接受还是拒绝新的测量值呢？这里的一个技巧便是使用Metropolis-Hastings准则，即接受新的测量值的概率正比于新旧两点的测量值之比。

事实上，理论保证了在这种情形下，如果我们能采样无数次，最终能得到完整的后验。幸运地是，实际上对于很多问题而言，我们只需要相对较少地采样就可以得到一个相当准确的近似。

数学表达：

让我们用数学语言来说明上述过程。如果 s=(s_1,s_2,….,s_M) 是我们想要的静态分布。我们想要创造一个马尔可夫链包含这个静态分布，那么这个马尔可夫链有 M 个状态，其转移矩阵 P 中的元素 p_{ij} 代表从状态 i 到状态 j 的概率。

那么我们要做的是：

（1）从一个随机的初始状态 i 开始

（2）随机选择一个目标状态 j （Proposal state），其选择概率与转移矩阵中第 i 行的概率有关

（3）计算一个测度（我们称之为接受概率，acceptance probability），计算方法为

a_{ij}=min(s_jp_{ji}/s_ip_{ij},1)\\

（4）从[0,1]中随机生成一个数，如果数字大于0.5，接受这个转移；如果小于0.5，拒绝这次的状态转移。

（5）重复上述过程

经过一段时间后，这个马尔可夫链就会收敛到我们想要的静态分布 s ，我们就可以用这个链上的状态作为这个分布的样本。

如果以beta分布为例，那么我们唯一用到PDF的时候就是求接受概率的时刻，而且由于是除法，常数C就在这个过程中消掉了。

假设在[0,1]上有一个具有无限状态的任意马尔可夫链，并且有一个过渡矩阵 P ，使得 p_{ij}=p_{ji}= 矩阵中的所有元素。应用上述方法，我们可以得到的是：

a_{ij}=min(s_jp_{ji}/s_ip_{ij},1)=min(s_j/s_i,1)\\ s_{i}=Ci^{\alpha-1}(1-i)^{\beta-1}\\ s_{j}=Cj^{\alpha-1}(1-j)^{\beta-1}\\

代码实现MCMC采样beta分布：

import random
# Lets define our Beta Function to generate s for any particular state. 
# We don't care for the normalizing constant here.

def beta_s(w,a,b): # beta distribution pdf
    return w**(a-1)*(1-w)**(b-1)

# This Function returns True if the coin with probability P of heads comes heads when flipped.
def random_coin(p):
    unif = random.uniform(0,1)
    if unif>=p:
        return False
    else:
        return True

# This Function runs the MCMC chain for Beta Distribution.
def beta_mcmc(N_hops,a,b):
    states = []
    cur = random.uniform(0,1)
    for i in range(0,N_hops):
        states.append(cur)
        next = random.uniform(0,1)
        ap = min(beta_s(next,a,b)/beta_s(cur,a,b),1) # Calculate the acceptance probability
        if random_coin(ap):
            cur = next
    return states[-1000:] # Returns the last 1000 states of the chain

import numpy as np
from matplotlib import pylab as pl
import scipy.special as ss
% matplotlib inline
plt.rcParams['figure.figsize'] = (17.0, 4.0)

# Actual Beta PDF.
def beta(a, b, i):
    e1 = ss.gamma(a + b)
    e2 = ss.gamma(a)
    e3 = ss.gamma(b)
    e4 = i ** (a - 1)
    e5 = (1 - i) ** (b - 1)
    return (e1/(e2*e3)) * e4 * e5

# Create a function to plot Actual Beta PDF with the Beta Sampled from MCMC Chain.
def plot_beta(a, b):
    Ly = []
    Lx = []
    i_list = np.mgrid[0:1:100j]
    for i in i_list:
        Lx.append(i)
        Ly.append(beta(a, b, i))
    pl.plot(Lx, Ly, label="Real Distribution: a="+str(a)+", b="+str(b))
    plt.hist(beta_mcmc(100000,a,b),density=True,bins =25, 
            histtype='step',label="Simulated_MCMC: a="+str(a)+", b="+str(b))
    pl.legend()
    pl.show()
    
plot_beta(0.1, 0.1)
plot_beta(1, 1)
plot_beta(2, 3)

Metropolis-Hastings Algorithm的进一步说明

从贝叶斯的角度看，Metropolis-Hastings算法使得我们能够从任意分布中以概率p(x)得到采样值，只要我们能算出某个与p(x)成比例的值。这一点很有用，因为在类似贝叶斯统计的许多问题中，最难的部分是计算归一化因子，也就是贝叶斯定理中的分母。Metropolis-Hastings算法的步骤如下：

（1）给参数 x_i 赋一个初始值，通常是随机初始化或者使用某些经验值。

（2）从某个简单的分布 Q(x_{i+1}|x_i) 中选一个新的值 x_{i+1} ，如高斯分布或者均匀分布。这一步可以看做是对状态的扰动。

（3）根据Metropolis-Hastings准则计算接受一个新的参数值的概率：

p_a(x_{i+1}|x_i)=min(1,\frac{p(x_{i+1})}{p(x_{i})}\frac{q(x_{i}|x_{i+1})}{q(x_{i+1}|x_i)})\\

（4）从位于区间[0,1]内的均匀分布中随机选一个值，如果第（3）步中得到的概率值比该值大，那么就接受新的值，否则仍保持原来的值。

（5）回到第（2）步重新迭代，直到我们有足够多的样本。

如果选取的分布 Q(x_{i+1}|x_i) 是对称的，那么可以得到 p_a(x_{i+1}|x_i)=min(1,\frac{p(x_{i+1})}{p(x_{i})}) ，这通常称为Metropolis准则。

步骤（3）和步骤（4）表明：我们总是会转移到一个比当前状态（或参数）概率更大的状态（或参数），对于概率更小的，则会以 x_{i+1} 与 x_{i} 之比的概率接受。该准则中的接受步骤使得采样过程相比网格近似方法更高效，同时保证了采样的准确性。

目标分布（贝叶斯统计中的后验分布）是通过记录下来的采样值来近似的。如果我们接受转移到新的状态 x_{i+1} ，那么我们就记录该采样值 x_{i+1} 。如果我们拒绝转移到 x_{i+1} ，那么我们就记录 x_i 。

最后，我们会得到一连串记录值，有时候也称采样链或者迹。如果一切都正常进行，那么这些采样值就是后验的近似。在采样链中出现次数最多的值就是对应后验中最可能的值。该过程的一个优点是：后验分析很简单，我们把对后验求积分的过程转化成了对采样链所构成的向量求和的过程。

强烈建议阅读以下博文，作者用一个简单的例子实现了metropolis方法，并将其用于求解后验分布，文中用非常好看的图展示了采样的过程，同时简单讨论了最初选取的步长是如何影响结果的。

一句话总结MCMC

MCMC generates samples from the posterior distribution by constructing a reversible Markov-chain that has as its equilibrium distribution the target posterior distribution.
https://twiecki.io/blog/2015/11/10/mcmc-sampling/

参考文献

1、 Rahul Agarwal，“MCMC Intuition for Everyone“， https://towardsdatascience.com/mcmc-intuition-for-everyone-5ae79fff22b1

2、“MCMC sampling for dummies”， https://twiecki.io/blog/2015/11/10/mcmc-sampling/

2、Wiki，“Markov Chain”， https://en.wikipedia.org/wiki/Markov_chain#Steady-state_analysis_and_limiting_distributions

3、Osvaldo Martin，“Bayesian Analysis with Python”（用贝叶斯做数据分析），代码仓库： https://github.com/PacktPublishing/Bayesian-Analysis-with-Python-Second-Edition

4、Wiki，“Markov chain Monte Carlo”， https://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo

编辑于 2020-04-15 17:44

马尔可夫链

马尔科夫过程

贝叶斯统计