复现经典:《统计学习方法》第19章 马尔可夫链蒙特卡罗法
第19章 马尔可夫链蒙特卡罗法
本文是李航老师的《统计学习方法》一书的代码复现。作者:黄海广
备注:代码都可以在github中下载。我将陆续将代码发布在公众号“机器学习初学者”,可以在这个专辑在线阅读。
蒙特卡罗法是通过基于概率模型的抽样进行数值近似计算的方法,蒙特卡罗法可以用于概率分布的抽样、概率分布数学期望的估计、定积分的近似计算。
随机抽样是蒙特卡罗法的一种应用,有直接抽样法、接受拒绝抽样法等。接受拒绝法的基本想法是,找一个容易抽样的建议分布,其密度函数的数倍大于等于想要抽样的概率分布的密度函数。按照建议分布随机抽样得到样本,再按要抽样的概率分布与建议分布的倍数的比例随机决定接受或拒绝该样本,循环执行以上过程。
马尔可夫链蒙特卡罗法数学期望估计是蒙特卡罗法的另一种应用,按照概率分布
抽取随机变量的个独立样本,根据大数定律可知,当样本容量增大时,函数的样本均值以概率1收敛于函数的数学期望
计算样本均值
,作为数学期望的估计值。
马尔可夫链是具有马尔可夫性的随机过程
通常考虑时间齐次马尔可夫链。有离散状态马尔可夫链和连续状态马尔可夫链,分别由概率转移矩阵
和概率转移核定义。
满足
或的状态分布称为马尔可夫链的平稳分布。
马尔可夫链有不可约性、非周期性、正常返等性质。一个马尔可夫链若是不可约、非周期、正常返的,则该马尔可夫链满足遍历定理。当时间趋于无穷时,马尔可夫链的状态分布趋近于平稳分布,函数的样本平均依概率收敛于该函数的数学期望。
可逆马尔可夫链是满足遍历定理的充分条件。
马尔可夫链蒙特卡罗法是以马尔可夫链为概率模型的蒙特卡罗积分方法,其基本想法如下:
(1)在随机变量
的状态空间上构造一个满足遍历定理条件的马尔可夫链,其平稳分布为目标分布;
(2)由状态空间的某一点
出发,用所构造的马尔可夫链进行随机游走,产生样本序列;
(3)应用马尔可夫链遍历定理,确定正整数
和$n(m<n)$,得到样本集合$\{ 1="" 2="" x="" _="" {="" m="" +="" }="" ,="" \cdots="" n="" \}$,进行函数$f(x)$的均值(遍历均值)估计:<="" p="">
Metropolis-Hastings算法是最基本的马尔可夫链蒙特卡罗法。假设目标是对概率分布
进行抽样,构造建议分布,定义接受分布进行随机游走,假设当前处于状态,按照建议分布机抽样,按照概率接受抽样,转移到状态 ,按照概率拒绝抽样,停留在状态,持续以上操作,得到一系列样本。这样的随机游走是根据转移核为的可逆马尔可夫链(满足遍历定理条件)进行的,其平稳分布就是要抽样的目标分布。吉布斯抽样(Gibbs sampling)用于多元联合分布的抽样和估计吉布斯抽样是单分量 Metropolis-Hastings算法的特殊情况。这时建议分布为满条件概率分布
吉布斯抽样的基本做法是,从联合分布定义满条件概率分布,依次从满条件概率分布进行抽样,得到联合分布的随机样本。假设多元联合概率分布为
,吉布斯抽样从一个初始样本出发,不断进行迭代,每一次迭代得到联合分布的一个样本,在第次迭代中,依次对第个变量按照满条件概率分布随机抽样,得到最终得到样本序列。
蒙特卡洛法(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)。
假设有随机变量
, 取值 , 其概率密度函数为 , 为定义在 上的函数, 目标是求函数 关于密度函数 的数学期望 。
针对这个问题,蒙特卡罗法按照概率分布
独立地抽取 个样本,比如用以上的抽样方法,之后计算函
数
的样本均值
作为数学期望
近似值。
根据大数定律可知, 当样本容量增大时, 样本均值以概率1收敛于数学期望:
这样就得到了数学期望的近似计算方法:
马尔可夫链
考虑一个随机变量的序列
这里 ,表示时刻 的随机变量, .
每个随机变量
的取值集合相同, 称为状态空间, 表示为. 随机变量可以是离散的, 也可以是连续的。
以上随机变量的序列构成随机过程(stochastic process)。
假设在时刻
的随机变量 遵循概率分布 ,称为初始状态分布。在某个时刻 的随机变量 与前
一个时刻的随机变量
之间有条件分布 如果 只依赖于 , 而不依赖于过去的随机变量
,
,
这一性质称为马尔可夫性,即
具有马尔可夫性的随机序列
称为马尔可夫链, 或马尔可夫过程(Markov process)。 条件概率分布
称为马尔可夫链的转移概率分布。 转移概率分布决定了马尔可夫裢的特性。
平稳分布
设有马尔可夫链
,其状态空间为 ,转移概率矩阵为 , 如果存在状态空间 上的一个分布
使得
则称丌为马尔可夫裢
的平稳分布。
直观上,如果马尔可夫链的平稳分布存在,那么以该平稳分布作为初始分布,面向未来进行随机状态转移,之后任何一个时刻的状态分布都是该平稳分布。
引理19.1
给定一个马尔可夫链
, 状态空间为, 移概率矩阵为, 则分布 为 的平稳分布的充要条件是是下列方程组的解:
吉布斯采样
输入: 目标概率分布的密度函数
, 函数;
输出:
的随机样本 ,函数样本均值 ;
参数: 收敛步数
, 迭代步数 .
初始化。给出初始样本
($x^{0}{1}, x^{0}{2},..., x^{0}_{k}^{T}$.对
循环执行
设第次迭代结束前的样本为($x^{i-1}{1}, x^{i-1}{2},..., x^{i-1}_{k}^{T}i$次迭代进行如下几步操作:
(1)由满条件分布
抽取...
(j)由满条件分布 抽取
(k)由满条件分布
抽取
得到第
次迭代值 .
得到样本集合
{
}
计算
网络资源:
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
import random
import math
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as nptransfer_matrix = np.array([[0.6, 0.2, 0.2], [0.3, 0.4, 0.3], [0, 0.3, 0.7]],dtype='float32')
start_matrix = np.array([[0.5, 0.3, 0.2]], dtype='float32')value1 = []
value2 = []
value3 = []
for i in range(30):start_matrix = np.dot(start_matrix, transfer_matrix)value1.append(start_matrix[0][0])value2.append(start_matrix[0][1])value3.append(start_matrix[0][2])
print(start_matrix)
[[0.23076935 0.30769244 0.46153864]]
#进行可视化
x = np.arange(30)
plt.plot(x,value1,label='cheerful')
plt.plot(x,value2,label='so-so')
plt.plot(x,value3,label='sad')
plt.legend()
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 + 1pi_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_normalsamplesource = 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 = 2rho = 0.5
y = m2for 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()
本章代码来源:https://github.com/hktxt/Learn-Statistical-Learning-Method
下载地址
https://github.com/fengdu78/lihang-code
参考资料:
[1] 《统计学习方法》: https://baike.baidu.com/item/统计学习方法/10430179
[2] 黄海广: https://github.com/fengdu78
[3] github: https://github.com/fengdu78/lihang-code
[4] wzyonggege: https://github.com/wzyonggege/statistical-learning-method
[5] WenDesi: https://github.com/WenDesi/lihang_book_algorithm
[6] 火烫火烫的: https://blog.csdn.net/tudaodiaozhale
[7] hktxt: https://github.com/hktxt/Learn-Statistical-Learning-Method
复现经典:《统计学习方法》第19章 马尔可夫链蒙特卡罗法相关推荐
- 李航《统计学习方法》第二章课后答案链接
李航<统计学习方法>第二章课后答案链接 李航 统计学习方法 第二章 课后 习题 答案 http://blog.csdn.net/cracker180/article/details/787 ...
- 李航《统计学习方法》第一章课后答案链接
李航<统计学习方法>第一章课后答案链接 李航 统计学习方法 第一章 课后 习题 答案 http://blog.csdn.net/familyshizhouna/article/detail ...
- 统计学习方法——第1章(个人笔记)
统计学习方法--第1章 统计学习及监督学习概论 <统计学习方法>(第二版)李航,学习笔记 1.1 统计学习 1.特点 (1)以计算机及网络为平台,是建立在计算机及网络上的: (2)以数据为 ...
- 统计学习方法笔记第二章-感知机
统计学习方法笔记第二章-感知机 2.1 感知机模型 2.2感知机学习策略 2.2.1数据集的线性可分型 2.2.2感知机学习策略 2.3感知机学习算法 2.3.1感知机算法的原始形式 2.3.2算法的 ...
- 统计学习方法 | 第7章 支持向量机
第7章 支持向量机 <统计学习方法>Python代码实现 [转载自Github开源项目]https://github.com/fengdu78/lihang-code 1.支持向量机最简单 ...
- 一篇详解带你再次重现《统计学习方法》——第二章、感知机模型
个性签名:整个建筑最重要的是地基,地基不稳,地动山摇. 而学技术更要扎稳基础,关注我,带你稳扎每一板块邻域的基础. 博客主页:七归的博客 专栏:<统计学习方法>第二版--个人笔记 创作不易 ...
- 统计学习方法 - 第1章 - 概论
全书章节 第1章 统计学习方法概论 第2章 感知机 第3章 k近邻法 第4章 朴素贝叶斯法 第5章 决策树 第6章 逻辑斯谛回归与最大熵模型 第7章 支持向量机 第8章 提升方法 第9章 EM算法及其 ...
- 机器学习——统计学习方法——第1章 统计学习及监督学习概论
监督学习是从标注数据中学习模型的机器学习问题,是统计学习的重要组成部分. 1.1 统计学习 统计学习的特点 统计学习是关于计算机基于数据构建概率统计模型并运用模型对数据进行预测与分析的一门学科. 特点 ...
- 《统计学习方法》第一章总结
统计学习是关于计算机基于数据构建概率统计模型并运用模型对数据进行预测与分析的一门学科. 统计学习的对象是数据.目的是对数据进行预测和分析. 统计学习关于数据的基本假设是同类数据具有一定的统计规律性,这 ...
最新文章
- [恢]hdu 2015
- 谷歌搜索:几乎所有的英文搜索都用上BERT了
- 最大连续子序列乘积(DP)
- 计算机原理及基础 —— 有符号类型和无符号类型
- wpf 开源框架_.NET Core跨平台基础框架:10 篇热文汇总
- OpenGL整体概念
- java中二进制怎么说_面试常用:说清楚Java中synchronized和volatile的区别
- stata导入数据问题
- idea 如何关闭 field injection is not recommended 警告
- 机器学习算法汇总:人工神经网络、深度学习及其它
- 拖拽图片到另一个div里
- 巴塞罗那IoT“首秀”归来,新华三成功展现物联网风采
- nlp中region_百度5年深耕NLP 他把聊天机器人变成你的“全科医生”
- c语言 枚举大小写,C语言枚举类型(Enum)深入理解
- 今日头条的个性化推荐
- 这些免费版音视频格式转换器哪个最好用
- Jquery引用在线CDN公共资源库
- Educoder---Java继承与接口、文件
- fore-each操作数组
- 随身Q代理服务器大升级