0

这是我第一次使用 rjags，我正在尝试拟合一些计数数据 Y。我使用的分层混合模型如下：

Y ~ p*Poisson(N*lambda1) + (1-p)*Poisson(N*lambda2)
lambda1 ~ Gamma(a,b)
lambda2 ~ Lognormal(c,d)
a ~ Gamma(1,1)
b ~ Gamma(1,1)
c ~ Normal(0,1)
d ~ Gamma(1,1)

这里，Y是我观察到的计数数据，N是已知的。

我写了一个我一直在玩的简单的 rjags 模型。但是，在对简单的模拟数据进行测试时，我得到了非常糟糕的结果。这是我生成模拟数据并运行模型的代码：

a <- 0.5
b <- 0.5
c <- -10
d <- 1

lambda1 <- rgamma(30,a,b)
lambda2 <- rlnorm(70,c,d)
counts <- rpois(100,1000*c(lambda1,lambda2))

model_string <- "model{
  # Likelihood 
  for (i in 1:n) {
    mu1[i] <- N*lambda1[i]
    mu2[i] <- N*lambda2[i]
    lambda1[i] ~ dgamma(a,b)
    lambda2[i] ~ dlnorm(c,d)
    m[i] ~ dcat(mprior[])
    mu[i] <- equals(m[i],1)*mu1[i] + equals(m[i],2)*mu2[i]
    Y[i] ~ dpois(mu[i])
  }
  # Prior
  mprior[1] <- 0.5
  mprior[2] <- 0.5
  a ~ dgamma(1,1)
  b ~ dgamma(1,1)
  c ~ dnorm(0,1)
  d ~ dgamma(1,1)
}"

model <- jags.model(textConnection(model_string), 
                    data = list(Y=counts,N=1000,n=100))
update(model,10000)
samp <- coda.samples(model, 
                     variable.names=c("a","b","c","d","m"), 
                     n.iter=20000)

print(colMeans(samp[[1]])[1:4])

运行此之后，a、b、c、d 的后验估计非常差，并且组件分配 m 也与真实分配不匹配。我还注意到绘制链看起来并不好，即使增加迭代次数也是如此。

有什么建议吗？我不确定我是否以最佳方式拟合混合物。如果还有其他更易于使用的发行版，我也绝对愿意更改我在 a、b、c、d 上的先验知识。

4

1 回答 1

0

在您的 JAGS 模型中，lambda 取决于 i，这可能不是您的意图。该模型无法使用这样的定义来估计参数，因为基本上有太多的 lambda 或者换一种方式，每次绘制都遵循它自己的分布。

也许模型应该看起来更像

lambda1 ~ dgamma(a,b)
lambda2 ~ dlnorm(c,d)

或者，如果您真的想拥有单独的 lambda，那么每个 lambda 需要更多的绘制才能使后验分布有意义。现在，您只是用每个 lambda 一个数据点进行估计。

这同样适用于您用于混合的分类分布。

于 2019-08-29T17:41:28.110 回答