拓端数据tecdat|R语言和Stan,JAGS：用rstan,rjags建立贝叶斯多元线性回归预测选举数据

原文链接：http://tecdat.cn/?p=21978 

本文将介绍如何在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）作为误差方差和截距先验分布，并比较不同先验的回归系数。

加载并标准化选举数据

1.  
   # 加载数据
2.  
    
3.  
    
4.  
    
5.  
   load("elec.RData")
6.  
    
7.  
   Y <- Y[!is.na(Y+rowSums(X))]
8.  
   X <- X[!is.na(Y+rowSums(X)),]
9.  
   n <- length(Y)
10.  
    p <- ncol(X)

## [1] 3111

 p

## [1] 15

1.  
   X <- scale(X)
2.  
    
3.  
   # 将模型拟合到大小为100的训练集，并对剩余的观测值进行预测
4.  
    
5.  
   test <- order(runif(n))>100
6.  
   table(test)

1.  
   ## test
2.  
   ## FALSE TRUE
3.  
   ## 100 3011

1.  
   Yo <- Y[!test] # 观测数据
2.  
   Xo <- X[!test,]
3.  
    
4.  
   Yp <- Y[test] # 为预测预留的地区
5.  
   Xp <- X[test,]
6.  
    

选举数据的探索性分析 

1.  
    
2.  
   boxplot(X, las = 3

image(1:p, 1:p, main = "预测因子之间的相关性")

rstan中实现

统一先验分布

如果模型没有明确指定先验分布，默认情况下，Stan将在参数的合适范围内发出一个统一的先验分布。注意这个先验可能是不合适的，但是只要数据创建了一个合适的后验值就可以了。

1.  
    
2.  
   data {
3.  
   int<lower=0> n; // 数据项数
4.  
   int<lower=0> k; // 预测变量数
5.  
   matrix[n,k] X; // 预测变量矩阵
6.  
   vector[n] Y; // 结果向量
7.  
   }
8.  
   parameters {
9.  
   real alpha; // 截距
10.  
    vector[k] beta; // 预测变量系数
11.  
    real<lower=0> sigma; // 误差
12.  
     

1.  
   rstan_options(auto_write = TRUE)
2.  
    
3.  
   #fit <- stan(file = 'mlr.stan', data = dat)

print(fit)

hist(fit, pars = pars)

dens(fit)

traceplot(fit)

 

rjags中实现

用高斯先验拟合线性回归模型

1.  
   library(rjags)
2.  
    
3.  
   model{
4.  
   # 预测
5.  
   for(i in 1:np){
6.  
   Yp[i] ~ dnorm(mup[i],inv.var)
7.  
   mup[i] <- alpha + inprod(Xp[i,],beta[])
8.  
    
9.  
   # 先验概率
10.  
     
11.  
    alpha ~ dnorm(0, 0.01)
12.  
    inv.var ~ dgamma(0.01, 0.01)
13.  
    sigma <- 1/sqrt(inv.var)
14.  
     

在JAGS中编译模型

1.  
   # 注意：Yp不发送给JAGS
2.  
   jags.model(model,
3.  
   data = list(Yo=Yo,no=no,np=np,p=p,Xo=Xo,Xp=Xp))

1.  
   coda.samples(model,
2.  
   variable.names=c("beta","sigma","Yp","alpha"),
3.  
    

从后验预测分布(PPD)和JAGS预测分布绘制样本

1.  
   #提取每个参数的样本
2.  
    
3.  
   samps <- samp[[1]]
4.  
   Yp.samps <- samps[,1:np]
5.  
    
6.  
   #计算JAGS预测的后验平均值
7.  
    
8.  
   beta.mn <- colMeans(beta.samps)
9.  
    
10.  
     
11.  
    # 绘制后验预测分布和JAGS预测
12.  
     
13.  
    for(j in 1:5)
14.  
    # JAGS预测
15.  
    y <- rnorm(20000,mu,sigma.mn)
16.  
    plot(density(y),col=2,xlab="Y",main="PPD")
17.  
     
18.  
    # 后验预测分布
19.  
    lines(density(Yp.samps[,j]))
20.  
     
21.  
    # 真值
22.  
    abline(v=Yp[j],col=3,lwd=2)

1.  
   # 95% 置信区间
2.  
   alpha.mn+Xp%*%beta.mn - 1.96*sigma.mn
3.  
   alpha.mn+Xp%*%beta.mn + 1.96*sigma.mn

## [1] 0.9452009

1.  
   # PPD 95% 置信区间
2.  
   apply(Yp.samps,2,quantile,0.025)
3.  
   apply(Yp.samps,2,quantile,0.975)
4.  
    

## [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模型实现