【机器学习基础】数学推导+纯Python实现机器学习算法21:马尔可夫链蒙特卡洛...
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 yT = 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]))) # 从均匀分布中随机抽取一个数uu = 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();
MCMC与贝叶斯推断
MCMC对高效贝叶斯推断有着重要的作用。假设有如下贝叶斯后验分布推导:
在概率分布多元且形式复杂的情形下,经过贝叶斯先验和似然推导后(即右边的分式),很难进行积分运算。具体包括以下三种积分运算:规范化、边缘化和数学期望。
首先是上式中的分母,即规范化计算:
如果是多元随机变量或者是包含隐变量
,后验分布还需要边缘化计算 :
另外如有一个函数
,可以计算该函数关于后验概率分布的数学期望:
当观测数据、先验分布和似然函数都比较复杂的时候,以上三个积分计算都会变得极为困难,这也是早期贝叶斯推断受到冷落的一个原因。后来MCMC方法兴起,Bayesian+MCMC的一套做法逐渐流行起来。
参考资料:
统计学习方法第二版
https://zhuanlan.zhihu.com/p/37121528
往期精彩:
数学推导+纯Python实现机器学习算法20:LDA线性判别分析
数学推导+纯Python实现机器学习算法19:PCA降维
数学推导+纯Python实现机器学习算法18:奇异值分解SVD
数学推导+纯Python实现机器学习算法17:XGBoost
数学推导+纯Python实现机器学习算法16:Adaboost
数学推导+纯Python实现机器学习算法15:GBDT
数学推导+纯Python实现机器学习算法14:Ridge岭回归
数学推导+纯Python实现机器学习算法13:Lasso回归
数学推导+纯Python实现机器学习算法12:贝叶斯网络
数学推导+纯Python实现机器学习算法11:朴素贝叶斯
数学推导+纯Python实现机器学习算法10:线性不可分支持向量机
数学推导+纯Python实现机器学习算法8-9:线性可分支持向量机和线性支持向量机
数学推导+纯Python实现机器学习算法7:神经网络
数学推导+纯Python实现机器学习算法6:感知机
数学推导+纯Python实现机器学习算法5:决策树之CART算法
数学推导+纯Python实现机器学习算法4:决策树之ID3算法
数学推导+纯Python实现机器学习算法3:k近邻
数学推导+纯Python实现机器学习算法2:逻辑回归
数学推导+纯Python实现机器学习算法1:线性回归
往期精彩回顾适合初学者入门人工智能的路线及资料下载机器学习及深度学习笔记等资料打印机器学习在线手册深度学习笔记专辑《统计学习方法》的代码复现专辑
AI基础下载机器学习的数学基础专辑获取一折本站知识星球优惠券,复制链接直接打开:https://t.zsxq.com/yFQV7am本站qq群1003271085。加入微信群请扫码进群:
【机器学习基础】数学推导+纯Python实现机器学习算法21:马尔可夫链蒙特卡洛...相关推荐
- 【机器学习基础】数学推导+纯Python实现机器学习算法30:系列总结与感悟
Python机器学习算法实现 Author:louwill Machine Learning Lab 终于到了最后的总结.从第一篇线性回归的文章开始到现在,已经接近有两年的时间了.当然,也不是纯写这3 ...
- 【机器学习基础】数学推导+纯Python实现机器学习算法24:HMM隐马尔可夫模型
Python机器学习算法实现 Author:louwill Machine Learning Lab HMM(Hidden Markov Model)也就是隐马尔可夫模型,是一种由隐藏的马尔可夫链随机 ...
- 【机器学习基础】数学推导+纯Python实现机器学习算法28:CRF条件随机场
Python机器学习算法实现 Author:louwill Machine Learning Lab 本文我们来看一下条件随机场(Conditional Random Field,CRF)模型.作为概 ...
- 【机器学习基础】数学推导+纯Python实现机器学习算法27:EM算法
Python机器学习算法实现 Author:louwill Machine Learning Lab 从本篇开始,整个机器学习系列还剩下最后三篇涉及导概率模型的文章,分别是EM算法.CRF条件随机场和 ...
- 【机器学习基础】数学推导+纯Python实现机器学习算法26:随机森林
Python机器学习算法实现 Author:louwill Machine Learning Lab 自从第14篇文章结束,所有的单模型基本就讲完了.而后我们进入了集成学习的系列,整整花了5篇文章的篇 ...
- 【机器学习基础】数学推导+纯Python实现机器学习算法25:CatBoost
Python机器学习算法实现 Author:louwill Machine Learning Lab 本文介绍GBDT系列的最后一个强大的工程实现模型--CatBoost.CatBoost与XGBoo ...
- 【机器学习基础】数学推导+纯Python实现机器学习算法24:LightGBM
Python机器学习算法实现 Author:louwill Machine Learning Lab 第17讲我们谈到了竞赛大杀器XGBoost,本篇我们来看一种比XGBoost还要犀利的Boosti ...
- 【机器学习基础】数学推导+纯Python实现机器学习算法23:kmeans聚类
Python机器学习算法实现 Author:louwill Machine Learning Lab 聚类分析(Cluster Analysis)是一类经典的无监督学习算法.在给定样本的情况下,聚类分 ...
- 【机器学习基础】数学推导+纯Python实现机器学习算法22:最大熵模型
Python机器学习算法实现 Author:louwill Machine Learning Lab 最大熵原理(Maximum Entropy Principle)是一种基于信息熵理论的一般原理,在 ...
最新文章
- LeetCode简单题之字符串转化后的各位数字之和
- 简单C++线程池包装类源码示例
- selenium grid2 使用远程机器的浏览器
- 10 张图打开 CPU 缓存一致性的大门
- 电大最全计算机应用技术基础答案,电大最新最全计算机应用技术基础答案100%通过率...
- ios刷android8.0,颤抖吧 iOS, Android 8.0正式发布!
- 母版中menu控件上传后出现脚本错误
- 软件行业选择大公司还是小公司
- 关于fineui中在gird中插入按钮的知识
- docker server 容器连接sql_借力 Docker ,三分钟搞定 MySQL 主从复制!
- java -classpath or -cp 的设置和解释
- 100天,Python从入门到精通!
- 过年了,没买到票的可以过来看一看,12306抢票助手
- ubuntu下研华工控机CAN卡驱动的安装与测试
- 移动App测试中的最佳做法
- 大数据学习基础知识总纲
- 炼数成金hadoop视频干货06-10
- linux增加预读缓存区大小,Linux使用blockdev命令调整文件预读大小的方法
- zjyxmdshoes
- 计算机会计u8实验报告,会计电算化用友实验报告-20210406233157.pdf-原创力文档