使用JAGS純粹的出生過程推斷

Question

我想用JAGS來推斷（隨機）純出生過程中的出生率。

在化學的語言，這種模式是相當於反應：X-> 2X與速率的α* X（也可以被看作是一個鏈式反應的模型）

這是R代碼我用於生成過程（在固定時間）以及用於對參數α進行推斷的鋸齒代碼。

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 

我試圖像把阿爾法* Y [I-1]在一個新的變量不同的東西（比如，拉姆達[我]]或更改新[i-1]的調用新[i]，但沒有任何工作。任何想法爲什麼這是失敗的？另一個更聰明的方法來做到這一點？

預先感謝您。

Answer 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說法。