R語言和STAN,JAGS：用RSTAN,RJAG創建貝葉斯多元線性迴歸預測選舉數據

原文連接：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模型實現