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

最近咱们被客户要求撰写对于BUGS/JAGS贝叶斯剖析的钻研报告，包含一些图形和统计输入。

在许多状况下，咱们没有足够的计算能力评估空间中所有n维像素的后验概率 。在这些状况下，咱们偏向于利用称为Markov-Chain Monte Carlo 算法的程序 。此办法应用参数空间中的随机跳跃来（最终）确定后验散布

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

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

相干视频

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

，时长08:47

马尔科夫链蒙特卡洛办法

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

################ ＃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=

点击题目查阅往期内容

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

左右滑动查看更多

01

02

03

04

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

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)

咱们须要估算最适宜此教训散布的伽马速率和形态参数。这是适宜此散布的Gamma的一个示例：

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

二维（对数）似然面：

############### 定义二维参数空间##############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)

这是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)

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

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-Hastingschain.length <- 11gth,ncol=2)colnames(guesses) <- names(startingvals)guesses[1,] <- startingvalscounter <- 2  post.rat <- PosteriorRatio2(oldguess,newguess)  prob.accept <- min(1,post    oldguess <- newguess    guesses[coun#可视化contour(x=shapevec,y=scal

咱们运行更长的工夫 

########### 获取更多MCMC示例chain.length <- 100oldgucounter <- 2while(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=suurface2D,levels=c(-30,-40,-80,-5lines(guesses,col="red")

更长的工夫

#############更长的工夫chain.length <- 1000oldguess <- startingvalschain.length,ncol=2)colnames(guesses) <- names(startingvals)guesses[1,] <- startingvaess)  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

看起来更好！搜索算法能够很好地找到参数空间的高似然局部！

当初，让咱们看一下“ shape”参数的链

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

对于比例参数 

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

咱们能够说这些链曾经收敛于形态参数的后验散布吗？

首先，链的终点“记住”起始值，因而不是固定散布。咱们须要删除链的第一局部。

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

但这看起来还不是很好。让咱们运行更长的工夫，看看是否失去的货色看起来更像是随机数生成器（白噪声）

########### 再试一次-运行更长的工夫chain.length <- 20000oldguess <- startingvo2(oldguess,newguess)  prob.accept <- mi

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

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

当初，让咱们再次看一下链条

在评估这些迹线图时，咱们心愿看到看起来像白噪声的“安稳散布”。该轨迹图看起来可能具备一些自相干。解决此问题的一种办法是稠密MCMC样本：

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

当初咱们能够查看咱们的后验散布！

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

咱们能够像以前一样可视化。

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

能够批改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))

############## 当初结构一个吉布斯采样器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)

在这里，马尔可夫链的样本中有很多显著的自相干。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.001file="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)

评估收敛

第一步是视觉查看-咱们寻找以下内容来评估收敛性：

  • 当视为“轨迹图”时，每个参数的链应看起来像白噪声（含糊的毛毛虫）或相似的噪声
  • 多个具备不同起始条件的链条看起来应该雷同

咱们可能在这里能够做得更好的一种办法是使链条运行更长的工夫，并抛弃初始样本咱们还能够。

通过细化链来尝试缩小自相干，咱们每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)

从视觉上看，这看起来更好。当初咱们能够应用更多的定量收敛指标。

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诊断。如果测试失败，则应尝试运行更长的链！

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

点击文末 “浏览原文”

获取全文残缺代码数据资料。

本文选自《R语言BUGS/JAGS贝叶斯剖析: 马尔科夫链蒙特卡洛办法（MCMC）采样》。

点击题目查阅往期内容

R语言MCMC:Metropolis-Hastings采样用于回归的贝叶斯预计
Python用MCMC马尔科夫链蒙特卡洛、回绝抽样和Metropolis-Hastings采样算法
R语言贝叶斯METROPOLIS-HASTINGS GIBBS 吉布斯采样器预计变点指数分布剖析泊松过程车站等待时间
R语言马尔可夫MCMC中的METROPOLIS HASTINGS，MH算法抽样（采样）法可视化实例
python贝叶斯随机过程：马尔可夫链Markov-Chain，MC和Metropolis-Hastings，MH采样算法可视化
Python贝叶斯推断Metropolis-Hastings（M-H）MCMC采样算法的实现
Metropolis Hastings采样和贝叶斯泊松回归Poisson模型
Matlab用BUGS马尔可夫区制转换Markov switching随机稳定率模型、序列蒙特卡罗SMC、M H采样剖析工夫序列R语言RSTAN MCMC：NUTS采样算法用LASSO 构建贝叶斯线性回归模型剖析职业声望数据
R语言BUGS序列蒙特卡罗SMC、马尔可夫转换随机稳定率SV模型、粒子滤波、Metropolis Hasting采样工夫序列剖析
R语言Metropolis Hastings采样和贝叶斯泊松回归Poisson模型
R语言贝叶斯MCMC：用rstan建设线性回归模型剖析汽车数据和可视化诊断
R语言贝叶斯MCMC：GLM逻辑回归、Rstan线性回归、Metropolis Hastings与Gibbs采样算法实例
R语言贝叶斯Poisson泊松-正态分布模型剖析职业足球比赛进球数
R语言用Rcpp减速Metropolis-Hastings抽样预计贝叶斯逻辑回归模型的参数
R语言逻辑回归、Naive Bayes贝叶斯、决策树、随机森林算法预测心脏病
R语言中贝叶斯网络（BN）、动静贝叶斯网络、线性模型剖析错颌畸形数据
R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归
Python贝叶斯回归剖析住房累赘能力数据集
R语言实现贝叶斯分位数回归、lasso和自适应lasso贝叶斯分位数回归剖析
Python用PyMC3实现贝叶斯线性回归模型
R语言用WinBUGS 软件对学术能力测验建设档次（分层）贝叶斯模型
R语言Gibbs抽样的贝叶斯简略线性回归仿真剖析
R语言和STAN,JAGS：用RSTAN,RJAG建设贝叶斯多元线性回归预测选举数据
R语言基于copula的贝叶斯分层混合模型的诊断准确性钻研
R语言贝叶斯线性回归和多元线性回归构建工资预测模型
R语言贝叶斯推断与MCMC：实现Metropolis-Hastings 采样算法示例
R语言stan进行基于贝叶斯推断的回归模型
R语言中RStan贝叶斯层次模型剖析示例
R语言应用Metropolis-Hastings采样算法自适应贝叶斯预计与可视化
R语言随机搜寻变量抉择SSVS预计贝叶斯向量自回归（BVAR）模型
WinBUGS对多元随机稳定率模型：贝叶斯预计与模型比拟
R语言实现MCMC中的Metropolis–Hastings算法与吉布斯采样
R语言贝叶斯推断与MCMC：实现Metropolis-Hastings 采样算法示例
R语言应用Metropolis-Hastings采样算法自适应贝叶斯预计与可视化
视频：R语言中的Stan概率编程MCMC采样的贝叶斯模型
R语言MCMC:Metropolis-Hastings采样用于回归的贝叶斯预计