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

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

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


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


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

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

O]@{RC]_DR)3D_CT46@T7AV.png

NPCCTMM(}`[`YOLXL3BXW47.png

相关视频

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

JJU%GREROW(}MGH`66O)V$S.png

马尔科夫链蒙特卡洛方法


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 软件对学术能力测验建立层次(分层)贝叶斯模型

左右滑动查看更多

01

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

02

1RTB0I]0}WY6RWN9_FVC@UJ.png

03

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

04

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


相关文章
|
2月前
|
数据采集 机器学习/深度学习 数据可视化
R语言从数据到决策:R语言在商业分析中的实践
【9月更文挑战第1天】R语言在商业分析中的应用广泛而深入,从数据收集、预处理、分析到预测模型构建和决策支持,R语言都提供了强大的工具和功能。通过学习和掌握R语言在商业分析中的实践应用,我们可以更好地利用数据驱动企业决策,提升企业的竞争力和盈利能力。未来,随着大数据和人工智能技术的不断发展,R语言在商业分析领域的应用将更加广泛和深入,为企业带来更多的机遇和挑战。
|
25天前
|
数据挖掘 C语言 C++
R语言是一种强大的统计分析工具,提供了丰富的函数和包用于时间序列分析。
【10月更文挑战第21天】时间序列分析是一种重要的数据分析方法,广泛应用于经济学、金融学、气象学、生态学等领域。R语言是一种强大的统计分析工具,提供了丰富的函数和包用于时间序列分析。本文将介绍使用R语言进行时间序列分析的基本概念、方法和实例,帮助读者掌握R语言在时间序列分析中的应用。
41 3
|
6月前
|
数据可视化 数据挖掘 API
【R语言实战】聚类分析及可视化
【R语言实战】聚类分析及可视化
|
2月前
|
数据采集 数据可视化 数据挖掘
R语言在金融数据分析中的深度应用:探索数据背后的市场智慧
【9月更文挑战第1天】R语言在金融数据分析中展现出了强大的功能和广泛的应用前景。通过丰富的数据处理函数、强大的统计分析功能和优秀的可视化效果,R语言能够帮助金融机构深入挖掘数据价值,洞察市场动态。未来,随着金融数据的不断积累和技术的不断进步,R语言在金融数据分析中的应用将更加广泛和深入。
|
3月前
|
机器学习/深度学习 数据采集 数据可视化
R语言在数据科学中的应用实例:探索与预测分析
【8月更文挑战第31天】通过上述实例,我们展示了R语言在数据科学中的强大应用。从数据准备、探索、预处理到建模与预测,R语言提供了完整的解决方案和丰富的工具集。当然,数据科学远不止于此,随着技术的不断发展和业务需求的不断变化,我们需要不断学习和探索新的方法和工具,以更好地应对挑战,挖掘数据的潜在价值。 未来,随着大数据和人工智能技术的普及,R语言在数据科学领域的应用将更加广泛和深入。我们期待看到更多创新的应用实例,为各行各业的发展注入新的动力。
|
3月前
|
数据采集 存储 数据可视化
R语言时间序列分析:处理与建模时间序列数据的深度探索
【8月更文挑战第31天】R语言作为一款功能强大的数据分析工具,为处理时间序列数据提供了丰富的函数和包。从数据读取、预处理、建模到可视化,R语言都提供了灵活且强大的解决方案。然而,时间序列数据的处理和分析是一个复杂的过程,需要结合具体的应用场景和需求来选择合适的方法和模型。希望本文能为读者在R语言中进行时间序列分析提供一些有益的参考和启示。
|
3月前
|
资源调度 数据挖掘
R语言回归分析:线性回归模型的构建与评估
【8月更文挑战第31天】线性回归模型是统计分析中一种重要且实用的工具,能够帮助我们理解和预测自变量与因变量之间的线性关系。在R语言中,我们可以轻松地构建和评估线性回归模型,从而对数据背后的关系进行深入的探索和分析。
|
3月前
|
机器学习/深度学习 数据采集
R语言逻辑回归、GAM、LDA、KNN、PCA主成分分类分析预测房价及交叉验证
上述介绍仅为简要概述,每个模型在实施时都需要仔细调整与优化。为了实现高度精确的预测,模型选择与调参是至关重要的步骤,并且交叉验证是提升模型稳健性的有效途径。在真实世界的房价预测问题中,可能还需要结合地域经济、市场趋势等宏观因素进行综合分析。
69 3
|
6月前
|
数据可视化 算法
【R语言实战】——kNN和朴素贝叶斯方法实战
【R语言实战】——kNN和朴素贝叶斯方法实战
|
6月前
|
算法 项目管理
R语言实现蒙特卡洛模拟算法
R语言实现蒙特卡洛模拟算法

热门文章

最新文章