2

我正在尝试使用 JAGS 来推断（随机）纯出生过程中的出生率。

用化学的语言来说，这个模型等价于反应：X->2X with rate alpha*X（也可以看成是链式反应的模型）

这是我用来生成过程的 R 代码（在固定时间）和用于推断参数 alpha 的 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-1] 更改对 New[i] 的调用，但没有任何效果。知道为什么会失败吗？另一种更聪明的方法是什么？

先感谢您。

4

1 回答 1

1

另一种解决方案是更改模拟数据的方式并在模型中使用链接功能。

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)

您可以在下面看到该模型也正确检索了参数（alpha = 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

于 2016-09-29T14:01:10.250 回答