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

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

看起来更好!搜索算法可以很好地找到参数空间的高似然部分!

现在,让我们看一下“ shape”参数的链

#############
# 评估MCMC样本的“轨迹图” ...
##### Shape 参数
plot(1:chain.length,guesses[,'sha

W_81)1IA$GZF@~{S6F1}$WY.png

对于比例参数

###### 比例参数 
plot(1:chain.length,guesses[,'scale'],type="l

`N@E(`XU{G]JOTPGD`RUZ30.png

我们可以说这些链已经收敛于形状参数的后验分布吗?

首先,链的起点“记住”起始值,因此不是固定分布。我们需要删除链的第一部分。

############
# 删除预烧期(允许MCMC有一段时间达到后验概率)
burn.in <- 100
MCMCsamples <- guesses[-c(1:burn.in),]

0J6}25P$)Z_D}(R0X20~5AG.png

{L@WNM]%R[`O`@MOG[_HM9V.png

但这看起来还不是很好。让我们运行更长的时间,看看是否得到的东西看起来更像是随机数生成器(白噪声)

##########
# 再试一次-运行更长的时间
chain.length <- 20000
oldguess <- startingv
o2(oldguess,newguess)
  prob.accept <- mi

YZT}9K7G4E3S[6%96M[W597.png

让我们首先删除前5000个样本作为预烧期

#############
# 使用更长的“预烧”
burn.in <- 5000
MCMCsamples <- guesses[-c(1:bur

现在,让我们再次看一下链条

%9WPFZZ`~P(1FRPW`TS(0F3.png

6CO)FLDRZC4WO1KR$Q9W32Y.png

在评估这些迹线图时,我们希望看到看起来像白噪声的“平稳分布”。该轨迹图看起来可能具有一些自相关。解决此问题的一种方法是稀疏MCMC样本:

##########
# “稀疏” MCMC样本
thinnedMCMC <- MCMCsamples[seq(1,chain.length,by=5),]

H]NAZKF3X(MXM~F28%08U26.png

XS`VY$3CO_O(N84C[Y9S%KL.png

现在我们可以检查我们的后验分布!

# 可视化后验分布
plot(density(thinnedMCMC[,'scale'])

OIQ[RPC]N6IP}_LPIQK{0$E.png

VV~F8$44DG4~VE4P[M0Q{)I.png

我们可以像以前一样可视化。

#########
# 更多后验概率检察
par(mfrow=c(3,2))
plot(thinnedMCMC,col=1:10000)
plot(thinnedMCMC,type="l")

(QDV7`F0W_O_A{SN~)}~)TQ.png

可以修改Metropolis-Hastings MCMC方法来拟合任意模型的任意数量的自由参数。但是,MH算法本身不一定是最有效和灵活的。在实验中,我们使用吉布斯采样,大多采用建模语言 BUGS

注意:BUGS实现(例如JAGS)实际上倾向于结合使用MH和Gibbs采样,MH和Gibbs采样器并不是唯一的MCMC例程。例如,“ stan”使用MH采样的一种改进形式,称为“ Hamiltonian Monte Carlo”。

吉布斯Gibbs采样器

Gibbs采样器非常简单有效。基本上,该算法从完整的条件 概率分布(即, 在模型中所有其他参数的已知值作为条件的条件下,对任意参数i的后验分布)中进行 连续采样 。

在很多情况下,我们不能直接制定出我们的模型后验分布,但我们 可以 分析出条件后验分布。尽管如此,即使它在分析上不易处理,我们也可以使用单变量MH程序作为最后方法。

:为什么Gibbs采样器通常比纯MH采样器效率更高?

二元正态例子

MCMC采样器只是随机数生成器的一种。我们可以使用Gibbs采样器来开发自己的随机数生成器,以实现相当简单的已知分布。在此示例中,我们使用Gibbs采样器从标准双变量正态概率分布生成随机数。注意,吉布斯采样器在许多方面都比MH算法更简单明了。

#############
#Gibbs采样器的简单示例
#############
########
# 首先,回顾一下我们简单的双变量正态采样器
rbvn<-function (n, rho){  #f函数用于从双变量标准正态分布中提取任意数量的独立样本。
        x <- rnorm(n, 0, 1)
        y <- rnorm(n, rho * x, sqrt(1 - rho^2))

AML1N@CHA`[4OWSNTEP8XDT.png

#############
# 现在构造一个吉布斯采样器
gibbs<-function (n, rho){    # 双变量随机数生成器的gibbs采样器实现
    mat <- matrix(ncol = 2, nrow = n)   # 用于存储随机样本的矩阵
    mat[1, ] <- c(x, y)        # initialize the markov chain
    for (i in 2:n) {
            x <- rnorm(1, rho * y, sqrt(1 - rho^2))        # 以y为条件的x中的样本
            y <- rnorm(1, rho * x, sqrt(1 - rho^2))        # 以x为条件的y中的样本
            mat[i, ] <- c(x, y)

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

##########
# 测试吉布斯采样器
plot(ts(bvn[,2]))
hist(bvn[,1],40)
hist(bvn[,2],40)

$%TF_MO1~(YZK]U1[XBU[Y6.png

在这里,马尔可夫链的样本中有很多明显的自相关。Gibbs采样器经常有此问题。

示例

BUGS语言

最后,让我们为我们最喜欢的粘瘤病示例创建一个Gibbs采样器,为此,我们将使用BUGS语言(在JAGS中实现)来帮助我们!

JAGS,全称是Just another Gibbs sampler,是基于BUGS语言开发的利用MCMC来进行贝叶斯建模的软件包。它没有提供建模所用的GUI以及MCMC抽样的后处理,这些要在其它的程序软件上来处理,比如说利用R包(rjags)来调用JAGS并后处理MCMC的输出。JAGS相对于WinBUGS/OpenBUGS的主要优点在于平台的独立性,可以应用于各种操作系统,而WinBUGS/OpenBUGS只能应用于windows系统;JAGS也可以在64-bit平台上以64-bit应用来进行编译。

BUGS语言看起来与R类似,但是有几个主要区别:

  • 首先,BUGS是一种编译语言,因此代码中的操作顺序并不重要
  • BUGS不是矢量化的-您需要使用FOR循环
  • 在BUGS中,几个概率分布的参数差异很大。值得注意的是,正态分布通过均值和精度(1 / Variance )进行参数化。

这是用BUGS语言实现的粘液病示例:

model {
  #############
  # 似然
  ############
  for(obs in 1:n.observations){
    titer[obs] ~ dgamma(shape,rate
  #############
  # 先验
  ############
  rate <- 1/scale   # 将BUGS的scale参数转换为“ rate”
}

我们可以使用R中的“ cat”函数将此模型写到您的工作目录中的文本文件中:

###########
# BUGS建模语言中的粘液瘤示例
##########
# 将BUGS模型写入文件
cat("
  model {
    #############
    # 似然
    ############
    for(obs in 1:n.observations){
      titer[obs] ~ dgamma(shape,rate)
    #############
    # 先验
    ############
    shape ~ dgamma(0.001,0.001
file="BUGSmodel.txt")

现在我们已经将BUGS模型打包为文本文件,我们将数据捆绑到一个列表对象中,该列表对象包含BUGS代码中引用的所有相关数据:

############
# 将数据封装到单个“列表”对象中
myx.data <- list(
  n.observations = length(Myx$titer
## $titer
##  [1] 5.207 5.734 6.613 5.997 6.612 6.810 5.930 6.501 7.182 7.292 7.819
## [12] 7.489 6.918 6.808 6.235 6.916 4.196 7.682 8.189 7.707 7.597 7.112
## [23] 7.354 7.158 7.466 7.927 8.499
## 
## $n.observations
## [1] 27

然后,我们需要为所有参数定义初始值。将其定义为一个函数很方便,因此可以使用不同的起始值来初始化每个MCMC链。

###########
# 用于为所有自由参数生成随机初始值的函数
init.vals.for.bugs()
## $shape
## [1] 51.80414
## 
## $scale
## [1] 0.2202549
init.vals.for.bugs()
## $shape
## [1] 89.85618
## 
## $scale
## [1] 0.2678733
init.vals.for.bugs()
## $shape
## [1] 69.22457
## 
## $scale
## [1] 0.1728908

现在我们可以调用JAGS!

###########
# 运行 JAGS 
##########
## Loading required package: rjags
## The following object is masked from 'package:coda':
## 
##     traceplot
jags.fit <- run.jags(model="BUGSmodel.txt",
## Compiling rjags model...
## Calling the simulation using the rjags method...
## Adapting the model for 100 iterations...
## Running the model for 5000 iterations...
## Simulation complete
## Finished running the simulation
jags.fit)   # 转换为“ MCMC”对象(CODA包)
## 
## Iterations = 101:5100
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##          Mean       SD  Naive SE Time-series SE
## shape 51.6013 14.36711 0.1173070       2.216657
## scale  0.1454  0.04405 0.0003597       0.006722
## 
## 2. Quantiles for each variable:
## 
##           2.5%     25%     50%     75%   97.5%
## shape 28.76333 40.9574 50.1722 60.2463 82.7532
## scale  0.08346  0.1148  0.1378  0.1687  0.2389
plot(jagsfit.mcmc)


L~J2U1A6NS%2)47TJ%1M%OX.png

ZX6}L50$35E4HXJ[7QMT{O2.png

评估收敛

第一步是视觉检查-我们寻找以下内容来评估收敛性:

  • 当视为“轨迹图”时,每个参数的链应看起来像白噪声(模糊的毛毛虫)或类似的噪声
  • 多个具有不同起始条件的链条看起来应该相同

我们可能在这里可以做得更好的一种方法是使链条运行更长的时间,并丢弃初始样本我们还可以。

通过细化链来尝试减少自相关,我们每20个样本中仅保留1个。

################
#运行链更长时间
jags.fit <- 
 sample = 10000,n.chains = 3,adapt = 1000,burnin = 1000,
                     summarise = FALSE,thin=20,method = "parallel" )
## Calling 3 simulations using the parallel method...
## Following the progress of chain 1 (the program will wait for all
## chains to finish before continuing):
## Welcome to JAGS 4.2.0 on Wed Oct 23 11:33:35 2019
## JAGS is free software and comes with ABSOLUTELY NO WARRANTY
## Loading module: basemod: ok
## Loading module: bugs: ok
## . . Reading data file data.txt
## . Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 27
##    Unobserved stochastic nodes: 2
##    Total graph size: 37
## . Reading parameter file inits1.txt
## . Initializing model
## . Adapting 1000
## -------------------------------------------------| 1000
## ++++++++++++++++++++++++++++++++++++++++++++++++++ 100%
## Adaptation successful
## . Updating 1000
## -------------------------------------------------| 1000
## ************************************************** 100%
## . . . Updating 200000
## -------------------------------------------------| 200000
## ************************************************** 100%
## . . . . Updating 0
## . Deleting model
## . 
## All chains have finished
## Simulation complete.  Reading coda files...
## Coda files loaded successfully
## Finished running the simulation
jagsfit.mcmc <- as.mcmc.list
  # 转换为“ MCMC”对象(CODA包)
## 
## Iterations = 2001:201981
## Thinning interval = 20 
## Number of chains = 3 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##          Mean      SD  Naive SE Time-series SE
## shape 47.1460 12.9801 0.0749404       0.292218
## scale  0.1591  0.0482 0.0002783       0.001075
## 
## 2. Quantiles for each variable:
## 
##           2.5%     25%     50%     75%   97.5%
## shape 25.14757 37.9988 45.9142 55.1436 75.5850
## scale  0.09147  0.1256  0.1509  0.1825  0.2773
plot(jagsfit.mcmc)

L~J2U1A6NS%2)47TJ%1M%OX.png

60F8T01Z6YO53(U_@AX1_$6.png

从视觉上看,这看起来更好。现在我们可以使用更多的定量收敛指标。

Gelman-Rubin诊断

一种简单而直观的收敛诊断程序是 Gelman-Rubin诊断程序,该诊断程序基于简单的蒙特卡洛误差来评估链之间是否比应有的链距更大:

##############
# 运行收敛诊断
gelman(jagsfit.mcmc)
## Potential scale reduction factors:
## 
##       Point est. Upper C.I.
## shape          1       1.00
## scale          1       1.01
## 
## Multivariate psrf
## 
## 1

通常,1.1或更高的值被认为收敛不佳。为模型中的所有可用参数计算GR诊断。如果测试失败,则应尝试运行更长的链!

所以这个模型看起来不错!

相关文章
|
2月前
|
数据采集 机器学习/深度学习 数据可视化
R语言从数据到决策:R语言在商业分析中的实践
【9月更文挑战第1天】R语言在商业分析中的应用广泛而深入,从数据收集、预处理、分析到预测模型构建和决策支持,R语言都提供了强大的工具和功能。通过学习和掌握R语言在商业分析中的实践应用,我们可以更好地利用数据驱动企业决策,提升企业的竞争力和盈利能力。未来,随着大数据和人工智能技术的不断发展,R语言在商业分析领域的应用将更加广泛和深入,为企业带来更多的机遇和挑战。
|
25天前
|
数据挖掘 C语言 C++
R语言是一种强大的统计分析工具,提供了丰富的函数和包用于时间序列分析。
【10月更文挑战第21天】时间序列分析是一种重要的数据分析方法,广泛应用于经济学、金融学、气象学、生态学等领域。R语言是一种强大的统计分析工具,提供了丰富的函数和包用于时间序列分析。本文将介绍使用R语言进行时间序列分析的基本概念、方法和实例,帮助读者掌握R语言在时间序列分析中的应用。
42 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语言实现蒙特卡洛模拟算法