[Stan] 2回測定のN-mixture model, bivariate Poisson distribution [統計]
ある生物の個体数を測定するとします。調査箇所数がRで、調査箇所iの個体数Niは、平均λのポアソン分布にしたがうとします。各調査地でT回の観測をおこないます。検出確率をpとすると、調査箇所iのt回目の観測数はbinominal(Ni, p)にしたがいます。すなわち
Ni ~ Poisson(λ)
Yi,t ~ Binomial(Ni, p)
こうしたデータからλとpを推定します。
Stanでは離散パラメーターがあつかえないので、Nを周辺化してモデリングしていましたが、multivariate Poisson分布でモデリングできるというので、そちらで実装してみました。ただし、ここではT=2とします。
参考にしたところ
- Dennis et al. (2015) Computational aspects of N-mixture models. Biometrics 71:237--246. DOI:10.1111/biom.12246
- Stan users mailing list https://groups.google.com/forum/#!topic/stan-users/9mMsp1oB69g
Stanコード(Nmix2.stan)
/** * N-mixture model with 2 replicated observations * * References * Dennis et al. (2015) Computational aspects of N-mixture models. * Biometrics 71:237--246. DOI:10.1111/biom.12246 * Stan users mailing list * https://groups.google.com/forum/#!topic/stan-users/9mMsp1oB69g */ functions { /** * Returns log likelihood of N-mixture model * with 2 replicated observations * * @param n Number of observed individuals * @param lambda Poisson mean of population size * @param p Detection probability */ real n_mixture2_lpmf(int[] n, real lambda, real p) { real lp; real s[min(n) + 1]; real theta_1 = lambda * p * (1 - p); real theta_0 = lambda * p * p; for (u in 0:min(n)) s[u + 1] = lchoose(n[1], u) + lchoose(n[2], u) + lgamma(u + 1) + u * (log(theta_0) - 2 * log(theta_1)); lp = -(2 * theta_1 + theta_0) + (n[1] + n[2]) * log(theta_1) - (lgamma(n[1] + 1) + lgamma(n[2] + 1)) + log_sum_exp(s); return lp; } } data { int R; // Number of records int Y[R, 2]; // Observed number of individuals } parameters { real<lower=0> lambda; real<lower=0,upper=1> p; } model { for (i in 1:R) Y[i] ~ n_mixture2(lambda, p); }
Rコード
library(rstan) rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores()) set.seed(123) R <- 160 T <- 2 lambda <- 4 p <- 0.7 N <- rpois(R, lambda) Y <- structure(rep(0, R * T), dim = c(R, T)) for (i in 1:R) Y[i, ] <- rbinom(T, N[i], p) stan_data <- list(R = R, Y = Y) fit <- stan("Nmix2.stan", data = stan_data) print(fit)
結果
> print(fit) Inference for Stan model: Nmix2. 4 chains, each with iter=2000; warmup=1000; thin=1; post-warmup draws per chain=1000, total post-warmup draws=4000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat lambda 4.30 0.01 0.31 3.78 4.08 4.26 4.48 5.00 730 1 p 0.66 0.00 0.04 0.57 0.64 0.67 0.69 0.74 753 1 lp__ -553.99 0.03 1.03 -556.68 -554.38 -553.66 -553.27 -552.99 1362 1 Samples were drawn using NUTS(diag_e) at Thu Jan 26 22:29:23 2017. For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).
うまく推定できているようです。
ただし、T=2のときは、1回目のみで観測された個体数、2回目のみで観測された個体数、1回目2回目両方で観測された個体数と場合分けするのですが、Tがおおきくなると、この場合分けの数がどんどんふえそうな気がします。
タグ:STAn
コメント 0