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

distance_sampling-2.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
  real B;                                   // Strip half-width
                                            // Larger than max observed distance
  int<lower = 1> N_D;                         // N intervals
  int<lower = 1, upper = N_D> D_class[N_ind]; // Distance class
  real<lower = 0, upper = B> Delta;           // Width of distance bins
  vector<lower = 0, upper = B>[N_D] Midpt;    // Interval mid-points
}
transformed data {
  simplex[N_D] pi = rep_vector(Delta / B, N_D); // Prob. of x in each interval
}
parameters {
  real<lower = 0> sigma;          // Half-normal scale
  real<lower = 0, upper = 1> psi; // DA parameter
}
transformed parameters {
  vector<lower = 0, upper = 1>[N_D] p = exp(-(Midpt .* Midpt) / (2 * square(sigma)));
}
model {
  sigma ~ normal(0, 1000);
  // 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);
  }
  // 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);
    }
    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);
    }
    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 / 60.0;
}