散文網(wǎng) » 生活 »日常 » 拓端tecdat|R語言和Stan,JAGS：用rstan,rjags建立貝葉斯多元線性回歸預(yù)測選舉數(shù)據(jù)

拓端tecdat|R語言和Stan,JAGS：用rstan,rjags建立貝葉斯多元線性回歸預(yù)測選舉數(shù)據(jù)

2021-07-20 11:39 作者:拓端tecdat 0人讀過 | 我要投稿

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

原文出處：拓端數(shù)據(jù)部落公眾號

本文將介紹如何在R中用rstan和rjags做貝葉斯回歸分析，R中有不少包可以用來做貝葉斯回歸分析，比如最早的（同時也是參考文獻和例子最多的）R2WinBUGS包。這個包會調(diào)用WinBUGS軟件來擬合模型，后來的JAGS軟件也使用與之類似的算法來做貝葉斯分析。然而JAGS的自由度更大，擴展性也更好。近來，STAN和它對應(yīng)的R包rstan一起進入了人們的視線。STAN使用的算法與WinBUGS和JAGS不同，它改用了一種更強大的算法使它能完成WinBUGS無法勝任的任務(wù)。同時Stan在計算上也更為快捷，能節(jié)約時間。

例子

設(shè)Yi為地區(qū)i=1，…，ni=1，…，n從2012年到2016年選舉支持率增加的百分比。我們的模型

式中，Xji是地區(qū)i的第j個協(xié)變量。所有變量均中心化并標(biāo)準(zhǔn)化。我們選擇σ2～InvGamma（0.01,0.01）和α～Normal（0100）作為誤差方差和截距先驗分布，并比較不同先驗的回歸系數(shù)。

加載并標(biāo)準(zhǔn)化選舉數(shù)據(jù)

# 加載數(shù)據(jù)
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的訓(xùn)練集，并對剩余的觀測值進行預(yù)測


test ?<- order(runif(n))>100
table(test)

## test
## FALSE ?TRUE
## ? 100 ?3011

Yo ? ?<- Y[!test] ? ?# 觀測數(shù)據(jù)
Xo ? ?<- X[!test,]


Yp ? ?<- Y[test] ? ? # 為預(yù)測預(yù)留的地區(qū)
Xp ? ?<- X[test,]

選舉數(shù)據(jù)的探索性分析?

boxplot(X, las = 3

image(1:p, 1:p, main = "預(yù)測因子之間的相關(guān)性")

rstan中實現(xiàn)

統(tǒng)一先驗分布

如果模型沒有明確指定先驗分布，默認(rèn)情況下，Stan將在參數(shù)的合適范圍內(nèi)發(fā)出一個統(tǒng)一的先驗分布。注意這個先驗可能是不合適的，但是只要數(shù)據(jù)創(chuàng)建了一個合適的后驗值就可以了。

data {
int<lower=0> n; // 數(shù)據(jù)項數(shù)
int<lower=0> k; // 預(yù)測變量數(shù)
matrix[n,k] X; // 預(yù)測變量矩陣
vector[n] Y; // 結(jié)果向量
}
parameters {
real alpha; // 截距
vector[k] beta; // 預(yù)測變量系數(shù)
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中實現(xiàn)

用高斯先驗擬合線性回歸模型

library(rjags)
model{
# ?預(yù)測
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不發(fā)送給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"),

從后驗預(yù)測分布(PPD)和JAGS預(yù)測分布繪制樣本

#提取每個參數(shù)的樣本
samps ? ? ? <- samp[[1]]
Yp.samps ? ?<- samps[,1:np]
#計算JAGS預(yù)測的后驗平均值
beta.mn ?<- colMeans(beta.samps)
# 繪制后驗預(yù)測分布和JAGS預(yù)測
for(j in 1:5)
# JAGS預(yù)測
y ?<- rnorm(20000,mu,sigma.mn)
plot(density(y),col=2,xlab="Y",main="PPD")
# 后驗預(yù)測分布
lines(density(Yp.samps[,j]))
# 真值
abline(v=Yp[j],col=3,lwd=2)

# 95% 置信區(qū)間
alpha.mn+Xp%*%beta.mn - 1.96*sigma.mn
alpha.mn+Xp%*%beta.mn + 1.96*sigma.mn

## [1] 0.9452009

# PPD 95% 置信區(qū)間
apply(Yp.samps,2,quantile,0.025)
apply(Yp.samps,2,quantile,0.975)

## [1] 0.9634673

請注意，PPD密度比JAGS預(yù)測密度略寬。這是考慮β和σ中不確定性的影響，它解釋了JAGS預(yù)測的covarage略低的原因。但是，對于這些數(shù)據(jù)，JAGS預(yù)測的覆蓋率仍然可以。

最受歡迎的見解

1.matlab使用貝葉斯優(yōu)化的深度學(xué)習(xí)

2.matlab貝葉斯隱馬爾可夫hmm模型實現(xiàn)

3.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真

4.R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸

5.R語言中的Stan概率編程MCMC采樣的貝葉斯模型

6.Python用PyMC3實現(xiàn)貝葉斯線性回歸模型

7.R語言使用貝葉斯層次模型進行空間數(shù)據(jù)分析

8.R語言隨機搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型

9.matlab貝葉斯隱馬爾可夫hmm模型實現(xiàn)

標(biāo)簽：