复现经典：《统计学习方法》第19章马尔可夫链蒙特卡罗法

原创

mb5fe18e9fef50b 2020-12-24 20:40:30 ©著作权

文章标签 Java 文章分类 Java 后端开发

©著作权归作者所有：来自51CTO博客作者mb5fe18e9fef50b的原创作品，请联系作者获取转载授权，否则将追究法律责任

第19章马尔可夫链蒙特卡罗法

本文是李航老师的《统计学习方法》一书的代码复现。作者：黄海广

备注：代码都可以在github中下载。我将陆续将代码发布在公众号“机器学习初学者”，可以在这个专辑在线阅读。

蒙特卡罗法是通过基于概率模型的抽样进行数值近似计算的方法，蒙特卡罗法可以用于概率分布的抽样、概率分布数学期望的估计、定积分的近似计算。

随机抽样是蒙特卡罗法的一种应用，有直接抽样法、接受拒绝抽样法等。接受拒绝法的基本想法是，找一个容易抽样的建议分布，其密度函数的数倍大于等于想要抽样的概率分布的密度函数。按照建议分布随机抽样得到样本，再按要抽样的概率分布与建议分布的倍数的比例随机决定接受或拒绝该样本，循环执行以上过程。

可逆马尔可夫链是满足遍历定理的充分条件。

3 . 马尔可夫链蒙特卡罗法是以马尔可夫链为概率模型的蒙特卡罗积分方法，其基本想法如下：

蒙特卡洛法（Monte Carlo method) ，也称为统计模拟方法 (statistical simulation method) ，是通过从概率模型的随机抽样进行近似数值计

算的方法。马尔可夫链陟特卡罗法 (Markov Chain Monte Carlo, MCMC)，则是以马尔可夫链 (Markov chain）为概率模型的蒙特卡洛法。

马尔可夫链蒙特卡罗法构建一个马尔可夫链，使其平稳分布就是要进行抽样的分布，首先基于该马尔可夫链进行随机游走，产生样本的序列，

之后使用该平稳分布的样本进行近似数值计算。

Metropolis-Hastings算法是最基本的马尔可夫链蒙特卡罗法，Metropolis等人在 1953年提出原始的算法，Hastings在1970年对之加以推广，

形成了现在的形式。吉布斯抽样(Gibbs sampling)是更简单、使用更广泛的马尔可夫链蒙特卡罗法，1984 年由S. Geman和D. Geman提出。

马尔可夫链蒙特卡罗法被应用于概率分布的估计、定积分的近似计算、最优化问题的近似求解等问题，特别是被应用于统计学习中概率模型的学习

与推理，是重要的统计学习计算方法。

一般的蒙特卡罗法有直接抽样法、接受-拒绝抽样法、重要性抽样法等。

接受-拒绝抽样法、重要性抽样法适合于概率密度函数复杂（如密度函数含有多个变量，各变量相互不独立，密度函数形式复杂），不能直接抽样的情况。

19.1.2 数学期望估计

一舣的蒙特卡罗法，如直接抽样法、接受·拒绝抽样法、重要性抽样法，也可以用于数学期望估计（estimation Of mathematical expectation)。

吉布斯采样

网络资源：

LDA-math-MCMC 和 Gibbs Sampling： https://cosx.org/2013/01/lda-math-mcmc-and-gibbs-sampling

MCMC蒙特卡罗方法： https://www.cnblogs.com/pinard/p/6625739.html

1. import random
2. import math
3. import matplotlib.pyplot as plt
4. import seaborn as sns
5. import numpy as np
6. 
7. transfer_matrix = np.array([[0.6, 0.2, 0.2], [0.3, 0.4, 0.3], [0, 0.3, 8.7]],
9.                            dtype='float32')
10. start_matrix = np.array([[0.5, 0.3, 0.2]], dtype='float32')
11. 
12. value1 = []
13. value2 = []
14. value3 = []
15. for i in range(30):
16.     start_matrix = np.dot(start_matrix, transfer_matrix)
17.     value1.append(start_matrix[0][0])
18.     value2.append(start_matrix[0][1])
19.     value3.append(start_matrix[0][2])
       print(start_matrix)

[[0.23076935 0.30769244 0.46153864]]

1.#进行可视化
2. x = np.arange(30)
3. plt.plot(x,value1,label='cheerful')
4. plt.plot(x,value2,label='so-so')
5. plt.plot(x,value3,label='sad')
6. plt.legend()
7. plt.show()

可以发现，从10轮左右开始，我们的状态概率分布就不变了，一直保持在 [0.23076934,0.30769244,0.4615386]

参考：https://zhuanlan.zhihu.com/p/37121528

M-H采样python实现

https://zhuanlan.zhihu.com/p/37121528

假设目标平稳分布是一个均值3，标准差2的正态分布，而选择的马尔可夫链状态转移矩阵的条件转移概率是以为均值,方差1的正态分布在位置

的值。

from scipy.stats import norm
def norm_dist_prob(theta):
    y = norm.pdf(theta, loc=3, scale=2)
    return y
T = 5000
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, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1])))  #alpha值

    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 = 50
plt.hist(pi,
         num_bins,
         density=1,
         facecolor='red',
         alpha=0.7,
         label='Samples Distribution')
plt.legend()
plt.show()

二维Gibbs采样实例python实现

from mpl_toolkits.mplot3d import Axes3D
from scipy.stats import multivariate_normal

samplesource = multivariate_normal(mean=[5,-1], cov=[[1,0.5],[0.5,2]])

def p_ygivenx(x, m1, m2, s1, s2):
    return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt(1 - rho ** 2) * s2))

def p_xgiveny(y, m1, m2, s1, s2):
    return (random.normalvariate(m1 + rho * s1 / s2 * (y - m2), math.sqrt(1 - rho ** 2) * s1))

N = 5000
K = 20
x_res = []
y_res = []
z_res = []
m1 = 5
m2 = -1
s1 = 1
s2 = 2

rho = 0.5
y = m2

for i in range(N):
    for j in range(K):
        x = p_xgiveny(y, m1, m2, s1, s2)   #y给定得到x的采样
        y = p_ygivenx(x, m1, m2, s1, s2)   #x给定得到y的采样
        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,density=1, facecolor='green', alpha=0.5,label='x')
plt.hist(y_res, num_bins, density=1, facecolor='red', alpha=0.5,label='y')
plt.title('Histogram')
plt.legend()
plt.show()