原文链接：http://tecdat.cn/?p=21978
原文出处：拓端数据部落公众号
最近咱们被客户要求撰写对于Stan,JAGS的钻研报告，包含一些图形和统计输入。
本文将介绍如何在R中用rstan和rjags做贝叶斯回归剖析，R中有不少包能够用来做贝叶斯回归剖析，比方最早的（同时也是参考文献和例子最多的）R2WinBUGS包。这个包会调用WinBUGS软件来拟合模型，起初的JAGS软件也应用与之类似的算法来做贝叶斯剖析。然而JAGS的自由度更大，扩展性也更好。近来，STAN和它对应的R包rstan一起进入了人们的眼帘。STAN应用的算法与WinBUGS和JAGS不同，它改用了一种更弱小的算法使它能实现WinBUGS无奈胜任的工作。同时Stan在计算上也更为快捷，能节约工夫。
例子
设Yi为地区i=1，…，ni=1，…，n从2012年到2016年选举支持率减少的百分比。咱们的模型
式中，Xji是地区i的第j个协变量。所有变量均中心化并标准化。咱们抉择σ2∼InvGamma（0.01,0.01）和α∼Normal（0100）作为误差方差和截距先验散布，并比拟不同先验的回归系数。
加载并标准化选举数据
# 加载数据
load("elec.RData")
Y <- Y[!is.na(Y+rowSums(X))]
X <- X[!is.na(Y+rowSums(X)),]
n <- length(Y)
p <- ncol(X)## [1] 3111 p## [1] 15 X <- scale(X)
# 将模型拟合到大小为100的训练集，并对残余的观测值进行预测
test <- order(runif(n))>100
table(test)## test
## FALSE TRUE
## 100 3011 Yo <- Y[!test] # 观测数据
Xo <- X[!test,]
Yp <- Y[test] # 为预测预留的地区
Xp <- X[test,]选举数据的探索性剖析
boxplot(X, las = 3image(1:p, 1:p, main = "预测因子之间的相关性")rstan中实现
对立先验散布
如果模型没有明确指定先验散布，默认状况下，Stan将在参数的适合范畴内收回一个对立的先验散布。留神这个先验可能是不适合的，然而只有数据创立了一个适合的后验值就能够了。
data {
int<lower=0> n; // 数据项数
int<lower=0> k; // 预测变量数
matrix[n,k] X; // 预测变量矩阵
vector[n] Y; // 后果向量
}
parameters {
real alpha; // 截距
vector[k] beta; // 预测变量系数
real<lower=0> sigma; // 误差rstan_options(auto_write = TRUE)
#fit <- stan(file = 'mlr.stan', data = dat)print(fit)
hist(fit, pars = pars)dens(fit)traceplot(fit)
rjags中实现
用高斯先验拟合线性回归模型
library(rjags)
model{
# 预测
for(i in 1:np){
Yp[i] ~ dnorm(mup[i],inv.var)
mup[i] <- alpha + inprod(Xp[i,],beta[])
# 先验概率
alpha ~ dnorm(0, 0.01)
inv.var ~ dgamma(0.01, 0.01)
sigma <- 1/sqrt(inv.var)在JAGS中编译模型
# 留神：Yp不发送给JAGS
jags.model(model,
data = list(Yo=Yo,no=no,np=np,p=p,Xo=Xo,Xp=Xp))coda.samples(model,
variable.names=c("beta","sigma","Yp","alpha"), 从后验预测散布(PPD)和JAGS预测散布绘制样本
#提取每个参数的样本
samps <- samp[[1]]
Yp.samps <- samps[,1:np]
#计算JAGS预测的后验平均值
beta.mn <- colMeans(beta.samps)
# 绘制后验预测散布和JAGS预测
for(j in 1:5)
# JAGS预测
y <- rnorm(20000,mu,sigma.mn)
plot(density(y),col=2,xlab="Y",main="PPD")
# 后验预测散布
lines(density(Yp.samps[,j]))
# 真值
abline(v=Yp[j],col=3,lwd=2)
# 95% 置信区间
alpha.mn+Xp%*%beta.mn - 1.96*sigma.mn
alpha.mn+Xp%*%beta.mn + 1.96*sigma.mn## [1] 0.9452009 # PPD 95% 置信区间
apply(Yp.samps,2,quantile,0.025)
apply(Yp.samps,2,quantile,0.975)## [1] 0.9634673请留神，PPD密度比JAGS预测密度略宽。这是思考β和σ中不确定性的影响，它解释了JAGS预测的covarage略低的起因。然而，对于这些数据，JAGS预测的覆盖率依然能够。
最受欢迎的见解
1.matlab应用贝叶斯优化的深度学习
2.matlab贝叶斯隐马尔可夫hmm模型实现
3.R语言Gibbs抽样的贝叶斯简略线性回归仿真
4.R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归
5.R语言中的Stan概率编程MCMC采样的贝叶斯模型
6.Python用PyMC3实现贝叶斯线性回归模型
7.R语言应用贝叶斯 层次模型进行空间数据分析
8.R语言随机搜寻变量抉择SSVS预计贝叶斯向量自回归（BVAR）模型
9.matlab贝叶斯隐马尔可夫hmm模型实现
发表回复