乐趣区

关于算法:R语言和STANJAGS用RSTANRJAG建立贝叶斯多元线性回归预测选举数据

jiezi

原文链接：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）作为误差方差和截距先验散布，并比拟不同先验的回归系数。

加载并标准化选举数据

# 加载数据



 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 = 3

​

image(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 模型实现

退出移动版