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

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

时间:2023-12-13 23:35:18浏览次数:52  
标签:Metropolis 后验 算法 MCMC 提议 我们 Hastings

全文链接:https://tecdat.cn/?p=34543

原文出处:拓端数据部落公众号

贝叶斯MCMC模拟是一个丰富的领域,涵盖了各种算法,共同目标是近似后验模型。例如,使用的rstan包采用了一个Hamiltonian Monte Carlo算法。用于贝叶斯建模的另一个rjags包采用了Gibbs sampling算法。尽管细节有所不同,但这两种算法都是基于基本的Metropolis-Hastings算法的变体。

主要思想

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

image.png

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

image.png

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

image.png

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

image.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")

image.png

Metropolis-Hastings算法

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

   
current <- 3

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

   
set.seed(8)
proposal

image.png

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

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

image.png

image.png

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

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

image.png

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

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

image.png

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

我们首先指定one_mh_iteration是两个参数的function():均匀分布的半宽度w和当前链值current

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

  r
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()可以复制上面的结果:

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

image.png

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

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

image.png

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

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

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

image.png

实施Metropolis-Hastings

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

  r
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的均匀提议模型:

  r
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")

QQ截图20231128144503.png

Beta-Binomial模型例子

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

image.png

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

image.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()来确定每个停止位置:

  R
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后验分布遍历

  R
set.seed(84735)
bebn_im <- betb_tour(N = 5000, a = 1, b = 1)

# 绘制结果
ggplot(beta

总结

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

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

deep-learning-in-product-design-533765077-standard-1536x1536.jpg

标签:Metropolis,后验,算法,MCMC,提议,我们,Hastings
From: https://www.cnblogs.com/tecdat/p/17900189.html

相关文章

  • R语言贝叶斯Metropolis-Hastings Gibbs 吉布斯采样器估计变点指数分布分析泊松过程车
    原文链接:http://tecdat.cn/?p=26578 原文出处:拓端数据部落公众号最近我们被客户要求撰写关于吉布斯采样器的研究报告,包括一些图形和统计输出。指数分布是泊松过程中事件之间时间的概率分布,因此它用于预测到下一个事件的等待时间,例如,您需要在公共汽车站等待的时间,直到下一班车......
  • R语言贝叶斯Metropolis-Hastings Gibbs 吉布斯采样器估计变点指数分布分析泊松过程车
    最近我们被客户要求撰写关于吉布斯采样器的研究报告,包括一些图形和统计输出。指数分布是泊松过程中事件之间时间的概率分布,因此它用于预测到下一个事件的等待时间,例如,您需要在公共汽车站等待的时间,直到下一班车到了。在本文中,我们将使用指数分布,假设它的参数λ,即事件之间的平均......
  • python机器学习算法原理实现——MCMC算法之gibbs采样
    【算法原理】Gibbs采样是一种用于估计多元分布的联合概率分布的方法。在MCNC(Markov Chain Monte Carlo)中,Gibbs采样是一种常用的方法。通俗理解Gibbs采样,可以想象你在一个多维空间中,你需要找到这个空间的某个特定区域(这个区域代表了你感兴趣的分布)。但是,你不能直接看到整个空间,只......
  • matlab用Logistic逻辑回归建模和马尔可夫链蒙特卡罗MCMC方法分析汽车实验数据
    原文链接:http://tecdat.cn/?p=24103原文出处:拓端数据部落公众号 此示例说明如何使用逻辑回归模型进行贝叶斯推断。统计推断通常基于最大似然估计(MLE)。MLE选择能够使数据似然最大化的参数,是一种较为自然的方法。在MLE中,假定参数是未知但固定的数值,并在一定的置信度下进......
  • R语言实现MCMC中的Metropolis–Hastings算法与吉布斯采样|附代码数据
     最近我们被客户要求撰写关于MCMC的研究报告,包括一些图形和统计输出。创建测试数据第一步,我们创建一些测试数据,用来拟合我们的模型。我们假设预测变量和因变量之间存在线性关系,所以我们用线性模型并添加一些噪音。  trueA<-5trueB<-0trueSd<-10sampleSize<-31......
  • R语言贝叶斯Metropolis-Hastings Gibbs 吉布斯采样器估计变点指数分布分析泊松过程车
    原文链接:http://tecdat.cn/?p=26578 原文出处:拓端数据部落公众号最近我们被客户要求撰写关于吉布斯采样器的研究报告,包括一些图形和统计输出。指数分布是泊松过程中事件之间时间的概率分布,因此它用于预测到下一个事件的等待时间,例如,您需要在公共汽车站等待的时间,直到下一班车......
  • Python贝叶斯MCMC:Metropolis-Hastings、Gibbs抽样、分层模型、收敛性评估
     全文链接:https://tecdat.cn/?p=33961原文出处:拓端数据部落公众号在常规的马尔可夫链模型中,我们通常感兴趣的是找到一个平衡分布。MCMC则是反过来思考——我们将平衡分布固定为后验分布:并寻找一种转移核,使其收敛到该平衡分布。岛屿示例首先提供一个示例,以具体展示Metropo......
  • Metropolis Algorithms for Representative Subgraph Sampling
    目录概主要内容MetropolisgraphsamplingH\¨{u}blerC.andKriegelH.,BorgwardtK.andGhahramaniZ.Metropolisalgorithmsforrepresentativesubgraphsampling.ICDM,2008.概提出了一种尽可能保持拓扑结构的子图采样方法.主要内容假设我们有一个大图\(G\),......
  • Python随机波动性SV模型:贝叶斯推断马尔可夫链蒙特卡洛MCMC分析英镑/美元汇率时间序列
    全文链接:https://tecdat.cn/?p=33885原文出处:拓端数据部落公众号本文描述了帮助客户使用马尔可夫链蒙特卡洛(MCMC)方法通过贝叶斯方法估计基本的单变量随机波动模型,就像Kim等人(1998年)所做的那样。定义模型以及从条件后验中抽取样本的函数的代码也在Python脚本中提供。  ......
  • R语言中的Stan概率编程MCMC采样的贝叶斯模型|附代码数据
    原文链接:http://tecdat.cn/?p=11161最近我们被客户要求撰写关于贝叶斯模型的研究报告,包括一些图形和统计输出。概率编程使我们能够实现统计模型,而不必担心技术细节。这对于基于MCMC采样的贝叶斯模型特别有用R语言中RStan贝叶斯层次模型分析示例stan简介Stan是用于贝叶斯推理......