R语言BUGS/JAGS贝叶斯分析: 马尔科夫链蒙特卡洛方法（MCMC）采样（1）

2024-04-28 50

版权

本文内容由阿里云实名注册用户自发贡献，版权归原作者所有，阿里云开发者社区不拥有其著作权，亦不承担相应法律责任。具体规则请查看《阿里云开发者社区用户服务协议》和《阿里云开发者社区知识产权保护指引》。如果您发现本社区中有涉嫌抄袭的内容，填写侵权投诉表单进行举报，一经查实，本社区将立刻删除涉嫌侵权内容。

简介： R语言BUGS/JAGS贝叶斯分析: 马尔科夫链蒙特卡洛方法（MCMC）采样（1）

全文链接：http://tecdat.cn/?p=17884

在许多情况下，我们没有足够的计算能力评估空间中所有n维像素的后验概率。在这些情况下，我们倾向于利用称为Markov-Chain Monte Carlo 算法的程序。此方法使用参数空间中的随机跳跃来（最终）确定后验分布（点击文末“阅读原文”获取完整代码数据）。

相关视频

马尔可夫链蒙特卡罗方法MCMC原理与R语言实现

马尔科夫链蒙特卡洛方法

MCMC的关键如下：

跳跃概率的比例与后验概率的比例成正比。

跳跃概率可以表征为：

概率（跳跃）*概率（接受）

从长远来看，该链将花费大量时间在参数空间的高概率部分，从而实质上捕获了后验分布。有了足够的跳跃，长期分布将与联合后验概率分布匹配。

MCMC本质上是一种特殊类型的随机数生成器，旨在从难以描述（例如，多元，分层）的概率分布中采样。在许多/大多数情况下，后验分布是很难描述的概率分布。

MCMC使您可以从实际上不可能完全定义的概率分布中进行采样！

令人惊讶的是，MCMC的核心并不难于描述或实施。让我们看一个简单的MCMC算法。

Metropolis-Hastings算法

该算法与模拟退火算法非常相似。

MH算法可以表示为：

Prob（acceptB | A）= min（1，Posterior（B）Posterior（A）⋅Prob（b→a）Prob（a→b））

请注意，从本质上讲，这与“ Metropolis”模拟退火算法相同，后验概率代替了概率，并且 k 参数设置为1。

二元正态例子

请记住，MCMC采样器只是随机数生成器的一种。我们可以使用Metropolis-Hastings采样器来开发自己的随机数生成器，生成进行简单的已知分布。在此示例中，我们使用MH采样器从标准双变量正态概率分布生成随机数。

对于这个简单的示例，我们不需要MCMC采样器。一种实现方法是使用以下代码，该代码从具有相关参数ρ的双变量标准正态分布中绘制并可视化任意数量的独立样本。

#################
＃MCMC采样的简单示例
#################
#########
# ＃首先，让我们构建一个从双变量标准正态分布生成随机数的函数
rbvn<-function (n, rho)   #用于从二元标准正态分布中提取任意数量的独立样本。
{
        x <- rnorm(n, 0, 1)
        y <- rnorm(n, rho * x, sqrt(1 - rho^2))
        cbind(x, y)
}
#########
# 现在，从该分布图中绘制随机抽样 
bvn<-rbvn(10000,0.98)
par(mfrow=c(3,2))
plot(bvn,col=1:10000

T1%F11QG~A8A_50UWGIC_98.png

###############
# ＃Metropolis-Hastings双变量正态采样器的实现...
library(mvtnorm)    # 加载一个包，该包使我们能够计算mv正态分布的概率密度
metropoli<- function (n, rho=0.98){    # 双变量随机数生成器的MCMC采样器实现
    mat <- matrix(ncol = 2, nrow = n)   # 用于存储随机样本的矩阵
    x <- 0   # 所有参数的初始值
    prev <- dmvnorm(c(x,y),mean=c(0,0),sig
 # 起始位置分布的概率密度
    mat[1, ] <- c(x, y)        # 初始化马尔可夫链
      newx <- rnorm(1,x,0.5)     # 进行跳转 
      newprob <- dmvnorm(c(newx,newy),sigma =
  # 评估跳转
      ratio <- newprob/prev   # 计算旧位置（跳出）和建议位置（跳到）的概率之比。
      prob.accept <- min(1,ratio)     # 决定接受新跳跃的概率！
      if(rand<=prob.accept){
        x=newx;y=newy    # 将x和y设置为新位置
        mat[counter,] <- c(x,y)    # 将其存储在存储阵列中
        prev <- newprob    # 准备下一次迭代

然后，我们可以使用MH采样器从该已知分布中获取随机样本…

###########
# 测试新的M-H采样器
bvn<-metropoli(10000,0.98)
par(mfrow=c(3,2))
plot(bvn,col=1:10000)
plot(bvn,type=

W$N_}20YGEV1K@(~XE53E}V.png

点击标题查阅往期内容

$72_Q@UEAB2W}$[3JHW{X$~Y.png$

R语言用WinBUGS 软件对学术能力测验建立层次（分层）贝叶斯模型

左右滑动查看更多

V}2H@)@MW7JOS@53[1PB5(6.png

1RTB0I]0}WY6RWN9_FVC@UJ.png

4[]O`4V}N%MUCKPIEBUY`H5.png

$ML{%`P5MBC1_){VL36FCQO8.png$

让我们尝试解决一个问题。

MCMC对粘液瘤病进行调查

############
#黏液病示例的MCMC实现
############

subset(MyxDat,grade==1

##   grade day titer
## 1     1   2 5.207
## 2     1   2 5.734
## 3     1   2 6.613
## 4     1   3 5.997
## 5     1   3 6.612
## 6     1   3 6.810

选择使用Gamma分布。这是经验分布：

###########
# 第100次可视化粘液病数据
hist(Myx$titer,freq=FALSE)

3I9)VJF(]C6YJO@CD6QY4GG.png

Z}DP`4V7T0D%S()DZ($HLKR.png

我们需要估算最适合此经验分布的伽马速率和形状参数。这是适合此分布的Gamma的一个示例：

#########
# ...覆盖生成模型的数据（伽玛分布）
curve(dgamma(x,shape=40,scale=0.15),add=T,col="red")

$)Z{QMT3YB%Q(O9T]7V{X_CY.png$

二维（对数）似然面：

##############
# 定义二维参数空间
##############
shapevec <- seq(3,100,by=0.1)   
scalevec <- seq(0.01,0.5,by=0.001)
##############
# ＃定义参数空间内此网格上的似然面
##############
GammaLogLikelihoodFunction <- function(par
}
surface2D <- matrix(nrow=length(shapevec),ncol=length(scalevec))   #初始化存储变量
newparams <- c(sha
    surface2D[i,j] <- GammaLogLikelihoodFunction(newparams) 
############
# 可视化似然面
############
contour(x=shapevec,y=scalevec,z=surface2D,levels=c(-30,-40,-80,-500),add=T)

HD7[W_H_C88H@ESU6A3_`RB.png

Z`RJCP50OE@1E%KY~4WXFZX.png

这是MH算法的实现，用于找到后验分布！

首先，我们需要一个似然函数–这次，我们将返回真实概率–而不是对数转换的概率

############
#编写非对数转换的似然函数
GammaLike- function(params){
  prod(dgamma(Myx$titer,shape=params['shape']
params <- c(shape=40,

## shape scale 
## 40.00  0.15

GammaLike(params)

## [1] 2.906766e-22

GammaLogLike(params)

## [1] -49.58983

然后，我们需要预先分配参数！在这种情况下，我们分配gamma（shape = 0.01，scale = 100）和gamma（shape = 0.1，scale = 10）的分布（均值为1且方差略低）：

#############
# 函数返回参数空间中任意点的先验概率密度
GammaPriorFunction <- function(params){
  prior <- c(shape=NA,scale=NA)
],3,100)        
  # prior['scale'] <- dunif(params['
GammaLogPriorFunction <- function(params){
  prior <- c(shape=NA,scale=NA)
'],shape=0.001,scale=1000,log=T)
  # prior['shape'] <- dunif(params['shape'],3,100)        
  # prior['scale'] <- dunif(params['
curve(dgamma(x,shape=0.01,scale=1000),3,100)

DH$F90U%W$WFJ$8}@ITNA%O.png

params

## shape scale 
## 40.00  0.15

GammaPrior(params)

## [1] 1.104038e-06

prior2D <- matrix(nrow=length(shapevec),ncol=length(scalevec))   # 初始化存储变量
newparams <- c(shape=50,scale=0 
  for(j in 1:length(scalevec)){
    newparams['scale'] <- sca 
############
# 可视化似然面
############
image(x=shapevec,y=scalevec,z=prior2D,zlim

LP}BTT1](OTZ}(0A1RG589A.png

contour(x=shapevec,y=scalevec,z=prior2D,levels=c(-30,-40,-80,-500),add=T)

我们假设形状和比例在先验中是独立的（联合先验的乘法概率）。然而，并没有对后验参数相关性提出相同的假设，因为概率可以反映在后验分布中。

然后，我们需要一个函数，该函数可以计算参数空间中任何给定跳转的后验概率比率。因为我们正在处理后验概率的比率，所以 我们不需要计算归一化常数。

无需归一化常数，我们只需要计算加权似然比（即先验加权的似然比）

############
# 函数用于计算参数空间中任意两点之间的后验密度比
PosteriorRatio <- function(oldguess,newguess
  oldLik <- max(1e-90,GammaLikelihoodFunction(oldguess))   # 计算旧猜测的可能性和先验密度
  newLik <- GammaLikelihoodFunction(newguess)             # 在新的猜测值下计算可能性和先验密度
  return((newLik*newPrior)/(oldLik*oldPrior))          # 计算加权似然比
}
PosteriorRatio2 <- function(oldguess,newguess){
  oldLogLik <- GammaLogLikelihoodFunction(oldguess)   # 计算旧猜测的可能性和先验密度
  newLogLik <- GammaLogLikelihoodFunction(newguess)             # 在新的猜测值下计算可能性和先验密度
  return(exp((newLogLik+newLogPrior)-(oldLogLik+oldLogPrior)))          # 计算加权似然比
}

## [1] 0.01436301

PosteriorRatio2(oldguess,newguess)

## [1] 0.01436301

然后，我们需要一个函数进行新的推测或在参数空间中跳转：

############
# 为参数空间中的跳转定义：使用正态分布函数进行新的推测
newGuess <- function(oldguess)
  jump <- c(shape=rnorm(1,mean=0,sd=sdshapejump),scale=rnorm(1,0,sdscalejump))
  newguess <- abs(oldguess + ju
}
  # 在原始推测附近设置新的推测
newGuess(oldguess=params)

##      shape      scale 
## 35.7132110  0.1576337

newGuess(oldguess=params)

##      shape      scale 
## 45.1202345  0.2094243

newGuess(oldguess=params)

##       shape       scale 
## 42.87840436  0.08152061

现在，我们准备实现Metropolis-Hastings MCMC算法：

我们需要一个初始点：

##########
# 在参数spacer中设置起点
startingvals <- c(shape=75,scale=0.28)    # 算法的起点

测试函数

###########
# 尝试我们的新函数
newguess <- newGuess(startingvals)    # 在参数空间中跳跃
newguess

##      shape      scale 
## 73.9663949  0.3149796

PosteriorRatio2(startingvals,newguess)   # 后验比例差异

## [1] 2.922783e-57

现在让我们看一下Metropolis-Hastings：

###############
#可视化Metropolis-Hastings
chain.length <- 11
gth,ncol=2)
colnames(guesses) <- names(startingvals)
guesses[1,] <- startingvals
counter <- 2
  post.rat <- PosteriorRatio2(oldguess,newguess)
  prob.accept <- min(1,post
    oldguess <- newguess
    guesses[coun
#可视化
contour(x=shapevec,y=scal

$_$(Q]{3F1V0EM3K7QI5`OAP.png$

我们运行更长的时间

##########
# 获取更多MCMC示例
chain.length <- 100
oldgu
counter <- 2
while(counter <= chain.length){
  newguess <- newGuess(oldguess)
  post.rat <- Posterio
  rand <- runif(1)
  if(rand<=prob.accept){
ewguess 
    counter=counte
#可视化
image(x=shapevec,y=scalevec,z=su
urface2D,levels=c(-30,-40,-80,-5
lines(guesses,col="red")

$[N~)(L{GMPI108HA_6PA5XN.png$

更长的时间

############
#更长的时间
chain.length <- 1000
oldguess <- startingvals
chain.length,ncol=2)
colnames(guesses) <- names(startingvals)
guesses[1,] <- startingva
ess)
  post.rat <- PosteriorRatio2(oldguess,newguess)
  prob.accept <- min(1,post.rat)
  rand <- runif(1)
    guesse
#可视化
image(x=shapevec,y=scalevec,
rface2D,levels=c(-30,-40,-80,-500),a

XBU5FH8HT0F468IZU5))V7B.png

R语言BUGS/JAGS贝叶斯分析: 马尔科夫链蒙特卡洛方法（MCMC）采样（1）

相关视频：马尔可夫链原理可视化解释与R语言区制转换Markov regime switching实例

相关视频

马尔科夫链蒙特卡洛方法

Metropolis-Hastings算法

二元正态例子

MCMC对粘液瘤病进行调查

热门文章

最新文章

相关课程

相关电子书

相关实验场景

热门

活动广场

任务中心

开发者评测

高校计划

乘风者计划

训练营

阿里云MVP

话题

直播

下载

镜像站

技术资料

插件

R语言BUGS/JAGS贝叶斯分析: 马尔科夫链蒙特卡洛方法（MCMC）采样（1）

相关视频：马尔可夫链原理可视化解释与R语言区制转换Markov regime switching实例

相关视频

马尔科夫链蒙特卡洛方法

Metropolis-Hastings算法

二元正态例子

MCMC对粘液瘤病进行调查

热门文章

最新文章

相关课程

相关电子书

相关实验场景