So-net無料ブログ作成
検索選択

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に相当する部分を周辺化消去しています。

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 = binomial_rng(N_obs + N_aug, omega);
}

実行する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.32   24.61   27.04   28.47   29.99   33.86  1035 1.01
omega    0.79    0.00  0.09    0.61    0.73    0.79    0.85    0.97  1122 1.01
N      219.05    0.75 25.89  167.00  201.00  219.00  237.00  270.00  1178 1.01
lp__  -122.07    0.04  1.13 -125.20 -122.47 -121.71 -121.28 -121.00   988 1.00

Samples were drawn using NUTS(diag_e) at Sat Dec 10 10:17:51 2016.
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).

タグ:STAn
nice!(2)  コメント(0)  トラックバック(0) 
共通テーマ:日記・雑感

nice! 2

コメント 0

コメントを書く

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

Facebook コメント

トラックバック 0