MCMC中的Metropolis–Hastings算法与吉布斯采样
来源:互联网 发布:淘宝如何查看收藏宝贝 编辑:程序博客网 时间:2024/05/21 06:47
Metropolis–Hastings算法是一种具体的MCMC方法,而吉布斯采样(Gibbs Sampling)是Metropolis–Hastings算法的一种特殊形式。二者在机器学习中具有重要作用,Bishop在他的机器学习经典之作PRML中也专门用了一章的篇幅来介绍随机采样方法。本文将结合R语言实例来探讨这两种算法的相关话题。本文是这个系列文章的最后一篇,主要介绍MCMC中的Metropolis–Hastings算法和吉布斯采样(Gibbs Sampling)方面的内容。
作为本系列文章的组成部分,也作为你阅读本文所必须的预备知识,希望各位读者确认已经对如下文章所谈之话题了然于心:
- 蒙特卡洛(Monte Carlo)法求定积分
- 蒙特卡洛采样之拒绝采样(Reject Sampling)
- 矩阵的极限与马尔科夫链(上)
- 矩阵的极限与马尔科夫链(下)
- Chapman-Kolmogorov Equation
- 蒙特卡洛马尔科夫链(MCMC)
欢迎关注白马负金羁的CSDN博客 http://blog.csdn.net/baimafujinji,为保证公式、图表得以正确显示,强烈建议你从该地址上查看原版博文。本博客主要关注方向包括:数字图像处理、算法设计与分析、数据结构、机器学习、数据挖掘、统计分析方法、自然语言处理。
Metropolis–Hastings算法
In statistics, the Metropolis–Hastings algorithm is a MCMC method for obtaining a sequence of random samples from a probability distribution for which direct sampling is difficult.
听到Metropolis–Hastings这个名字,有没有点似曾相识,或者在其他地方也听过。如果你做过数学建模,或者学过优化算法方面的内容,那么在研究模拟退火算法时就应该听过Metropolis–Hastings算法这个名字。本质上来说,模拟退火中的Metropolis–Hastings和MCMC中的Metropolis–Hastings确实是一回事。The original algorithm used in simulated annealing and MCMC’s is due to Metropolis. Later generalized by Hastings. Hastings showed that it is not necessary to use a symmetric proposal distribution, and proposed that the proposed new state can be generated from any
Metropolis–Hastings算法的执行步骤如下,注意其中的q就是前面文字描述中的P,即后验概率:
Metropolis–Hastings算法的执行步骤是先随便指定一个初始的样本点
既然已经从理论上证明Metropolis–Hastings算法的确实可行,下面我们举一个简单的例子来看看实际中Metropolis–Hastings算法的效果。稍微有点不一样的地方是,我们这里并未出现
柯西分布也叫作柯西-洛伦兹分布,它是以奥古斯丁·路易·柯西与亨德里克·洛伦兹名字命名的连续概率分布,其概率密度函数为
其中
下面是在R中进行基于Metropolis–Hastings算法对柯西分布进行采样的示例代码。
rm(list=ls()) ## 清除全部对象set.seed(201912)reps=40000# target density: the cauchy distributioncauchy<-function(x, x0=0, gamma=1){ out<-1/(pi*gamma*(1+((x-x0)/gamma)^2)) return(out)}chain<-c(0)for(i in 1:reps){ proposal<-chain[i]+runif(1, min=-1, max=1) accept<-runif(1)<cauchy(proposal)/cauchy(chain[i]) chain[i+1]<-ifelse(accept==T, proposal, chain[i])}# plot over timeplot(chain, type="l")plot(density(chain[1000:reps]), xlim=c(-5,5), ylim=c(0,0.4), col="red")den<-cauchy(seq(from=-5, to=5, by=0.1), x=0, gamma=1)lines(den~seq(from=-5, to=5, by=0.1), lty=2, col="blue")
下图是每次采样点的数值分布情况。
为了更清晰明确的看到采样结果符合预期分布,我们也绘制了分布的密度图,如下图所示,红色实线是采样分布的密度图,而蓝色虚线则是实际柯西分布的密度图,可见二者吻合地相当好。
当然,作为一个简单的开始,上面的例子中我们并没有涉及到
rm(list=ls()) ## 清除全部对象f <- function(x, sigma) { if (any(x < 0)) return (0) stopifnot(sigma > 0) return((x / sigma^2) * exp(-x^2 / (2*sigma^2)))}m <- 40000sigma <- 4x <- numeric(m)x[1] <- rchisq(1, df=1)k <- 0u <- runif(m)for (i in 2:m) { xt <- x[i-1] y <- rchisq(1, df = xt) num <- f(y, sigma) * dchisq(xt, df = y) den <- f(xt, sigma) * dchisq(y, df = xt) if (u[i] <= num/den) x[i] <- y else { x[i] <- xt k <- k+1 #y is rejected }}
然后我们要验证一下,我的生产的采样数据是否真的符合瑞利分布。注意在R中使用瑞利分布的相关函数,需要加装VGAM包。下面的代码可以让我们直观地感受到采样结果的分布情况。
> curve(drayleigh(x, scale = 4, log = FALSE), from = -1, to = 15, xlim=c(-1,15), ylim=c(0,0.2), + lty = 2, col="blue",xlab = "", ylab="")> par(new=TRUE)> plot(density(x[1000:m]) , xlim=c(-1,15), ylim=c(0,0.2),col="red")
从下图中不难看出,我们的采样点分布确实符合预期。当然,这仅仅是一个演示用的小例子。显然,它并不高效。Œ因为我们采用自由度为
吉布斯采样(Gibbs Sampling)
In statistics, Gibbs sampling is a MCMC algorithm for obtaining a sequence of observations which are approximated from a specified multivariate probability distribution, when direct sampling is difficult. 吉布斯采样 can be seen as a special case of the Metropolis-Hastings algorithm. 而这个特殊之处就在于,我们在使用吉布斯采样时,通常是对多变量分布进行采样的。比如说,现在我们有一个分布
当然如果要从理论上证明吉布斯采样确实可以得到预期的分布,我们就应该考察它是否满足detailed balance。但是前面也讲过,吉布斯采样是Metropolis–Hastings算法的一个特例。所以其实我们甚至无需费力去考察其是否满足detailed balance,如果能够证明吉布斯采样就是Metropolis–Hastings算法,而Metropolis–Hastings算法是满足detailed balance,其实也就得到了我们想要的。下面是证明的过程,其中
吉布斯采样的例子网上有很多资料可以参考,这里就不再赘述了。最后说明一点就是,以吉布斯采样为基础实现的MCMC在自然语言处理中的LDA里有重要应用。我们花费了很长的篇幅,用了一个系列的文章最后才把MCMC介绍完毕,如果你正在或者后续准备研究LDA的话,这些内容将是非常必要的基础,也是相当不错的开端。
参考文献
[1] 本文中的英文定义来自于维基百科
[2] http://zhfuzh.blog.163.com/blog/static/1455393872012822854853/
[3] 悉尼科技大学徐亦达博士的机器学习公开课授课材料
[4] 莱斯大学Justin Esarey助理教授(http://jee3.web.rice.edu/index.htm)的公开课资料(https://www.youtube.com/watch?v=j4nEAqUUnVw)
[5] http://staff.ustc.edu.cn/~zwp/teach/Stat-Comp/Lec8.R
[6] Christopher Bishop, Pattern Recognition and Machine Learning, Springer, 2007
- MCMC中的Metropolis–Hastings算法与吉布斯采样
- MCMC与Metropolis-Hastings 算法
- MCMC马尔可夫蒙特卡洛方法总结-Gibss吉布斯采样, Metropolis-Hastings方法
- 采样之Metropolis Hastings算法
- Metropolis–Hastings算法
- Metropolis-Hastings MCMC algorithms
- 随机过程–Metropolis-Hastings算法
- Metropolis-Hastings 算法
- Metropolis-Hastings算法理解
- MCMC: Metropolis-Hastings, Gibbs and slice sampling
- 抽样,mcmc, Metropolis-Hastings,Gibbs Sampling
- Metropolis 采样与蒙特卡洛算法
- MCMC(马尔科夫蒙特卡洛方法)和MH(Metropolis-Hastings)算法的一个python实现
- 随机过程--Metropolis-Hastings算法
- 马尔科夫链蒙特卡洛算法(MCMC)中的Gibbs采样
- MCMC采样算法理解
- PRML读书会第十一章 Sampling Methods(MCMC, Markov Chain Monte Carlo,Metropolis-Hastings,Gibbs Sampling)
- 初探密码破译:Metropolis-Hastings算法破解密文
- MFC简易音乐播放器
- 每天对着点电脑坐久了
- SSH框架部署到服务器上js中url中文乱码问题
- RERAN:安卓系统的定时和点击的录制和回放——(5)
- spring boot 配合前端实现跨域
- MCMC中的Metropolis–Hastings算法与吉布斯采样
- LintCode 52-下一个排列
- OpenGL入门学习
- Ubuntu 16.04安装/更新VMware tools&设置共享文件夹方法
- USACO 2.1 castle <图搜索>*****
- 互斥锁mutex
- 100. Same Tree 难度:easy
- 数据结构与算法(C语言版)__红黑树05
- 51nod 1489 蜥蜴和地下室(dfs or dp)