R语言贝叶斯Metropolis-Hastings采样 MCMC算法理解和应用可视化案例

简介: R语言贝叶斯Metropolis-Hastings采样 MCMC算法理解和应用可视化案例

贝叶斯MCMC模拟是一个丰富的领域,涵盖了各种算法,共同目标是近似后验模型点击文末“阅读原文”获取完整代码数据

相关视频

image.png

例如,使用的rstan包采用了一个Hamiltonian Monte Carlo算法。用于贝叶斯建模的另一个rjags包采用了Gibbs sampling算法。尽管细节有所不同,但这两种算法都是基于基本的Metropolis-Hastings算法的变体。

主要思想

考虑以下数值结果为Y的正态-正态模型,其围绕未知均值μ的标准差为0.75:

ab4f331cf424a57e73f4587988539dd9.png

相应的似然函数L(μ|y)和先验概率密度函数f(μ)对于y∈(−∞,∞)和μ∈(−∞,∞)是:

50f1fb0ec77de4fa0fbbb3b6de0bbf32.png

假设我们观察到一个结果Y=6.25。μ的后验模型是具有均值4和标准差0.6的正态分布:

8163e930f9f456d936f342e1d74a5658.png

如果我们无法指定μ的后验模型(假装一下),我们可以使用MCMC模拟来近似它。为了了解这是如何工作的,考虑一个潜在的N=5000次迭代的MCMC模拟结果。将这个马尔可夫链{μ(1),μ(2),…,μ(N)}视为μ的后验可能值范围的一次游览,你可以将自己看作是导游。左侧的轨迹图显示了游览路线或游览站点的顺序,μ(i)。右侧的直方图显示了你在每个μ区域中停留的相对时间。

93887157dd47625bd4ab375f55a14378.png

我们可以通过简单地使用rnorm()从N(4,0.62)后验中直接抽样来实现这个算法。结果是来自后验的一个很好的独立样本,这反过来又产生了一个准确的后验近似:

set.seed(84375)
mc_tour <- data.frame(mu = rnorm(5000, mean = 4, sd = 0.6))
ggplot(mc_tour, aes(x = mu)) + 
  geom_histogram(aes(y = ..density..), color = "white", bins = 15) + 
  stat_function(fun = dnorm, args = list(4, 0.6), color = "blue")

5e52f9c7d19d0a98b7ea93708086e411.png


点击标题查阅往期内容


R语言贝叶斯MCMC:GLM逻辑回归、Rstan线性回归、Metropolis Hastings与Gibbs采样算法实例


5ae4a05b3da6628dd3a28ea9a7b47831.png


Metropolis-Hastings算法

用于构建马尔可夫链{μ(1),μ(2),...,μ(N)}的Metropolis-Hastings算法在这里得到了规范。我们将始终乐意将我们移动到更合理的后验区域。例如,假设我们的马尔可夫游览当前位于位置“3”:

current <- 3

然后,要确定下一个游览站点,我们首先通过从Unif(current - 1, current + 1)模型中随机抽样来提出一个位置(步骤1):

set.seed(8)
proposal

6a737659fd2ad61ee6ace1865421a904.png

重新查看图,我们观察到所提议的点(2.93)的(未归一化的)后验可能性略微低于当前点(3)。我们可以使用dnorm()计算这两个μ的未归一化后验可能性,即f(μ|y=6.25)∝f(μ)L(μ|y=6.25)。

prolaus <- dnorm(osal, 0, 1) * dnorm(6.25, proal, 0.75)
curaus  <- dnorm(curnt, 0, 1) * dnorm(6.25, nt, 0.75)

45602c7fbc4bfed06a255dc904c97900.png

de54340d44393e1dd88d2a38e3a7403d.png

由此可见,尽管不确定,接受并随后移动到所提议位置的概率α的相对高:

alpha <- min(1, prplaus / curr_plaus)
alpha

9ed9eb7518bac8228082fcdb9c11e569.png

为了做出最终决定,我们设置了一个加权硬币,它以概率α(0.824)接受提议,并以概率1−α(0.176)拒绝提议。通过sample()函数随机抛掷这个硬币,我们接受了提议,意味着下一点是2.933:

nex_sop <- sape(c(roposl, curet),
                    size = 1, ob = c(alha, 1-apha))
net_top

4d4b536b2722187b5b9bb60745d2c46c.png

这仅仅是对我们的正态后验进行一次Metropolis-Hastings算法迭代的无数可能结果之一。为了简化这个过程,我们将编写自己的R函数one_mh_iteration(),该函数实现从任何给定的当前点开始的单个Metropolis-Hastings迭代,并利用具有任意半宽度w的均匀提议模型。

我们首先指定one_mhiteation是两个参数的function():均匀分布的范围w和当前链值current

在函数内部,我们执行与上面相同的步骤,并return()三个信息:proposal、接受概率alphanext_stop

one_meration <- fntion(w, crret){
 # 第一步:提议下一个链位置
 prpsal <- runif(1, nt - w, max = crrent+ w)
  
 # 第二步:决定是否移动到新位置
 prop_plau <- dnorm(prpoal, 0, 1) * dnorm(625, proosl, 0.75)
 curetplaus  <- dnorm(cunt, 0, 1) * dnorm(6.25, urent, 0.75)
 alpha <- min(1, propoalplaus / current_plaus)
 
 # 返回结果
 return(dta.fame(proposa, alph, next_stp))
}

让我们试一试。在种子为8的情况下从当前点3运行onemh_ieration()可以复制上面的结果:

set.seed(8)
one_h_itraton(w = 1, current = 3)

16465c04f2cef3627f918acfd469c6be.png

如果我们使用83这个种子,所提议的下一个点是2.018,对应的接受概率很低,只有0.017:

set.eed(83)
one_mh_itraton(w =1, crrent = 3)


fc5002b030089e43d00b2eb6fe9feaf4.png

这是有道理的。从图中可以看出,2.018的后验可能性远低于我们当前所在地3的后验可能性。虽然我们确实想要探索这样极端的值,但我们不想经常这样做。事实上,在我们投掷硬币时,提议被拒绝,并且旅游将再次访问位置3。

我们可以确认当所提议的下一个旅游站点(这里是3.978)的后验可能性大于我们当前位置时,接受概率为1,提议将自动接受:

set.seed(7)
onemh_ieratonw = 1, urrnt = 3)


5e11cf1ed714d8655cc5a29ebb6d938c.png

实施Metropolis-Hastings

现在,我们需要一遍又一遍地重复上一个单个迭代的过程来构建一个由N(4,0.62)后验分布组成的Metropolis-Hastings。下面的mh_tour()函数可以构建一个给定长度N的Metropolis-Hastings,利用任何给定半宽度w的均匀提议模型:

mh_tour <- function(N, w){
  # 1. 在位置3开始
  curet <- 3
  # 2. 初始化模拟
  mu <- rep(0, N)
  # 3. 模拟N个马尔科夫链停止
  for(i in 1:N){    
    # 模拟一个迭代
    sim <- onemh_itrationw = w,currnt  crrent)
    
    # 记录下一个位置
    mu[i] <- si$xt_top
    
    # 重置当前位置
    current <- sim$nt_tp
  }
  
  # 4. 返回位置
  retun(datfrae(tratio= c(:N,mu))
}

在调用此函数时:

  1. 在位置3处开始,这是一个基于我们对μ的先前了解而做出的相当任意的选择。
  2. 通过设置“空”向量来初始化模拟,我们最终将在其中存储N次停止的位置(mu)。
  3. 利用for循环,在1到N的每个停留点i中运行on_m_iteaion(),并将结果的next_stop存储在mu向量的第i个元素中。在关闭for循环之前,更新current停止以作为下一次迭代的起点。
  4. 返回具有迭代号和相应停留点mu的数据框。

要查看此函数的实际应用,请使用m_our()模拟长度为N = 5000的Markov链,利用半宽度w=1的均匀提议模型:

set.seed(84735)
mh_sulio_1 <- m_our(N = 5000, w = 1)

下面显示了结果的跟踪图和直方图。值得注意的是,产生了对N(4,0.62)后验分布的极其准确的近似。通过严格过程,我们利用了均匀分布的相关采样来近似正态分布。

ggplo= 20) + 
  stat_fuctin(fun = dnorm,args = ist(4,0.6), color = "blue")

c4eb15710d562d2bb7d5d2874d09bf71.png

Beta-Binomial模型例子

让我们实现Metropolis-Hastings算法来处理一个Beta-Binomial模型,其中我们观察到在2次尝试中有1次成功,即 Y=1,n=2: 。

603a23a565d828b21bce43deb11d5ee6.png

同样,假设我们只能将后验概率密度上定义到某些缺失的归一化常数,

23cf8cf5f2ea6b8637152c8813ed1938.png

下面的oneiertion()函数实现了该独立采样算法的单次迭代,从任何给定的当前值π开始,并对给定的a和b使用Beta(a,b)建议模型。在计算接受概率α时,请注意我们使用dbeta()来评估先验概率密度函数和建议概率密度函数,以及使用dbinom()来评估具有数据Y=1,n=2,π的二项式似然函数:

one_terton <- function(a, , curnt){
# 第 1 步:提出下一个链位置
 proposal <- rbeta(1, a, b)
  
 # 第 2 步:决定
 
 popoal_plas <- deta(prposa, 2, 3) * dbno(1, 2, prosal)
 proposlq     <- dbeta(prpsal, a, b)
 curentplus  <- dbea(curet, 2, 3) * dbnom(1, 2, )
 curnt_q      <- dbeta(curent, a, b)

接下来,我们编写一个betbn_tour()函数,为任何Beta(a,b)建议模型构建一个N长度的Markov链遍历,利用one_iertion()来确定每个停止位置:


beai_tour <- fucton(N, a, b){
  # 1. 位置0.5开始链
  current <- 0.5
  # 2. 初始化模拟
  pi <- rp(0, N)
  
  # 3. 模拟N次Markov链停止
  for(i in 1:N){    
    # 模拟一次迭代
    sim <- on_ieaion(a = a, b = b, crent = curent)
    
    # 记录下一个位置
    pi[i] <- simnet_sop
    
    # 重置当前位置
    current <- simnetstop
  }
}

我们尝试不同调整的Beta(a,b)建议模型。在这里,为了简单起见,我们使用Beta(1,1),即Unif(0,1)的建议模型进行了5000步的Beta-Binomial后验分布遍历


set.seed(84735)
bebn_im <- betb_tour(N = 5000, a = 1, b = 1)
# 绘制结果
ggplot(beta

总结

本文建立了对基本的Metropolis-Hastings MCMC算法的概念理解。还实现了该算法来研究常见的正态-正态和Beta-Binomial模型。无论是在这些相对简单的单参数模型设置中,还是在更复杂的模型设置中,Metropolis-Hastings算法通过两个步骤之间的迭代产生了后验分布的近似样本:

  1. 通过从提议概率密度函数中抽取一个新的链位置来提出一个新的位置,这可能取决于当前位置。
  2. 确定是否接受提议。简单地说,我们是否接受提议取决于其后验可能性相对于当前位置的后验可能性有多有利。

相关文章
|
1月前
|
机器学习/深度学习 数据采集 算法
基于贝叶斯算法的手机垃圾短信过滤
基于贝叶斯算法的手机垃圾短信过滤
|
1月前
|
机器学习/深度学习 存储 算法
用kNN算法诊断乳腺癌--基于R语言
用kNN算法诊断乳腺癌--基于R语言
|
1月前
|
算法 Python
利用贝叶斯算法对简单应用实现预测分类
利用贝叶斯算法对简单应用实现预测分类
24 0
|
1月前
|
算法 项目管理
R语言实现蒙特卡洛模拟算法
R语言实现蒙特卡洛模拟算法
|
1月前
|
机器学习/深度学习 数据采集 算法
数据分享|R语言机器学习预测案例合集:众筹平台、机票折扣、糖尿病患者、员工满意度
数据分享|R语言机器学习预测案例合集:众筹平台、机票折扣、糖尿病患者、员工满意度
|
1月前
|
算法 搜索推荐
R语言混合SVD模型IBCF协同过滤推荐算法研究——以母婴购物平台为例
R语言混合SVD模型IBCF协同过滤推荐算法研究——以母婴购物平台为例
|
1月前
|
数据可视化
R语言离散时间马尔可夫链(Markov chain)模型分类案例可视化分析
R语言离散时间马尔可夫链(Markov chain)模型分类案例可视化分析
|
1天前
|
机器学习/深度学习 算法 数据可视化
m基于PSO-LSTM粒子群优化长短记忆网络的电力负荷数据预测算法matlab仿真
在MATLAB 2022a中,应用PSO优化的LSTM模型提升了电力负荷预测效果。优化前预测波动大,优化后预测更稳定。PSO借鉴群体智能,寻找LSTM超参数(如学习率、隐藏层大小)的最优组合,以最小化误差。LSTM通过门控机制处理序列数据。代码显示了模型训练、预测及误差可视化过程。经过优化,模型性能得到改善。
14 6
|
1天前
|
算法 调度
基于变异混合蛙跳算法的车间调度最优化matlab仿真,可以任意调整工件数和机器数,输出甘特图
**摘要:** 实现变异混合蛙跳算法的MATLAB2022a版车间调度优化程序,支持动态调整工件和机器数,输出甘特图。核心算法结合SFLA与变异策略,解决Job-Shop Scheduling Problem,最小化总完成时间。SFLA模拟蛙群行为,分组进行局部搜索和全局信息交换。变异策略增强全局探索,避免局部最优。程序初始化随机解,按规则更新,经多次迭代和信息交换后终止。
|
1天前
|
机器学习/深度学习 算法 计算机视觉
基于ADAS的车道线检测算法matlab仿真
**摘要:** 基于ADAS的车道线检测算法利用Hough变换和边缘检测在视频中识别车道线,判断车道弯曲情况,提供行驶方向信息,并高亮显示。在MATLAB2022a中实现,系统包括图像预处理(灰度化、滤波、边缘检测)、车道线特征提取(霍夫变换、曲线拟合)和车道线跟踪,确保在实时场景中的准确性和稳定性。预处理通过灰度转换减少光照影响,滤波去除噪声,Canny算法检测边缘。霍夫变换用于直线检测,曲线拟合适应弯道,跟踪则增强连续帧的车道线检测。

热门文章

最新文章