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

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

Stan: distance sampling (2)　[統計] [編集]

前回のつづきです。今回は、距離データが実際の測定値ではなく、クラスわけした値となっている場合です。

前回のデータを加工します。観測幅を幅deltaで分割したnD段階のクラスとして、距離データdclassが得られるとしています。

delta <- 10
xg <- seq(0, width, delta)[-1];
dclass <- x %/% delta + 1
nD <- length(xg)
midpt <- xg - delta / 2

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

> dclass
 [1] 6 2 1 2 1 1 2 4 4 1 2 5 1 6 2 2 3 1 5 2 2 3 5 2 2 6 3 1 3 3 4 2 2 4 4 4 6 1
[39] 1 2 4 1 1 3 2 1 2 3 7 4 3 5 1 4 2 6 1 3 2 2 6 1 3 3 3 1 2 3 1 4 5 2 1 2 5 1
[77] 2 4

Stanコード(ds2.stan)です。前回と同様にData augmentationを使用しています。Data augmentationで追加したゼロデータ部分で、dclassに相当する部分を周辺化消去しています。（2018-09-16 追記: Nを推定する部分を修正しました。）

data {
  int<lower=0> N_obs;                 // Number of observed individuals
  int<lower=0> N_aug;                 // Number of augmentated zeros
  int<lower=0> ND;                    // Number of bins
  real<lower=0> B;                    // Transect width
  real<lower=0,upper=B> D;            // Bin width
  vector<lower=0,upper=B>[ND] M;      // Mid point for each cell
  int<lower=0,upper=ND> DC[N_obs];    // Observed distance class
}

transformed data {
  simplex[ND] p_i; // Probability of x in each interval

  p_i = rep_vector(D / B, ND);
}

parameters {
  real<lower=0,upper=10^3> sigma; // Scale parameter of detection func.
  real<lower=0,upper=1> omega;    // Inclusion parameter
}

transformed parameters {
  vector[ND] p;   // Detection probability

  p = exp(-(M .* M) / (2 * sigma^2));
}

model {
  sigma ~ uniform(0, 10^3);
  // Observed data
  for (i in 1:N_obs) {
    1 ~ bernoulli(omega);
    DC[i] ~ categorical(p_i);
    1 ~ bernoulli(p[DC[i]]);
  }
  // Augmented zero data
  for (i in 1:N_aug) {
    vector[ND] lp;
    real lp1;
    real lp2;

    for (j in 1:ND) {
      lp1 = bernoulli_lpmf(0 | omega);
      lp2 = bernoulli_lpmf(1 | omega)
          + bernoulli_lpmf(0 | p[j]);
      lp[j] = categorical_lpmf(j | p_i)
            + log_sum_exp(lp1, lp2);
    }
    target += log_sum_exp(lp);
  }
}

generated quantities {
  int N = N_obs;

  for (i in 1:N_aug) {
    real lp;
    real lpnz[ND];
    for (j in 1:ND) {
      lpnz[j] = bernoulli_lpmf(0 | p[j]) + categorical_lpmf(j | p_i);
    }
    lp = bernoulli_lpmf(1 | omega) + log_sum_exp(lpnz)
         - log_sum_exp(bernoulli_lpmf(1 | omega) + log_sum_exp(lpnz),
                       bernoulli_lpmf(0 | omega));
    N = N + bernoulli_rng(exp(lp));
  }
}

実行するRコードです。

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

stan_data <- list(N_obs = length(dclass), N_aug = 200,
                  ND = nD, B = 100, D = delta,
                  M = midpt, DC = dclass)
params = c("sigma", "omega", "N", "lp__")
n.chains <- 4
inits <- lapply(1:n.chains, function(i) {
    list(sigma = runif(1, 5, 20),
         omega = runif(1, 0.3, 0.9))
    })
fit <- stan("ds2.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: ds2.
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.63    0.07  2.37   24.64   26.92   28.42   30.05   34.00  1194 1.00
omega    0.79    0.00  0.09    0.61    0.72    0.79    0.86    0.97  1235 1.00
N      219.39    0.70 25.24  172.00  201.00  219.00  238.00  268.03  1292 1.00
lp__  -122.11    0.04  1.10 -125.06 -122.54 -121.79 -121.33 -121.01   786 1.01

タグ：STAn