4

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

在下面的代码中，我比较了两个模型。这些模型在形式上是相同的。每个都有一个需要估计的参数，mE。可以看出，这些模型的先验有所不同。两个先验都分布为众数为 0.5 的 beta 分布。然而，模型 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 次迭代的链，这对于像这样的简单模型来说应该足够长了。

如果有人对如何解决此问题有任何想法，我将不胜感激。

干杯，蒂姆

4

3 回答 3

1

诀窍不是为模型分配固定概率，而是phi基于统一的先验估计它（如下）。然后，您需要后验分布，phi因为它告诉您选择模型 2 的概率（即，“成功”意味着 m=1；Pr(model 1) = 1-phi）。

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

  m ~ dbern(phi)
  phi ~ dunif(0,1)

  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("phi")
于 2016-01-13T00:51:46.920 回答
1

我无法弄清楚这一点，但我认为从事此工作的其他任何人都可能会欣赏以下代码，该代码将使用 rjags 从 R 开始从头到尾重现问题（必须安装 JAGS）。

请注意，由于此示例中只有两个竞争模型，因此我将其更改m ~ dcat()为m ~ dbern()，然后在代码中的其他任何地方m替换。m+1我希望这可能会改善这种行为，但事实并非如此。另请注意，如果我们为 m 指定初始值，无论我们选择哪个值，它都会停留在该值上，因此 m 无法正确更新（而不是奇怪地被一个模型或另一个模型吸引）。让我头疼；可能值得在http://sourceforge.net/p/mcmc-jags/discussion/上张贴马丁的眼睛

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)
于 2015-11-17T01:59:03.283 回答
0

请参阅我上面对 Mark S 的回答的评论。

这个答案是为了举例说明为什么我们要推断 m 而不是 phi。

想象一下，我们有一个模型

data <- c(-1, 0, 1, .5, .1)

m~dbern(phi)
data[i] ~ m*dnorm(0, 1) + (1-m)*dnorm(100, 1)

现在，很明显 m 的真实值为 1。但我们对 phi 的真实值了解多少？显然更高的 phi 值更有可能，但我们实际上并没有很好的证据来排除更低的 phi 值。例如，phi=0.1 仍然有 10% 的机会产生 m=1；并且 phi=0.5 仍然有 50% 的机会产生 m=1。所以我们没有很好的证据来反对相当低的 phi 值，即使我们有确凿的证据证明 m=1。我们想要对 m 进行推断。

于 2016-04-18T04:10:59.697 回答