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

きのうのモデルですが、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 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] = poisson_lpmf(n[1] - u | theta_1)
               + poisson_lpmf(n[2] - u | theta_1)
               + poisson_lpmf(u | theta_0);
    return log_sum_exp(s);
  }
}

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);
}

タグ：STAn