首页
学习
活动
专区
工具
TVP
发布
精选内容/技术社群/优惠产品,尽在小程序
立即前往

R语言和STAN,JAGS:用RSTAN,RJAG建立贝叶斯多元线性回归预测选举数据

这个包会调用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

X

n

p

## [1] 3111

p

## [1] 15

X

# 将模型拟合到大小为100的训练集,并对剩余的观测值进行预测

test 100

table(test)

## test

## FALSE TRUE

## 100 3011

Yo

Xo

Yp

Xp

选举数据的探索性分析 

boxplot(X, las = 3

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

rstan中实现

统一先验分布

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

data {

int n; // 数据项数

int k; // 预测变量数

matrix[n,k] X; // 预测变量矩阵

vector[n] Y; // 结果向量

}

parameters {

real alpha; // 截距

vector[k] beta; // 预测变量系数

real sigma; // 误差

rstan_options(auto_write = TRUE)

#fit

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 ~ dnorm(0, 0.01)

inv.var ~ dgamma(0.01, 0.01)

sigma

在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

Yp.samps

#计算JAGS预测的后验平均值

beta.mn

# 绘制后验预测分布和JAGS预测

for(j in 1:5)

# JAGS预测

y

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预测的覆盖率仍然可以。

  • 发表于:
  • 原文链接https://page.om.qq.com/page/Oc30fBC3cU4JcGNtQiR5i7yw0
  • 腾讯「腾讯云开发者社区」是腾讯内容开放平台帐号(企鹅号)传播渠道之一,根据《腾讯内容开放平台服务协议》转载发布内容。
  • 如有侵权,请联系 cloudcommunity@tencent.com 删除。

相关快讯

扫码

添加站长 进交流群

领取专属 10元无门槛券

私享最新 技术干货

扫码加入开发者社群
领券