winbugs – JAGS – 层次模型比较,即使是伪造的,也不会在模型之间跳跃

2023年1月4日 182次阅读

我正在使用
Kruschke描述的分层建模框架来设置JAGS中两个模型之间的比较.此框架中的想法是通过将每个版本指定为分类变量的一个级别来运行和比较模型的多个版本.然后,该分类变量的后验分布可以被解释为各种模型的相对概率.

在下面的代码中,我正在比较两个模型.模型形式相同.每个都有一个需要估算的参数,mE.可以看出,这些模型的先验不同.两个先验都作为beta分布分布,模式为0.5.然而,模型2的先前分布更加集中.还要注意我已经使用了伪先生,我曾希望它会阻止链条卡在其中一个模型上.但无论如何,这个模型似乎都被卡住了.

这是模型：

 model {

  m ~ dcat( mPriorProb[] )
  mPriorProb[1] <- .5
  mPriorProb[2] <- .5

  omegaM1[1] <- 0.5      #true prior
  omegaM1[2] <- 0.5      #psuedo prior 
  kappaM1[1] <- 3        #true prior for Model 1
  kappaM1[2] <- 5        #puedo prior for Model 1

  omegaM2[1] <- 0.5      #psuedo prior
  omegaM2[2] <- 0.5      #true prior
  kappaM2[1] <- 5        #puedo  prior for Model 2
  kappaM2[2] <- 10        #true prior for Model 2

  for ( s in 1:Nsubj ) {

    mE1[s] ~ dbeta(omegaM1[m]*(kappaM1[m]-2)+1 , (1-omegaM1[m])*(kappaM1[m]-2)+1 )
    mE2[s] ~ dbeta(omegaM2[m]*(kappaM2[m]-2)+1 , (1-omegaM2[m])*(kappaM2[m]-2)+1 )

    mE[s] <- equals(m,1)*mE1[s] + equals(m,2)*mE2[s]

    z[s] ~ dbin( mE[s] , N[s] )

  }
}

这是相关数据的R代码：

dataList = list(
  z = c(81, 59, 36, 18, 28, 59, 63, 57, 42, 28, 47, 55, 38, 
        30, 22, 32, 31, 30, 32, 33, 32, 26, 13, 33, 30), 
  N = rep(96, 25),
  Nsubj = 25
)

当我运行这个模型时,MCMC将每次迭代都花费在m = 1,并且永远不会跳到m = 2.我已经尝试了许多不同的先验和伪先验的组合,并且似乎找不到组合哪个MCMC会考虑m = 2.我甚至尝试为模型1和2指定相同的先验和伪先验,但这没有帮助.在这种情况下,我希望MCMC在模型之间相当频繁地跳跃,花费大约一半时间考虑一个模型,一半时间考虑另一个模型.然而,JAGS仍然在m = 1时花费了整个时间.我已经运行链长达6000次迭代,这对于像这样的简单模型应该足够长.

如果有人对如何解决这个问题有任何想法,我将非常感激.

干杯,
蒂姆

最佳答案 我无法弄明白这一点,但我认为任何其他人都可能会喜欢以下代码,这将从R与Rjgs(必须安装JAGS)重现问题从头到尾完成.

请注意,由于此示例中只有两个竞争模型,我将m~dcat()更改为m~dbern(),然后在代码中的其他位置将m替换为m 1.我希望这可以改善这种行为,但事实并非如此.另请注意,如果我们指定m的初始值,则无论我们选择哪个值,它都会一直停留在该值,因此我无法正确更新(而不是被一个模型或另一个模型所吸引).对我来说是一个令人头疼的事;可能值得在http://sourceforge.net/p/mcmc-jags/discussion/发布Martyn的眼睛

library(rjags)
load.module('glm')

dataList = list(
  z = c(81, 59, 36, 18, 28, 59, 63, 57, 42, 28, 47, 55, 38, 
        30, 22, 32, 31, 30, 32, 33, 32, 26, 13, 33, 30), 
  N = rep(96, 25),
  Nsubj = 25
)

sink("mymodel.txt")
cat("model {

  m ~ dbern(.5)

  omegaM1[1] <- 0.5      #true prior
  omegaM1[2] <- 0.5      #psuedo prior 
  kappaM1[1] <- 3        #true prior for Model 1
  kappaM1[2] <- 5        #puedo prior for Model 1

  omegaM2[1] <- 0.5      #psuedo prior
  omegaM2[2] <- 0.5      #true prior
  kappaM2[1] <- 5        #puedo  prior for Model 2
  kappaM2[2] <- 10        #true prior for Model 2

    for ( s in 1:Nsubj ) {

    mE1[s] ~ dbeta(omegaM1[m+1]*(kappaM1[m+1]-2)+1 , (1-omegaM1[m+1])*(kappaM1[m+1]-2)+1 )
    mE2[s] ~ dbeta(omegaM2[m+1]*(kappaM2[m+1]-2)+1 , (1-omegaM2[m+1])*(kappaM2[m+1]-2)+1 )


    z[s] ~ dbin( (1-m)*mE1[s] + m*mE2[s] , N[s] )

    }
    }
    ", fill=TRUE)
sink()
inits <- function(){list(m=0)}

params <- c("m")

nc <- 1
n.adapt <-100
n.burn <- 200
n.iter <- 5000
thin <- 1
mymodel <- jags.model('mymodel.txt', data = dataList, inits=inits, n.chains=nc, n.adapt=n.adapt)
update(mymodel, n.burn)
mymodel_samples <- coda.samples(mymodel,params,n.iter=n.iter, thin=thin)
summary(mymodel_samples)
点赞