bayesian - 混合先验在 JAGS 中不起作用，仅当包含似然项时

Question

底部的代码将复制问题，只需将其复制并粘贴到 R 中即可。

我想要的是均值和精度在 30% 的时间内为 (-100, 100)，在 70% 的时间内为 (200, 1000)。将其视为排列在 a、b 和 p 中。

所以“选择”应该是 1 30% 的时间，以及 2 70% 的时间。

实际发生的是在每次迭代中，pick 为 2(如果 p 的第一个元素较大，则为 1)。您可以在摘要中看到这一点，其中“pick”、“testa”和“testb”的分位数始终保持不变。最奇怪的是，如果您删除似然循环，那么 pick then 将完全按预期工作。

我希望这能解释问题，如果不能让我知道。这是我第一次发帖，所以我一定会把事情搞砸的。

library(rjags)
n = 10
y <- rnorm(n, 5, 10)
a = c(-100, 200)
b = c(100, 1000)
p = c(0.3, 0.7)

## Model

mod_str = "model{
    # Likelihood
    for (i in 1:n){
            y[i] ~ dnorm(mu, 10)
    }

    # ISSUE HERE: MIXTURE PRIOR
    mu ~ dnorm(a[pick], b[pick])
    pick ~ dcat(p[1:2])

    testa = a[pick]
    testb = b[pick]
}"

model = jags.model(textConnection(mod_str), data = list(y = y, n=n, a=a, b=b, p=p), n.chains=1)
update(model, 10000)
res = coda.samples(model, variable.names = c('pick', 'testa', 'testb', 'mu'), n.iter = 10000)
summary(res)

Answer 1

我认为您遇到问题有几个原因。首先，您提供给模型的数据(即 y)不是正态分布的混合。因此，模型本身没有必要混合。相反，我会生成这样的数据:

set.seed(320)

# number of samples
n <- 10

# Because it is a mixture of 2 we can just use an indicator variable.
#  here, pick (in the long run), would be '1' 30% of the time.
pick <- rbinom(n, 1, p[1])

# generate the data. b is in terms of precision so we are converting this
#  to standard deviations (which is what R wants).
y_det <- pick * rnorm(n, a[1], sqrt(1/b[1])) + (1 - pick) * rnorm(n, a[2], sqrt(1/b[2]))

# add a small amount of noise, can change to be more as necessary.
y <- rnorm(n, y_det, 1)

这些数据看起来更像是您想要提供给混合模型的数据。

在此之后，我将以与数据生成过程类似的方式对模型进行编码。我想要一些指标变量在两个正态分布之间跳跃。因此，对于 y 中的每个标量，mu 可能会发生变化。

mod_str = "model{
    # Likelihood
    for (i in 1:n){
            y[i] ~ dnorm(mu[i], 10)
            mu[i] <- mu_ind[i] * a_mu + (1 - mu_ind[i]) * b_mu
            mu_ind[i] ~ dbern(p[1])

    }
    a_mu ~ dnorm(a[1], b[1])
    b_mu ~ dnorm(a[2], b[2])

}"

model = jags.model(textConnection(mod_str), data = list(y = y, n=n, a=a, b=b, p=p), n.chains=1)
update(model, 10000)
res = coda.samples(model, variable.names = c('mu_ind', 'a_mu', 'b_mu'), n.iter = 10000)
summary(res)

             2.5%    25%    50%    75% 97.5%
a_mu       -100.4 -100.3 -100.2 -100.1  -100
b_mu        199.9  200.0  200.0  200.0   200
mu_ind[1]     0.0    0.0    0.0    0.0     0
mu_ind[2]     1.0    1.0    1.0    1.0     1
mu_ind[3]     0.0    0.0    0.0    0.0     0
mu_ind[4]     1.0    1.0    1.0    1.0     1
mu_ind[5]     0.0    0.0    0.0    0.0     0
mu_ind[6]     0.0    0.0    0.0    0.0     0
mu_ind[7]     1.0    1.0    1.0    1.0     1
mu_ind[8]     0.0    0.0    0.0    0.0     0
mu_ind[9]     0.0    0.0    0.0    0.0     0
mu_ind[10]    1.0    1.0    1.0    1.0     1

如果您提供更多数据，您将(从长远来看)让指标变量 mu_ind 在 30% 的时间内取值 1。如果您有 2 个以上的发行版，您可以改用 dcat。因此，另一种更通用的方法是 ( and I am borrowing heavily from this post by John Kruschke ):

mod_str = "model {
  # Likelihood:
  for( i in 1 : n ) {
    y[i] ~ dnorm( mu[i] , 10 ) 
    mu[i] <- muOfpick[ pick[i] ]
    pick[i] ~ dcat( p[1:2] )
  }
  # Prior:
  for ( i in 1:2 ) {
    muOfpick[i] ~ dnorm( a[i] , b[i] )
  }
}"
model = jags.model(textConnection(mod_str), data = list(y = y, n=n, a=a, b=b, p=p), n.chains=1)
update(model, 10000)
res = coda.samples(model, variable.names = c('pick', 'muOfpick'), n.iter = 10000)
summary(res)

              2.5%    25%    50%    75% 97.5%
muOfpick[1] -100.4 -100.3 -100.2 -100.1  -100
muOfpick[2]  199.9  200.0  200.0  200.0   200
pick[1]        2.0    2.0    2.0    2.0     2
pick[2]        1.0    1.0    1.0    1.0     1
pick[3]        2.0    2.0    2.0    2.0     2
pick[4]        1.0    1.0    1.0    1.0     1
pick[5]        2.0    2.0    2.0    2.0     2
pick[6]        2.0    2.0    2.0    2.0     2
pick[7]        1.0    1.0    1.0    1.0     1
pick[8]        2.0    2.0    2.0    2.0     2
pick[9]        2.0    2.0    2.0    2.0     2
pick[10]       1.0    1.0    1.0    1.0     1

上面的链接包含更多先验(例如，关于并入分类分布的概率的狄利克雷先验)。