Как определить усеченную модель Пуассона в JAGS? - PullRequest
Новая
       123

Как определить усеченную модель Пуассона в JAGS?

0 голосов
/ 24 октября 2019

Я хочу использовать ту же модель для подсчета данных с изменяющимися нижними порогами измерения. Усеченная модель Пуассона позволяет мне сделать это.

Модель может быть успешно определена в R. Но я не могу перевести ту же модель в JAGS. Кажется, проблема в том, что невозможно разделить PDF на CDF при определении вероятности.

Это то, что я пробовал:

Данные 

set.seed(42)
N <- 10000
x1 <- runif(N, -1, 1)
lower_limit = 7
upper_limit = Inf
b0 <- 0.5
b1 <- 1.3
b2 <- 0.1
true_lambda <- exp(b0 + b1 * x1 + b2 * x1^2)
# extraDistr::dtpois behaves exactly like the truncated poisson defined below
y <- extraDistr::rtpois(N, lambda = true_lambda, a = lower_limit, b = upper_limit)
summary(y)

Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
8.000   8.000   8.000   8.354   8.000  16.000 

hist(y)

гистограмма данных

Для модели я использую усеченное распределение Пуассона, которое можно определить с основанием R следующим образом: 

trunc <- function(x,s, lambda, log){
  out <- dpois(x, lambda)/(1-ppois(s, lambda))
  if(log){
    out <- log(out)
  }
  return(out)
}


optfun <- function(data, par){
  b0 <- par[1]
  b1 <- par[2]
  b2 <- par[3]
  loglikelihood <- sum(trunc(
    x = data,
    lambda = exp(b0 + b1 * x1 + b2 * x1^2),
    s = lower_limit,
    log = T
  ))
  return(-loglikelihood)
}

output <- optim(par = rep(0, 3),
                fn = optfun2,
                data = y)
output$par

[1] 0.4808197 1.3141579 0.1097217

Однако я хочуопределить ту же модель в JAGS: 

Model.Data <- list(
  y = y,
  N = N,
  s = rep(lower_limit,N),
  x1 = x1
)

sink("trunc_pois.txt")
cat(
  "
  model{

    # Likelihood data model:
    for (i in 1:N) {
      # truncated poisson in JAGS:
      y[i] ~ dpois(lambda_fit[i])/(1-ppois(s[i],lambda_fit[i]))
      log(lambda_fit[i]) <- b0 + b1 * x1[i] + b2 * pow(x1[i],2)
    }

    # Priors:
    b0 ~ dnorm(0, 0.0001) # mean, precision = N(0, 10^4)
    b1 ~ dnorm(0, 0.0001)
    b2 ~ dnorm(0, 0.0001)
  }
"
  )

sink()


inits.fn <- function() list(
  b0 = rnorm(1),
  b1 = rnorm(1),
  b2 = rnorm(1)
)

jagsModel <- rjags::jags.model(file = "trunc_pois.txt", data = Model.Data,
                        inits = inits.fn,n.chains = 3, n.adapt = 100)

Это дает мне следующее сообщение об ошибке: 

Error in rjags::jags.model(file = "trunc_pois.txt", data = Model.Data,  : 

Error parsing model file:
syntax error on line 7 near "/"

Есть ли способ определить усеченную модель Пуассона в JAGS? Должен ли я использовать T () для усечения модели Пуассона?

...