共计 14002 个字符，预计需要花费 36 分钟才能阅读完成。
全文链接：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.15GammaLike(params)## [1] 2.906766e-22GammaLogLike(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.15GammaPrior(params)## [1] 1.104038e-06prior2D <- 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,zlimcontour(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.01436301PosteriorRatio2(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.1576337newGuess(oldguess=params)## shape scale
## 45.1202345 0.2094243newGuess(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.3149796PosteriorRatio2(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咱们运行更长的工夫
##########
# 获取更多 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")更长的工夫
############
#更长的工夫
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看起来更好！搜索算法能够很好地找到参数空间的高似然局部！
当初，让咱们看一下“shape”参数的链
#############
# 评估 MCMC 样本的“轨迹图”...
##### Shape 参数
plot(1:chain.length,guesses[,'sha对于比例参数
###### 比例参数
plot(1:chain.length,guesses[,'scale'],type="l咱们能够说这些链曾经收敛于形态参数的后验散布吗？
首先，链的终点“记住”起始值，因而不是固定散布。咱们须要删除链的第一局部。
############
# 删除预烧期（容许 MCMC 有一段时间达到后验概率）burn.in <- 100
MCMCsamples <- guesses[-c(1:burn.in),]但这看起来还不是很好。让咱们运行更长的工夫，看看是否失去的货色看起来更像是随机数生成器（白噪声）
##########
# 再试一次 - 运行更长的工夫
chain.length <- 20000
oldguess <- startingv
o2(oldguess,newguess)
prob.accept <- mi让咱们首先删除前 5000 个样本作为预烧期
#############
# 应用更长的“预烧”burn.in <- 5000
MCMCsamples <- 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.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.2202549init.vals.for.bugs()## $shape
## [1] 89.85618
##
## $scale
## [1] 0.2678733init.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':
##
## traceplotjags.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 simulationjags.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.2389plot(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 simulationjagsfit.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.2773plot(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 采样用于回归的贝叶斯预计
