首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >Rjags中的分层混合模型

Rjags中的分层混合模型
EN

Stack Overflow用户
提问于 2019-08-09 06:39:28
回答 1查看 255关注 0票数 0

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

代码语言:javascript
复制
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模型,我一直在玩它。然而,当在简单的模拟数据上进行测试时,我得到了非常糟糕的结果。下面是生成模拟数据并运行模型的代码:

代码语言:javascript
复制
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上的经验。

EN

回答 1

Stack Overflow用户

发布于 2019-08-30 01:41:28

在你的JAGS模型中,lambda依赖于i,这可能不是你的本意。该模型不能使用这样的定义来估计参数,因为基本上有太多的lambda,或者换句话说,每次绘制都遵循它自己的分布。

也许模型应该看起来更像

代码语言:javascript
复制
lambda1 ~ dgamma(a,b)
lambda2 ~ dlnorm(c,d)

或者,如果您真的想要单独的lambda,那么您将需要更多的每个lambda的绘制,以便后验分布有意义。现在,您仅使用每个lambda一个数据点进行估计。

同样的情况也适用于你对混合物使用的分类分布。

票数 0
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/57421481

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档