全部版块 我的主页
› 论坛 › 数据科学与人工智能 › 数据分析与数据科学 › R语言论坛
Bayesian model average 模型后验概率计算
mercuro
3811 1
2016-12-19
悬赏 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 )

二维码

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

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

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

相关推荐
栏目导航
热门文章
推荐文章

分享

微信QQ空间QQ微博
扫码加好友，拉您进群
各岗位、行业、专业交流群