楼主: mercuro
3886 1

[问答] Bayesian model average 模型后验概率计算 [推广有奖]

  • 3关注
  • 0粉丝

已卖:95份资源

硕士生

62%

还不是VIP/贵宾

-

威望
0
论坛币
740 个
通用积分
2.4644
学术水平
0 点
热心指数
0 点
信用等级
0 点
经验
4694 点
帖子
53
精华
0
在线时间
301 小时
注册时间
2014-5-10
最后登录
2024-3-22

楼主
mercuro 学生认证  发表于 2016-12-19 02:08:15 |AI写论文
100论坛币
如题,小的在学习贝叶斯模型平均相关内容,在模型后验概率计算上了,用MCMC引入模型指示变量m,在简单数据可以做,回归数据就做不了了,jags提示m变量不匹配,求大神指点!!!
[mad]

程序:

rm(list=ls())
library(rjags)
# 数据 丢9次钢镚,正面6次
N=9
z=6
y = c( rep(0,N-z) , rep(1,z) )
dataList = list(
  y = y ,
  N = N
)
# 模型
modelString = "
model {
for ( i in 1:N ) {
y ~ dbern( theta )  
}
# 做的是钢镚偏向头,或者花的两种先验模型的比较  设置两个0.25   0.75 两个均值
theta ~ dbeta( omega[m]*(kappa-2)+1 , (1-omega[m])*(kappa-2)+1 )
omega[1] <- .25
omega[2] <- .75
kappa <- 12
# 设两个模型先验取值等可能
m ~ dcat( mPriorProb[] )
mPriorProb[1] <- .5
mPriorProb[2] <- .5
}
"
writeLines( modelString , con="TEMPmodel.txt" )

parameters = c("m")
adaptSteps = 1000             # Number of steps to "tune" the samplers.
burnInSteps = 1000           # Number of steps to "burn-in" the samplers.
nChains = 4                   # Number of chains to run.
numSavedSteps=50000          # Total number of steps in chains to save.
thinSteps=1                   # Number of steps to "thin" (1=keep every step).
nPerChain = ceiling( ( numSavedSteps * thinSteps ) / nChains ) # Steps per chain.
# Create, initialize, and adapt the model:
jagsModel = jags.model( "TEMPmodel.txt" , data=dataList , # inits=initsList ,
                        n.chains=nChains , n.adapt=adaptSteps )

update( jagsModel , n.iter=burnInSteps )

codaSamples = coda.samples( jagsModel , variable.names=parameters ,
                            n.iter=nPerChain , thin=thinSteps )


mcmcMat = as.matrix( codaSamples , chains=TRUE )
m = mcmcMat[,"m"]

pM1 = sum( m == 1 ) / length( m )
pM2 = 1 - pM1


#************************************************************上面是书的例子  能跑出结果
#
rm(list = ls())
library(BMS)
data(attitude)
#借用这个数据,选择两个自变量,组合的话能得到4种模型
bma4c = bms(rating ~complaints + privileges , data = attitude)
pmp.bma(bma4c)#得到4种模型的模型后验概率

#我想用上面例子的层次模型得到相似结果
dataList <- list( N = nrow(attitude),
                  y =attitude$rating,
                  x = attitude$complaints,
                  z = attitude$privileges
)

modelString = "
   model {
   for ( i in 1:N ) {
   y ~ dnorm( mu[i,m] , tau )
   mu[i,1] <- beta0 + beta1 *x
   mu[i,2] <- beta0  
   mu[i,3] <- beta0 + beta2 *z
   mu[i,4] <- beta0 + beta2 *z +beta1 *x

   }
   tau <- 1/100
   #也想照着例子写,但是做回归时候均值与i有关,怎么搞都不好使,每次报错说mu子集不匹配
   mu[] <- equals(m,1)*mu[,1] + equals(m,2)*mu[,2] + equals(m,3)*mu[,3] + equals(m,4)*mu[,4]

   beta0 ~ dnorm(0 , 1.0E-12)
   beta1 ~ dnorm(0 , 1.0E-12)
   beta2 ~ dnorm(0 , 1.0E-12)

   m ~ dcat(mPrior[])
   mPrior[1] <- 0.25
   mPrior[2] <- 0.25
   mPrior[3] <- 0.25
   mPrior[4] <- 0.25
}"


parameters = c("m","theta0","theta1","theta2")
adaptSteps = 1000             # Number of steps to "tune" the samplers.
burnInSteps = 1000           # Number of steps to "burn-in" the samplers.
nChains = 4                   # Number of chains to run.
numSavedSteps=50000          # Total number of steps in chains to save.
thinSteps=1                   # Number of steps to "thin" (1=keep every step).
nPerChain = ceiling( ( numSavedSteps * thinSteps ) / nChains ) # Steps per chain.
# Create, initialize, and adapt the model:
jagsModel = jags.model( "TEMPmodel.txt" , data=dataList , # inits=initsList ,
                        n.chains=nChains , n.adapt=adaptSteps )

update( jagsModel , n.iter=burnInSteps )

codaSamples = coda.samples( jagsModel , variable.names=parameters ,
                            n.iter=nPerChain , thin=thinSteps )


mcmcMat = as.matrix( codaSamples , chains=TRUE )
m = mcmcMat[,"m"]

pM1 = sum( m == 1 ) / length( m )
pM2 =  sum( m == 2 ) / length( m )
pM3 =  sum( m == 3 ) / length( m )
pM4 =  sum( m == 4 ) / length( m )

关键词:Bayesian average Bayes model 后验概率 贝叶斯 JAGS

沙发
courage20 发表于 2020-6-13 11:25:05
你好,我也遇到了类似的问题,希望前辈指点迷津,方便的话加个QQ好友一起讨论。1322267659   期待你的回复

您需要登录后才可以回帖 登录 | 我要注册

本版微信群
加好友,备注cda
拉您进交流群
GMT+8, 2025-12-30 17:48