Stan: distance sampling：Taglibro de H：SSブログ

	ブログをはじめるログイン

Stan: distance sampling　[統計] [編集]

Applied Hierarchical Modeling in Ecology vol.1の8章にあるDistance samplingデータを解析するStanモデルをつくってみました。

Applied Hierarchical Modeling in Ecology: Analysis of distribution, abundance and species richness in R and BUGS: Volume 1:Prelude and Static Models

作者: Marc Kery
出版社/メーカー: Academic Press
発売日: 2015/12/10
メディア: ハードカバー

ある直線上を移動しながら、両側一定幅内で観測された目的の生物の数をかぞえます。このとき、発見した生物との距離も記録しておきます。発見率は距離に応じて減少するとします。このようなデータから、発見されなかった個体もふくめた個体数を推定します。

データをシミュレートすると、以下のようになります。

set.seed(123)

## Half-normal detection function
g <- function(x, sigma) {
  exp(-x^2 / (2 * sigma^2))
}

N <- 200     # Number of individuals
sigma <- 30  # Scale parameter of half-normal detection function
width <- 100 # Half-width of the transect

xall <- runif(N, -width, width) # Distance
p <- g(xall, sigma)             # Detection probability
y <- rbinom(N, 1, p)            # Detection
x <- abs(xall[y == 1])          # Distance to detected individuals

このようなデータとなります。

> x
 [1] 57.6610271 18.2046156  5.6210976 10.2870029  8.6770529  9.3331688
 [7] 14.5266804 38.5606812 31.1411598  8.8132049 18.8284041 42.1680525
[13]  4.4408058 53.6748429 17.0907328 17.2551347 26.2309098  6.8075099
[19] 46.8054719 11.5599852 12.1895968 23.2060724 45.1232711 10.2967317
[25] 12.0336625 50.8950317 25.8442263  4.9366852 22.5542007 29.6404182
[31] 33.6111175 16.4706441 13.0214517 30.6203850 31.2967055 35.9253515
[37] 56.4588603  6.6441917  2.3010920 19.9977919 33.4352919  2.2773933
[43]  3.4195206 21.7469965 17.8620447  9.8569312 17.0966706 26.1022268
[49] 69.1595398 38.0014203 23.8512967 47.4155476  4.2271452 31.9676899
[55] 18.1050095 52.1800089  0.4945466 22.4181940 14.3870628 11.0463996
[61] 56.4018663  0.4599127 29.2190856 25.0572087 28.9109238  6.7375891
[67] 17.4507763 25.9946107  5.9671372 37.5103670 46.9964388 18.8686388
[73]  3.7420399 12.9180870 45.1666757  6.6934595 18.6334813 31.8460648

発見されなかった個体数を推定するため、Data augmentationをつかって、発見されなかった個体のデータを追加します。発見されなかった個体までの距離は欠測値としてあつかいます。以下がStanコード(ds1.stan)です。（2016-12-10 追記: Nを推定する部分のコードを修正しました。2018-09-16 追記: Nを推定する部分のコードをさらに修正しました。）

data {
  int<lower=0> N_obs; // Number of observed individuals
  int<lower=0> N_aug; // Number of augmentated zeros
  real<lower=0> B;    // Transect width
  vector<lower=0,upper=B>[N_obs] X; // Distance data
}

parameters {
  real<lower=0,upper=10^3> sigma; // Scale parameter of detection func.
  real<lower=0,upper=1> omega;    // Inclusion parameter
  vector<lower=0,upper=B>[N_aug] x_aug; // Unobserved distance data
}

transformed parameters {
  vector[N_obs+N_aug] p; // Detection probability

  for (i in 1:N_obs)
    p[i] = exp(-X[i]^2 / (2 * sigma^2));
  for (i in (N_obs + 1):(N_obs + N_aug))
    p[i] = exp(-x_aug[i - N_obs]^2 / (2 * sigma^2));
}

model {
  sigma ~ uniform(0, 10^3);
  x_aug ~ uniform(0, B);
  // Observed data
  1 ~ bernoulli(omega * p[1:N_obs]);
  // Augmentated zero data
  for (i in (N_obs + 1):(N_obs + N_aug)) {
    real lp[2];

    lp[1] = bernoulli_lpmf(0 | omega);
    lp[2] = bernoulli_lpmf(1 | omega)
          + bernoulli_lpmf(0 | p[i]);
    target += log_sum_exp(lp[1], lp[2]);
  }
}

generated quantities {
  int N = N_obs;
  
  for (i in 1:N_aug) {
    real lpnz =  bernoulli_lpmf(1 | omega) + bernoulli_lpmf(0 | p[N_obs + i]);
    real lp = lpnz - log_sum_exp(lpnz, bernoulli_lpmf(0 | omega));
    N = N + bernoulli_rng(exp(lp));
  }
}

Stanによるはてはめを実行するRコードです。Data augmentationで200個のゼロデータを追加することにしています。

## Stan
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

stan_data <- list(X = x, N_obs = length(x),
                  N_aug = 200, B = 100)
params = c("sigma", "omega", "N", "lp__")
n.chains <- 4
inits <- lapply(1:n.chains, function(i) {
    list(sigma = runif(1, 0, 10),
         omega = runif(1, 0, 1))
    })
fit <- stan("ds1.stan", data = stan_data,
            pars = params, init = inits,
            chains = n.chains, iter = 2000, warmup = 1000, thin = 1,
            seed = 1)
print(fit)

結果です。だいたい良い値を推定できました。

> print(fit)
Inference for Stan model: ds1.
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
sigma  28.79    0.04  2.35  24.81  27.14  28.58  30.30  33.97  4000 1.00
omega   0.78    0.00  0.09   0.61   0.72   0.78   0.85   0.96   993 1.00
N     218.54    0.78 24.41 171.00 202.00 218.00 235.00 268.00   981 1.00
lp__  431.08    0.63 15.45 400.03 420.76 431.56 441.19 461.32   596 1.01

タグ：STAn