SSブログ

[Stan] 2回測定のN-mixture model, bivariate Poisson distribution [統計]

ある生物の個体数を測定するとします。調査箇所数がRで、調査箇所iの個体数Niは、平均λのポアソン分布にしたがうとします。各調査地でT回の観測をおこないます。検出確率をpとすると、調査箇所it回目の観測数はbinominal(Ni, p)にしたがいます。すなわち
Ni ~ Poisson(λ)
Yi,t ~ Binomial(Ni, p)
こうしたデータからλpを推定します。

Stanでは離散パラメーターがあつかえないので、Nを周辺化してモデリングしていましたが、multivariate Poisson分布でモデリングできるというので、そちらで実装してみました。ただし、ここではT=2とします。

参考にしたところ

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
nice!(2)  コメント(0)  トラックバック(1) 
共通テーマ:日記・雑感

nice! 2

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

Facebook コメント

トラックバック 1