|
如题,小的在学习贝叶斯模型平均相关内容,在模型后验概率计算上了,用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 )
|