您如何在 JAGS/BUGS 中编写 Tweedie 分布？

Question

我想运行使用 JAGS 通过 R 在 Tweedie 分布式变量上建立模型。我知道 JAGS 没有标准的 Tweedie 分布，但可以指定一个作为化合物 Gamma/Poisson。不幸的是，我无法弄清楚如何在 JAGS 中对其进行编码。我根据从各种来源收集的代码编写了以下代码，以简单地尝试从 Tweedie 随机变量中恢复均值、功率和 phi 参数。由于 y 上的无效父值，它目前没有运行，大概是因为 y[i] 出现在表达式的右侧和左侧。这是在源代码中写的，但我显然滥用了它。任何关于如何正确指定此分布的指示都将不胜感激，并且可能会得到更广泛的使用，因为我无法找到关于如何在 JAGS 中设置 Tweedie 模型的示例的任何简单编码。

y = mgcv::rTweedie(mu=rep(2,100),p=1.33,phi=1)

jags_data = list(y=y,n=length(y))

jags_model = 
 "model{
    
    for (i in 1:n) {
      lambda[i] <- pow(mu,2-p)/(phi *(2-p))
      num[i] ~ dpois(lambda[i])
      shape[i,1] <- num[i]*((2-p)/(p-1))
      rate[i,1] <- 1/(phi*(p-1)*pow(mu,p-1))
      shape[i,2] <- 1
      rate[i,2] <- exp(-lambda[i])
      # Takes shape/rate parameter 1 if y > 0 and 2 if y = 0
      y[i] ~ dgamma(shape[i,1+equals(y[i],0)],rate[i,1+equals(y[i],0)]) 
    }
    
    mu    ~ dunif(0,100)
    p     ~ dunif(1,2)  ## Tweedie power parameter
    phi   ~ dunif(0,30) ## Dispersion parameter
    
  }
  
  "

model_file = tempfile(fileext = 'txt')
writeLines(jags_model,model_file)

jm = rjags::jags.model(
  file = model_file,
  data = jags_data,
  n.chains = 3,
  n.adapt = 1500
)

Answer 1

我能够让它发挥作用。它是否“有效”似乎更像是一个悬而未决的问题，但参数的后验均值非常接近真实值。我在 R 中使用 runjags 完成了此操作，但此处的所有部分都将在独立于 R 的 JAGS 中完成。

首先，我生成了数据，并将形状矩阵放入数据中，其中有一列 NA 值，这样它们就可以在模拟的每次迭代中被覆盖。我还添加了一个名为 yind 的变量，如果 y == 0 则为 2，否则为 1。这将用作索引 shape 和 rate 矩阵的值。在数据中执行此操作可能更有效，而不是让 JAGS 在每次迭代中针对开始时固定的 y 的每个值计算和 equals() 函数几次。

y = mgcv::rTweedie(mu=rep(2,100),p=1.33,phi=1)
shape_mat <- matrix(NA, nrow=length(y), ncol=2)
shape_mat[,2] <- 1
jags_data = list(y=y,n=length(y), yind = 2-(y > 0), shape = shape_mat)

接下来，模型是一样的，只是改变了形状矩阵的处理方式和添加了 yind。

jags_model = 
  "model{
    
    for (i in 1:n) {
      lambda[i] <- pow(mu,2-p)/(phi *(2-p))
      num[i] ~ dpois(lambda[i])
      shape[i,1] <- num[i]*((2-p)/(p-1))
      rate[i,1] <- 1/(phi*(p-1)*pow(mu,p-1))
      ## moved to data
      # shape[i,2] <- 1
      rate[i,2] <- exp(-lambda[i])
      # Takes shape/rate parameter 1 if y > 0 and 2 if y = 0
      y[i] ~ dgamma(shape[i,yind[i]],rate[i,yind[i]]) 
    }
    
    mu    ~ dunif(0,100)
    p     ~ dunif(1,2)  ## Tweedie power parameter
    phi   ~ dunif(0,30) ## Dispersion parameter
    
  }
  "

无效的父值可能来自初始值。 num 泊松图必须与 lambda 一致，后者是 mu、p 和 phi 的函数。因此，初始值对于确保这些都是一致的很重要。我将三个模型参数设置为下面的值，然后计算lambda。我根据这三个参数的值将num设置为最接近lambda的整数值

inits <- list(mu = 5, p=1.5, phi=5,  
              num = rep(1, length(y)))

接下来，我运行模型和总结：

library(runjags)
out <- run.jags(model=jags_model, data = jags_data, monitor=c("mu", "p", "phi"), burnin=5000, sample=10000, 
                inits = list(inits, inits), n.chains = 2, keep.jags.files=TRUE)
summary(out)
# Lower95   Median Upper95     Mean         SD Mode       MCerr MC%ofSD SSeff        AC.10     psrf
# mu  1.606100 1.921855 2.24793 1.927986 0.16489395   NA 0.001534324     0.9 11550 -0.006514202 1.000242
# p   1.252230 1.377415 1.50807 1.379773 0.06563517   NA 0.002049511     3.1  1026  0.354243967 1.013941
# phi 0.871354 1.098075 1.36870 1.109978 0.12874006   NA 0.003880090     3.0  1101  0.322202621 1.004731

您如何在 JAGS/BUGS 中编写 Tweedie 分布？

How do you code a Tweedie distribution in JAGS/BUGS?

r

jags

tweedie