Distance sampling model translated from Sec. 8.3.1 of "Advanced Hierarchical Modeling in Ecology, vol.1"

distance_sampling-1.stan

data {
  int<lower = 1> N_ind;                     // Number of individuals
  int<lower = 1> N_z;                       // Number of augment observed data
  vector<lower = 0>[N_ind] X;               // Observed distance
  int<lower = 0, upper = 1> Y[N_ind + N_z]; // Augumented inds. have y=0 by
                                            // definition
  real B;                                   // Strip half-width
                                            // Larger than max observed distance
}
parameters {
  real<lower = 0> sigma;          // Half-normal scale
  real<lower = 0, upper = 1> psi; // DA parameter
  vector<lower = 0, upper = B>[N_z] x_new;
}
transformed parameters {
  vector[N_ind + N_z] p;
  for (i in 1:N_ind) {
    p[i] = exp(-square(X[i]) / (2 * square(sigma)));
  }
  for (i in (N_ind + 1):(N_ind + N_z)) {
    p[i] = exp(-square(x_new[i - N_ind]) / (2 * square(sigma)));
  }
}
model {
  x_new ~ uniform(0, B);
  // Y == 1
  for (i in 1:N_ind) {
    target += bernoulli_lpmf(1 | psi) + bernoulli_lpmf(1 | p[i]);
  }
  // Y == 0
  for (i in (N_ind + 1):(N_ind + N_z)) {
    target += log_sum_exp(bernoulli_lpmf(0 | psi),
                          bernoulli_lpmf(1 | psi)
                          + bernoulli_lpmf(0 | p[i]));
  }
}
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 zp = psi * (1 - p[i]) / (psi * (1 - p[i]) + (1 - psi));
    N = N + bernoulli_rng(zp);
  }
  D = N / 60.0;
}