発見率をくみこんだ状態空間モデル

Kéry and Schaub (2011) Bayesian Population Analysis using WinBUGS: A hierarchical perspective にある状態空間モデルの例を参考に、発見率を推定するモデルをつくってみた。

rjagsとggplot2を使用。

library(rjags)
library(ggplot2)

データを生成。

set.seed(1234)
n.t <- 50                      # 観察回数
N.lat <- rep(50, n.t)          # 真の個体数
p <- 0.7                       # 発見確率
N.obs <- rbinom(n.t, N.lat, p) # 観察個体数

plt <- ggplot(data.frame(t = 1:n.t, N.lat = N.lat, N.obs = N.obs))
plt + geom_line(aes(x = t, y = N.lat)) +
  geom_line(aes(x = t, y = N.obs), colour = "red") +
  xlab("Time") + ylab("Population") +
  theme_classic(12, "Helvetica")

黒線が「真の個体数」、赤線が観測された個体数

モデル。

model.bugs <- "
var
  N,
  y[N],
  y_hat[N],
  lambda[N],
  p, tau, sigma;
model {
  lambda[1] ~ dnorm(0, 1.0E-4);
  for (t in 2:N) {
    lambda[t] ~ dnorm(lambda[t - 1], tau);
  }
  for (t in 1:N) {
    y_hat[t] <- trunc(exp(lambda[t]));
    y[t] ~ dbin(p, y_hat[t]);
  }
  ## priors
  p ~ dbeta(2, 2);
  sigma ~ dunif(0, 100);
  tau <- 1 / (sigma * sigma);
}"
write(model.bugs, file = "ks51.bug.txt")

初期値を設定して、JAGSを実行する。収束は非常におそい。

inits <- list()
inits[[1]] <- list(p = 0.9, sigma = 1, lambda = rep(log(max(N.obs) + 1), n.t))
inits[[2]] <- list(p = 0.7, sigma = 3, lambda = rep(log(max(N.obs) + 1), n.t))
inits[[3]] <- list(p = 0.8, sigma = 5, lambda = rep(log(max(N.obs) + 1), n.t))

model <- jags.model("ks51.bug.txt", data = list(N = n.t, y = N.obs),
                    inits = inits, n.chains = 3, n.adapt = 100000)
samp <- coda.samples(model, variable.names = c("y_hat", "sigma", "p"),
                     n.iter = 3000000, thin = 3000)

plt + geom_line(aes(x = t, y = N.lat)) +
  geom_line(aes(x = t, y = N.obs), colour = "red") +
  xlab("Time") + ylab("Population") +
  theme_classic(12, "Helvetica") +
  geom_line(aes(x = t, y = summary(samp)$statistics[3:(n.t + 2)]),
            colour = "blue")

結果

青線が、推定された「真の個体数」

> summary(samp[, c("p", "sigma")])

Iterations = 103000:3100000
Thinning interval = 3000 
Number of chains = 3 
Sample size per chain = 1000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean       SD  Naive SE Time-series SE
p     0.754660 0.054874 0.0010019      0.0024738
sigma 0.008927 0.008417 0.0001537      0.0001597

2. Quantiles for each variable:

           2.5%      25%      50%     75%   97.5%
p     0.6277649 0.725689 0.761888 0.79332 0.83910
sigma 0.0003017 0.003044 0.006592 0.01196 0.03154

タグ：jags R