这篇文章首先以一个最简单的例子+code，带你体验什么是MCMC，随后会推荐一些实用资源。[不知怎么的，这篇文章本来要写R+JAGS教程，最后硬是写成了一篇杂文😓，不过也还算是一篇过得去的低阶入门文章吧。]
贝叶斯模型的后验分布可以使用共轭先验进行推导，但是更为复杂的问题无法使用closed formed的形式找到后验分布，因为The Devil is in the Denominator! 这时可以使用基于MCMC采样的方法获得后验分布的估计--BUSGS(Bayesian Analysis Using Gibbs Sampling)。另外我了解到还有一种不需要采样（pakage: r-inla）的方法，但不适合新手并且适用的问题有限。
WinBugs，OpenBUGS和JAGS是实现BUGS的三个软件。OpenBUGS和WinBUGS在Windows操作系统上运行，有图形界面。JAGS(Just another Gibbs sampler)可在win/macos/linux上使用，无图形界面。几个软件使用了基本相同的语法，都可以用R调用。
一个最简单的栗子
假设观测到了1000个值来自正态分布的值，但是并不知道其总体的信息，目标是使用BUGS找出总体参数的后验分布。
1
安装所需软件
安装JAGS到电脑
R中安装：
rjags
coda
MCMCvis
2
生成数据
set.seed(432104)
n <- 1000
x <- rnorm(n, 0, 5)
3
建立JAGS模型
因为不知道总体的信息，我们假设观测值来自于均值为mu，presicion为tau的正态总体。其中tau=1/sigma^2。对于这两个未知的参数，我们的先验信息有限，如果觉得mu大概率为0的话，可以假设mu采样于正态N(0,0.0001)， sigma采样于均匀分布Unif(0,100)。[可以尝试使用其他的先验看一下对结果有什么影响]
JAGS模型如下
model {
for (i in 1:N){
x[i] ~ dnorm(mu, tau)
}
mu ~ dnorm(0,.0001)
tau <- pow(sigma, -2)
sigma ~ dunif(0,100)
}一般情况下，我们可以将JAGS的模型单独保存为一个文件，使用R调用JAGS调用JAGS运行该文件。这里可以偷懒使用textConnection，将这个模型保存为成一个string，假装我们已经已经保存了一个模型的文件。
model1.string <-"
model {
for (i in 1:N){
x[i] ~ dnorm(mu, tau)
}
mu ~ dnorm(0,.0001)
tau <- pow(sigma, -2)
sigma ~ dunif(0,100)
}"
model1.spec<-textConnection(model1.string)调用JAGS建立模型并采样
jags <- jags.model(model1.spec,
data = list('x' = x,
'N' = n),
n.chains=4,
n.adapt=100)
model=coda.samples(jags,
c('mu', 'tau','sigma'),
1000)➡jags.model
1. JAGS模型文件
这里填入之前的模型model1.spec
如果你将模型保存为了文件，这里则需改为该文件的路径+文件名。
2. data
list中变量的名称应该和JAGS所需要的数据名称一致。
此外需要注意数据是否有NA，数据的格式不能是charater
3. inits
每一个链的初始值，默认会随机生成。但是推荐手动设置，因为你在设置该参数的时候会想一下在估计什么参数。其形式是list of list
比如2个chain的初始值可以这么定义：
inits=list(
list(mu=0, sigma=6),
list(mu=10, sigma=50))4. n.chains
生成多少个chain
5. n.adapt
在一个chain中获得最优的sampler的迭代次数
➡coda.samples
这里用了coda.samples没用jags.samples是为了方便使用之后的可视化工具MCMCvis，使用coda.samples生成的结果的数据类型为mcmc.list。
n.iter: 每条链的迭代次数
n.thin：保存的比例，比如n.thin=10，意味着每隔10个数据保存一个数据，可有效降低数据量。
4
模型诊断
可以看到估计的mu和sigma比较接近真实值0和5。生成的summary中Rhat和n.eff是衡量MCMC是否收敛和稳健的指标。一般要求Rhat为1，n.eff在10000以上。
MCMCsummary(model)
mean sd 2.5% 50% 97.5% Rhat n.eff
mu 0.04226580 0.164487718 -0.28048634 0.04222556 0.37103327 1 4000
sigma 5.09866516 0.113608046 4.88476288 5.09587009 5.32551275 1 2745
tau 0.03852413 0.001715037 0.03525958 0.03850909 0.04190955 1 2766做个trace plot，看看有没有哪一条链被“困”在某一个地方，生成的posterior的density distribution看起来正不正常。
MCMCtrace(model,ISB = FALSE, ind = TRUE, pdf = FALSE)
MCMC不止于此
JAGS用的是Gibbs sampling的算法，其实MCMC的算法大概包括了:
Metropolis
Metropolis-Hastings
Gibbs （JAGS)
Hamiltonian MCMC （Stan)
No-U-Turn Sampler (NUTS）
具体几个sampler之间的异同，MCMC及原理，这些都是比较底层的东西，不一定会有很多人感兴趣，有机会再做介绍。
贝叶斯统计工具
如果做贝叶斯统计的话，现有的工具都非常容易上手，也不要求你理解MCMC的原理，为什么NUTS的算法会比Gibbs好。比如基于R+Stan的brms，可以看到用法和LM差不多，上手毫无难度。
另外还有比较流行的JASP，大概算是贝叶斯版的免费SPSS吧。今年有个Webinar，可以在官网找到视频教程。
贝叶斯教材推荐
最后推荐一本贝叶斯数据分析的教材：
John Kruschke的Doing Bayesian Data Analysis
没有复杂的公式和晦涩的内容，只不过每章开头一首蹩脚的小诗，看得人不知所云。能坚持看完700多页的内容，其他的三四百页的书都会觉得不在话下（这也许只是个幻觉）。印象最深的还是在power一章讨论两种统计目标，拒绝虚无假设（rejecting the null value)还是达到准确度（achieving precision）时作者说，前者是基于恐惧，害怕如果无法拒绝虚无假设，结果无法发表，后者是基于对真理的追求。像甘地说过的
Power is of two kinds. One is obtained by the fear of punishment and the other by acts of love. Power based on love is a thousand times more effective and permanent than the one derived from fear of punishment.
Mahatma Gandhi
如果觉得搞不定700多页的内容，可以考虑了解一下这篇浓缩版的paper
btw, Ben Lambert的A student's guide to Bayesian Statistics其实应该也还行，没来得及看。作者在油管上有很多视频教程，对公式推导或者英音感兴趣的人不妨前往观摩。
—END—
扫一扫
4138
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
