Как вы кодируете дистрибутив Tweedie в JAGS/BUGS?

Я хотел бы запустить модель с распределенной переменной Tweedie, используя JAGS через R. Я знаю, что JAGS не имеет стандартного распределения Tweedie, но его можно указать как составную гамму/Пуассона. К сожалению, я не могу понять, как это закодировать в JAGS. Я написал ниже на основе кода, извлеченного из различных источников, чтобы просто попытаться восстановить параметры среднего значения, мощности и фи из случайной величины Tweedie. В настоящее время он не выполняется из-за недопустимых родительских значений y, предположительно потому, что y[i] появляется в правой и левой частях выражения. Это так, как было написано в исходном коде, но я явно неправильно его использую. Любые указания о том, как правильно указать это распределение, были бы высоко оценены и, вероятно, имели бы более широкое применение, поскольку я не смог найти каких-либо простых закодированных примеров того, как настроить модели Tweedie в JAGS.

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
)

r jags tweedie
person kernowsam    schedule 18.08.2020    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Я смог заставить это работать. Сработало это или нет, кажется открытым вопросом, но апостериорные средние значения параметров довольно близки к истинным значениям. Я сделал это в R с помощью runjags, но все части здесь должны быть выполнены в JAGS независимо от R.

Сначала я сгенерировал данные, а также поместил матрицу формы в данные со столбцом значений NA, чтобы их можно было перезаписывать при каждой итерации моделирования. Я также добавил переменную с именем yind, которая равна 2, если y == 0, и 1 в противном случае. Это будет служить значением, которое будет индексировать матрицы shape и rate. Вероятно, более эффективно делать это в данных, а не заставлять JAGS оценивать и выполнять функцию equals() несколько раз на каждой итерации для каждого значения 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, которые являются функцией 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
person DaveArmstrong    schedule 02.09.2020