Distance sampling for point transect data with data argumentation, translated from Sec. 8.3.4 of "Applied Hierarchical Modeling in Ecology, vol. 1"

distance_sampling-4.stan

data {
  int<lower = 1> N_ind;                     // Number of individuals
  int<lower = 1> N_z;                       // Number of augment observed data
  int<lower = 0, upper = 1> Y[N_ind + N_z]; // Augumented inds. have y=0 by
                                            // definition
  int<lower = 1> N_D;
  real B;
  int <lower = 1, upper = N_D> D_class[N_ind];
  real<lower = 0, upper = B> Delta;           // Width of distance bins
  vector<lower = 0, upper = B>[N_D] Midpt;    // Interval mid-points
}
transformed data {
  vector<lower = 0, upper = 1>[N_D] pi = 2 * Midpt / square(B) * Delta;
  simplex[N_D] pi_probs = pi / sum(pi);
}
parameters {
  real<lower = 0> sigma;
  real<lower = 0, upper = 1> psi;
}
transformed parameters {
  vector<lower = 0, upper = 1>[N_D] p = exp(-(Midpt .* Midpt)
                                            / (2 * square(sigma)));
}
model {
  sigma ~ normal(0, 10);
  // Y == 1
  for (i in 1:N_ind) {
    target += bernoulli_lpmf(1 | psi)
              + bernoulli_lpmf(1 | p[D_class[i]])
              + categorical_lpmf(D_class[i] | pi_probs);
  }
  // Y == 0
  for (i in (N_ind + 1):(N_ind + N_z)) {
    real lp[N_D + 1];
    lp[1] = bernoulli_lpmf(0 | psi);
    for (g in 1:N_D) {
      lp[g + 1] = bernoulli_lpmf(1 | psi)
                  + bernoulli_lpmf(0 | p[g])
                  + categorical_lpmf(g | pi_probs);
    }
    target += log_sum_exp(lp);
  }
}
generated quantities {
  int<lower = 0, upper = N_ind + N_z> N = N_ind;
  real D;
  for (i in (N_ind + 1):(N_ind + N_z)) {
    real lp;
    real lp2[N_D];
    for (g in 1:N_D) {
      lp2[g] = bernoulli_lpmf(0 | p[g]) + categorical_lpmf(g | pi_probs);
    }
    lp = bernoulli_lpmf(1 | psi) + log_sum_exp(lp2)
         - log_sum_exp(bernoulli_lpmf(1 | psi) + log_sum_exp(lp2),
                       bernoulli_lpmf(0 | psi));
    N = N + bernoulli_rng(exp(lp));
  }
  D = N / (pi() * square(B));
}