Байесовское моделирование в пространстве состояний JAGS

Я пытаюсь использовать модель пространства состояний для оценки демографии населения (плодовитость, выживаемость, прирост населения, размер населения). У нас есть 4 разных возрастных состояния.

# J0 = number of individuals 0-1
# surv1 = survivorship from 0-1 

# J1 = number of individuals 0-1 
# surv2 = survivorship from 1-2 

# J2= = number of individuals 0-1
# surv3 = survivorship from 2-3 

# J3= number of individuals 0-1
# survad = survivorship >3 "adult")

# Data as vectors (Talek clan from 1988-2013)
# X0 = individuals 0-1 in years
# X1 = individuals 1-2 in years
# X2 = individuals 2-3 in years
# X3 = individuals 3+ in years
# Total = group size

X0 <- c(7,  9, 4, 8, 9, 5, 8, 5, 7, 5, 5, 8, 10, 3, 5, 7, 2, 6, 6, 11, 14, 12, 15, 9, 10)
X1 <- c( 4, 4, 3, 4, 8, 5, 2, 4, 3, 4, 4, 5, 3, 7, 0, 5, 6, 3, 3, 5, 10, 12, 10, 13, 8)
X2 <- c(3, 2, 3, 3, 3, 8, 4, 1, 1, 2, 2, 4, 2, 2, 5, 0, 5, 5, 4, 3, 3, 10, 12, 7, 10)
X3 <- c(18, 16, 13, 16, 29, 29, 26, 22, 21, 18, 16, 15, 16, 15, 11, 14, 9, 12, 16, 18, 21, 23, 33, 32, 31)
Total <- c(32, 31, 23, 31, 49, 47, 40, 32, 32, 29, 27, 32, 31, 27, 21, 26, 22, 26, 29, 37, 48, 57, 70, 61, 59)

Вот код ОШИБКИ:

sink(file = "HyenaIPM_all.txt")
cat("
model {

# Specify the priors for all parameters in the model
  N.est[1] ~ dnorm(50, tau.proc)T(0,)                               # Initial abundance
  mean.lambda ~ dunif(0, 5)
  sigma.proc ~ dunif(0, 50)
  tau.proc <- pow(sigma.proc, -2)
  for (t in 1:TT) {
    fec[t] ~ dunif(0, 5)                # per capita fecundidty
    surv1[t] ~ dunif(0, 1)          # survivorship from 0-1
    surv2[t] ~ dunif(0, 1)          # survivorship from 1-2
    surv3[t] ~ dunif(0, 1)          # survivorship from 2-3
    survad[t] ~ dunif(0, 1)         # adult survivorship 
    }

# Estimate fecundity and survivorship
    for (t in 2:TT) {

        # Fecundity
        J0[t+1] ~ dpois(survad[t]*fec[t])
        J0[t+1] <- J3[t] * fec[t]

        # Survivorship
        J1[t+1] ~ dbin(surv1[t], J0[t])
        J1[t+1] <- J0[t]*surv1[t]

        J2[t+1] ~ dbin(surv2[t], J1[t])
        J2[t+1] <- J1[t]*surv2[t]

        J3[t+1] ~ dbin(surv3[t], J2[t-1])
        J3[t+1] <- J2[t]*surv3[t] + J3[t]*survad[t]


        A[t+1] ~ dbin(survad[t], A[t])
        A[t+1] <- J3[t]*surv3[t] + A[t]*survad[t]

        # Lambda
        lambda[t+1] ~ dnorm(mean.lambda, tau.proc)
        N.est[t+1] <- N.est[t]*lambda[t]
        }

        # Population size  
        for (t in 1:TT){
            N[t] ~ dpois(N.est[t])
        }
}     
    ", fill = T)
sink()

# Parameters monitored
sp.params <- c("fec", "surv1", "surv2", "surv3", "survad", "lambda")

# MCMC settings
ni <- 200
nt <- 10
nb <- 100
nc <- 3

# Initial values
sp.inits <- function()list(mean.lambda = runif(1, 0, 1))

#Load all the data
sp.data <- list(N = Total, TT = length(Total), J0 = X0, J1 = X1, J2 = X2, J3     = X3)
library(R2jags)

hyena_model <- jags(sp.data, sp.inits, sp.params, "HyenaIPM_all.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)

К сожалению, при запуске кода я получаю следующую ошибку.

Error in jags.model(model.file, data = data, inits = init.values, n.chains = n.chains,  : 
  RUNTIME ERROR:
Index out of range for node J0

У кого-нибудь есть предложения, почему мы получаем эту ошибку? Не уверен, почему распределение было бы неправильным для J0.


person PendaFisi    schedule 05.06.2015    source источник


Ответы (1)


Это очень информативное сообщение об ошибке. Индекс для J0 равен t+1, который находится в диапазоне от 2+1 до TT+1, но J0 имеет длину TT. Поэтому, когда индекс равен TT+1, он выходит за пределы допустимого диапазона, поскольку он больше, чем TT.

person jaradniemi    schedule 07.06.2015