数学推导+纯Python实现马尔可夫链蒙特卡洛

最新推荐文章于 2024-10-08 09:18:30 发布

文文学霸

最新推荐文章于 2024-10-08 09:18:30 发布

阅读量2k

点赞数 4

Python机器学习算法实现

Author：louwill

Machine Learning Lab

蒙特卡洛（Monte Carlo，MC）方法作为一种统计模拟和近似计算方法，是一种通过对概率模型随机抽样进行近似数值计算的方法。马尔可夫链（Markov Chain，MC）则是一种具备马尔可夫性的随机序列。将二者结合起来便有了马尔可夫链蒙特卡洛方法（Markov Chain Monte Carlo，MCMC），即是以构造马尔可夫链为概率模型的蒙特卡洛方法。

限于篇幅，本文不对MCMC的前置知识进行详细介绍，关于蒙特卡洛方法和马尔可夫链的大量内容，请各位读者自行查阅相关材料进行学习。本文聚焦于MCMC方法本身原理和常用实现方法。

MCMC简介

一般来说，对目标概率模型进行随机抽样能够帮助我们得到该分布的近似数值解。但如果随机变量的多元的，或者所要抽样的概率密度函数形式是复杂的非标准式时，直接应用蒙特卡洛方法就会很困难。

MCMC方法的基本思路是：在随机变量的状态空间上定义一个马尔可夫链，使其平稳分布就是我们要抽样的目标分布。然后基于该马尔可夫链进行随机游走产生对应的样本序列进行数值计算。当时间足够长时，抽样所得到的分布就会趋近于平稳分布，基于该马尔可夫链的抽样结果就是目标概率分布的抽样结果。对抽样结果的函数均值计算就是目标数学期望值。

我们将上述流程梳理一下，完整的MCMC方法包括以下3个步骤：

在随机变量的状态空间上定义一个满足遍历定理的马尔可夫链，使其平稳分布为目标分布；
从状态空间种某一点出发，在构造的马尔可夫链上进行随机游走，可以得到样本序列；
根据马尔可夫链的遍历定理，确定正整数和，可得样本集合，最后可得函数的遍历均值：

所以MCMC的关键问题在于如何构造满足条件的马尔可夫链。常用的MCMC构建算法包括Metropolis-Hasting算法和Gibbs抽样。

Metropolis-Hasting算法

Metropolis-Hasting算法也可以简称为M-H采样，该算法由Metropolis提出后经Hasting改进，故而得名。假设目标抽样分布为，M-H算法采用的状态转移核为，则其马尔可夫链可定义为：

其中和分别为建议分布和接受分布。建议分布可以是另一个马尔可夫链的转移核，其形式也可能有多种。包括对称形式和独立抽样形式。假设建议分布是对称的，即对任意的和，有：

接受分布形式如下：

转移核为的马尔可夫链的随机游走过程如下：如果在时刻处于状态，即有，则先按照建议分布抽样产生一个候选状态，然后按照接受分布抽样决定是否接受状态。以概率接受，以概率拒绝。完整的Metropolis-Hasting算法流程如下：

在状态空间上任意选择一个初始值；
对遍历执行
- 设状态，按照建议分布抽取一个候选状态。
- 计算接受概率
- 从区间上按均匀分布随机抽取一个数，若，则状态；否则。
最后得到样本集合，计算：

借助Python的高级计算库scipy，下面给出Metropolis-Hasting算法的基本实现过程。假设目标平稳分布是一个正态分布，基于M-H的采样过程如下。

from scipy.stats import norm
import random
import matplotlib.pyplot as plt


# 定义平稳分布
def smooth_dist(theta):
    y = norm.pdf(theta, loc=3, scale=2)
    return y


T = 10000
pi = [0 for i in range(T)]
sigma = 1
# 设置初始值
t = 0
# 遍历执行
while t < T-1:
    t = t + 1
    # 状态转移进行随机抽样
    pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1, random_state=None)  
    # 计算接受概率
    alpha = min(1, (smooth_dist(pi_star[0]) / smooth_dist(pi[t - 1])))   
    # 从均匀分布中随机抽取一个数u
    u = random.uniform(0, 1)
    # 拒绝-接受采样
    if u < alpha:
        pi[t] = pi_star[0]
    else:
        pi[t] = pi[t - 1]


# 绘制采样分布
plt.scatter(pi, norm.pdf(pi, loc=3, scale=2), label='Target Distribution')
num_bins = 100
plt.hist(pi, 
         num_bins, 
         normed=1, 
         facecolor='red', 
         alpha=0.6, 
         label='Samples Distribution')
plt.legend()
plt.show();

经过10000次遍历迭代后，采样效果如下图所示。

Gibbs抽样

相较于M-H抽样，Gibbs抽样是目前更常用的MCMC抽样算法，Gibbs抽样可以视为特殊的M-H抽样方法。Gibbs抽样适用于多元随机变量联合分布的抽样和估计，其基本思路是从联合概率分布种定义满条件概率分布，依次对满条件概率分布进行抽样得到目标样本序列。

这里简单提一下满条件分布。假设MCMC的目标分布为多元联合概率分布，如果条件概率分布中所有个变量全部出现，其中，，这种条件概率分布即为满条件分布。

假设在第步得到样本，在第步先对第一个随机变量按照如下满条件分布进行随机抽样得到 :
之后依次对第个变量按照满条件概率分布进行抽样，最后可得整体样本。

Gibbs抽样可视为单分量M-H抽样的特殊情形。即Gibbs抽样对每次抽样结果都以1的概率进行接受，从不拒绝，这跟M-H采样有本质区别。Gibbs抽样的完整流程如下：

给定多随机变量初始值。
对遍历执行：假设第步得到样本，则第次迭代进行如下步骤：
- 由满条件分布抽取得到
- 由满条件分布抽取
- 由满条件分布抽取
最后得到样本集合，计算：

Gibbs抽样与单分量的M-H算法不同之处在于Gibbs抽样不会在一些样本上做停留，即抽样不会被拒绝。Gibbs抽样适用于满条件分布容易抽样的情况，而单分量的M-H算法适用于满条件分布不易抽样的情况，此时可选择容易抽样的条件分布作为建议分布。

下面以二维正态分布为例给出多元随机变量的Gibbs采样Python实现过程。

import math
# 导入多元正态分布函数
from scipy.stats import multivariate_normal
# 指定二元正态分布均值和协方差矩阵
samplesource = multivariate_normal(mean=[5,-1], cov=[[1,0.5],[0.5,2]])


# 定义给定x的条件下y的条件状态转移分布
def p_yx(x, m1, m2, s1, s2):
    return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt(1 - rho ** 2) * s2))
# 定义给定y的条件下x的条件状态转移分布
def p_xy(y, m1, m2, s1, s2):
    return (random.normalvariate(m1 + rho * s1 / s2 * (y - m2), math.sqrt(1 - rho ** 2) * s1))


# 指定相关参数
N, K= 5000, 20
x_res = []
y_res = []
z_res = []
m1, m2 = 5, -1
s1, s2 = 1, 2
rho, y = 0.5, m2


# 遍历迭代
for i in range(N):
    for j in range(K):
        # y给定得到x的采样
        x = p_xy(y, m1, m2, s1, s2)  
        # x给定得到y的采样
        y = p_yx(x, m1, m2, s1, s2)   
        z = samplesource.pdf([x,y])
        x_res.append(x)
        y_res.append(y)
        z_res.append(z)
# 绘图
num_bins = 50
plt.hist(x_res, num_bins, normed=1, facecolor='green', alpha=0.5,label='x')
plt.hist(y_res, num_bins, normed=1, facecolor='red', alpha=0.5,label='y')
plt.title('Histogram')
plt.legend()
plt.show();

采样效果如下：

二维效果展示：

from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = Axes3D(fig, rect=[0, 0, 1, 1], elev=30, azim=20)
ax.scatter(x_res, y_res, z_res, marker='o')
plt.show();