Чистый вывод процесса рождения с использованием JAGS

Я пытаюсь использовать JAGS, чтобы вывести уровень рождаемости в (стохастическом) чистом процессе рождения.

На языке химии эта модель эквивалентна реакции: Х->2Х со скоростью альфа*Х (также может рассматриваться как модель цепной реакции)

Это код R, который я использую для генерации процесса (в фиксированное время), и код jags для вывода альфа-параметра.

library(rjags)

y <- 1; # Starting number of "individuals"
N <- 25 # number of time samplings
alpha <- 0.2 # per-capita birth rate
# Generate the time series
for(i in 2:N) {
  y<- c(y,y[i-1]+rpois(1,alpha*y[i-1]))
};

# The jags code
model_string <- "model{
  for(i in 2:N) {
    New[i] ~ dpois(alpha*y[i-1])
    y[i] <- y[i-1] + New[i]
  }
  alpha ~ dunif(0, 2)
}"

# Create and run the jags model
model <- jags.model(textConnection(model_string), data = list(y = y,N = N), n.chains = 3, n.adapt= 10000)
update(model, 5000); # Burnin for 10000 samples
mcmc_samples <- coda.samples(model, variable.names=c("alpha"), n.iter=5000)

Когда я запускаю код, я получаю следующую ошибку:

Error in jags.model(textConnection(model_string), data = list(y = y, N = N),  : 
  RUNTIME ERROR:
Compilation error on line 4.
y[2] is a logical node and cannot be observed

Я пробовал разные вещи, например помещал alpha*y[i-1] в новую переменную (скажем, lambda[i]) или менял вызовы New[i] на New[i-1], но ничего не получалось. Любая идея, почему это терпит неудачу? Другой более умный способ сделать это?

Заранее спасибо.


r bayesian jags stochastic-process
person Donald Lockwood    schedule 29.09.2016    source источник
comment
Я нашел ответ в другом месте. Вы можете прочитать это здесь: sourceforge.net/p/mcmc-jags /обсуждение/610036/поток/6b159634   -  person Donald Lockwood    schedule 29.09.2016

Ответы (1)


arrow_upward
1
arrow_downward

Другим решением было бы изменить способ имитации данных и использовать функцию связи с моделью.

N <- 25 # number of time samplings
alpha <- 0.2 # log per-capita birth rate

# Generate the time series
steps <- 1:N # simulating 25 steps
log.y<- alpha*steps # the log-scale expected count
expected.y <- exp(log.y) # back to the real scale
y <- rpois(N, expected.y) # add Poisson noise to your expected.y

# The jags code
model_string <- "model{
  for(i in 1:N) {
    y[i] ~ dpois(lambda[i])
    log(lambda[i]) <- log.lambda[i]
    log.lambda[i] <- alpha * i
  }
  alpha ~ dunif(-10, 10)
}"

# Create and run the jags model
model <- jags.model(textConnection(model_string),inits = list(alpha = 1), data = list(y = y,N = N), n.chains = 1, n.adapt= 10000)
update(model, 5000); # Burnin for 10000 samples
mcmc_samples <- coda.samples(model, variable.names=c("alpha"), n.iter=5000)

Ниже вы можете видеть, что эта модель также правильно извлекает параметр (альфа = 0,2).

введите здесь описание изображения

Взяв экспоненту, вы получите уровень рождаемости (т.е. exp(0.2) = 1.22), или вы можете сделать это в рамках модели и отследить производный параметр, который является просто экспонентой alpha. Тогда модель будет:

model_string <- "model{
  for(i in 1:N) {
    y[i] ~ dpois(lambda[i])
  log(lambda[i]) <- log.lambda[i]
  log.lambda[i] <- alpha * i
  }
  alpha ~ dunif(-10, 10)
  birth.rate <- exp(alpha)
}"

И вы бы просто отследили birth.rate в аргументе variable.names.

person mfidino    schedule 29.09.2016
comment
Это определенно лучшее решение, чем я ожидал. Но я не совсем понимаю, почему лямбда-скорость является экспонентой alphai, а не просто alphai. - person Donald Lockwood; 30.09.2016
comment
Поскольку эта модель использует ссылку на журнал. Инверсией логарифмической связи является экспоненциальная, которая возвращает ее к шкале, которую легче интерпретировать (коэффициент). На странице Википедии, посвященной регрессии Пуассона, достаточно информации об этом. Так, например, если вы хотите смоделировать скорость увеличения 1,2 (увеличение на 20% от одного временного шага к другому), вы должны взять его журнал alpha <- log(1.2) и использовать его в моделировании. Эта модель также позволяет включать ковариаты, в то время как другое решение этого не позволяет. - person mfidino; 30.09.2016
comment
Теперь ясно. Спасибо еще раз. - person Donald Lockwood; 30.09.2016