R語(yǔ)言和Stan,JAGS：用rstan,rjags建立貝葉斯多元線性回歸預(yù)測(cè)選舉數(shù)據(jù)|附代碼數(shù)據(jù)
原文鏈接：http://tecdat.cn/?p=21978?
原文出處：拓端數(shù)據(jù)部落公眾號(hào)
?最近我們被客戶要求撰寫(xiě)關(guān)于Stan,JAGS的研究報(bào)告，包括一些圖形和統(tǒng)計(jì)輸出。
本文將介紹如何在R中用rstan和rjags做貝葉斯回歸分析，R中有不少包可以用來(lái)做貝葉斯回歸分析，比如最早的（同時(shí)也是參考文獻(xiàn)和例子最多的）R2WinBUGS包。這個(gè)包會(huì)調(diào)用WinBUGS軟件來(lái)擬合模型，后來(lái)的JAGS軟件也使用與之類(lèi)似的算法來(lái)做貝葉斯分析。然而JAGS的自由度更大，擴(kuò)展性也更好。近來(lái)，STAN和它對(duì)應(yīng)的R包rstan一起進(jìn)入了人們的視線。STAN使用的算法與WinBUGS和JAGS不同，它改用了一種更強(qiáng)大的算法使它能完成WinBUGS無(wú)法勝任的任務(wù)。同時(shí)Stan在計(jì)算上也更為快捷，能節(jié)約時(shí)間。
例子
設(shè)Yi為地區(qū)i=1，…，ni=1，…，n從2012年到2016年選舉支持率增加的百分比。我們的模型
式中，Xji是地區(qū)i的第j個(gè)協(xié)變量。所有變量均中心化并標(biāo)準(zhǔn)化。我們選擇σ2～I(xiàn)nvGamma（0.01,0.01）和α～Normal（0100）作為誤差方差和截距先驗(yàn)分布，并比較不同先驗(yàn)的回歸系數(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)練集，并對(duì)剩余的觀測(cè)值進(jìn)行預(yù)測(cè) test ?<- order(runif(n))>100 table(test)
## test## FALSE ?TRUE ## ? 100 ?3011
Yo ? ?<- Y[!test] ? ?# 觀測(cè)數(shù)據(jù) Xo ? ?<- X[!test,] Yp ? ?<- Y[test] ? ? # 為預(yù)測(cè)預(yù)留的地區(qū) Xp ? ?<- X[test,]
選舉數(shù)據(jù)的探索性分析?
boxplot(X, las = 3
image(1:p, 1:p, main = "預(yù)測(cè)因子之間的相關(guān)性")
rstan中實(shí)現(xiàn)
統(tǒng)一先驗(yàn)分布
如果模型沒(méi)有明確指定先驗(yàn)分布，默認(rèn)情況下，Stan將在參數(shù)的合適范圍內(nèi)發(fā)出一個(gè)統(tǒng)一的先驗(yàn)分布。注意這個(gè)先驗(yàn)可能是不合適的，但是只要數(shù)據(jù)創(chuàng)建了一個(gè)合適的后驗(yàn)值就可以了。
data { ?int<lower=0> n; // 數(shù)據(jù)項(xiàng)數(shù) ?int<lower=0> k; // 預(yù)測(cè)變量數(shù) ?matrix[n,k] X; // 預(yù)測(cè)變量矩陣 ?vector[n] Y; // 結(jié)果向量}parameters { ?real alpha; // 截距 ?vector[k] beta; // 預(yù)測(cè)變量系數(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中實(shí)現(xiàn)
用高斯先驗(yàn)擬合線性回歸模型
library(rjags)model{# ?預(yù)測(cè) ?for(i in 1:np){ ? ?Yp[i] ?~ dnorm(mup[i],inv.var) ? ?mup[i] <- alpha + inprod(Xp[i,],beta[]) ?# 先驗(yàn)概率 ?alpha ? ? ~ dnorm(0, 0.01) ?inv.var ? ~ dgamma(0.01, 0.01) ?sigma ? ? <- 1/sqrt(inv.var)
在JAGS中編譯模型
# 注意：Yp不發(fā)送給JAGSjags.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àn)預(yù)測(cè)分布(PPD)和JAGS預(yù)測(cè)分布繪制樣本
#提取每個(gè)參數(shù)的樣本 samps ? ? ? <- samp[[1]] Yp.samps ? ?<- samps[,1:np] #計(jì)算JAGS預(yù)測(cè)的后驗(yàn)平均值 beta.mn ?<- colMeans(beta.samps)# 繪制后驗(yàn)預(yù)測(cè)分布和JAGS預(yù)測(cè) for(j in 1:5) ? ?# JAGS預(yù)測(cè) ? ?y ?<- rnorm(20000,mu,sigma.mn) ? ?plot(density(y),col=2,xlab="Y",main="PPD") ? ?# 后驗(yàn)預(yù)測(cè)分布 ? ?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
請(qǐng)注意，PPD密度比JAGS預(yù)測(cè)密度略寬。這是考慮β和σ中不確定性的影響，它解釋了JAGS預(yù)測(cè)的covarage略低的原因。但是，對(duì)于這些數(shù)據(jù)，JAGS預(yù)測(cè)的覆蓋率仍然可以。
最受歡迎的見(jiàn)解
1.matlab使用貝葉斯優(yōu)化的深度學(xué)習(xí)
2.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)
3.R語(yǔ)言Gibbs抽樣的貝葉斯簡(jiǎn)單線性回歸仿真
4.R語(yǔ)言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
5.R語(yǔ)言中的Stan概率編程MCMC采樣的貝葉斯模型
6.Python用PyMC3實(shí)現(xiàn)貝葉斯線性回歸模型
7.R語(yǔ)言使用貝葉斯 層次模型進(jìn)行空間數(shù)據(jù)分析
8.R語(yǔ)言隨機(jī)搜索變量選擇SSVS估計(jì)貝葉斯向量自回歸（BVAR）模型
9.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)
