0

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

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
)
4

1 回答 1

1

我能够让它发挥作用。它是否“有效”似乎更像是一个悬而未决的问题，但参数的后验均值非常接近真实值。我在 Rrunjags中使用 .

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

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这是和的函数。因此，初始值对于确保这些都是一致的很重要。我将三个模型参数设置为以下值，然后计算。我根据这三个参数的值设置为最接近的整数值。mupphilambdanumlambda

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
于 2020-09-02T19:37:12.150 回答