Hierarchical Distance Sampling model using data augmentation, translated from Sec. 8.5.2 of "Advance Hierarchical Modeling in Ecology, vol. 1"

hds1.stan

data {
  int<lower = 1> N_sites;
  int<lower = 0> N_ind;
  int<lower = 0> N_z;
  real<lower = 0> B;
  vector[N_sites] Habitat;
  vector[N_sites] Wind;
  int<lower = 0, upper = 1> Y[N_ind + N_z];
  vector<lower = 0, upper = B>[N_ind] D;
  int<lower = 1, upper = N_sites> Site[N_ind];
}
transformed data {
  real area = N_sites * 1 * 2 * B;   // Unit length == 1, half-width = B
}
parameters {
  real beta0;  // Intercept of lambda-habitat regression
  real beta1;  // Slope of log(lambda) on habitat
  real alpha0; // Intercept of log(sigma) (half-normal scale)
  real alpha1; // Slope of log(sigma) on wind
  vector<lower = 0, upper = B>[N_z] d_new;
}
transformed parameters {
  vector[N_sites] lambda = exp(beta0 + beta1 * Habitat);
  vector[N_sites] sigma = exp(alpha0 + alpha1 * Wind);
  simplex[N_sites] site_probs = lambda / sum(lambda);
  real<lower = 0, upper = 1> psi = sum(lambda) / (N_ind + N_z);
}
model{
  beta0 ~ normal(0, 10);
  beta1 ~ normal(0, 10);
  alpha0 ~ normal(0, 10);
  alpha1 ~ normal(0, 10);

  d_new ~ uniform(0, B);

  // Y == 1
  for (i in 1:N_ind) {
    real p = exp(-square(D[i]) / (2 * square(sigma[Site[i]])));
    target += bernoulli_lpmf(1 | psi)
              + bernoulli_lpmf(1 | p)
              + categorical_lpmf(Site[i] | site_probs);
  }
  // Y == 0
  for (i in (N_ind + 1):(N_ind + N_z)) {
    vector[N_sites] lp;
    for (s in 1:N_sites) {
      real z[2];
      real p = exp(-square(d_new[i - N_ind]) / (2 * square(sigma[s])));
      z[1] = bernoulli_lpmf(0 | psi);
      z[2] = bernoulli_lpmf(1 | psi)
             + bernoulli_lpmf(0 | p);
      lp[s] = categorical_lpmf(s | site_probs) + log_sum_exp(z);
    }
    target += log_sum_exp(lp);
  }
}
generated quantities {
  int<lower = 0, upper = N_ind + N_z> N_total = N_ind;
  real dens;
  for (i in (N_ind + 1):(N_ind + N_z)) {
    int s = categorical_rng(site_probs);
    real p = exp(-square(d_new[i - N_ind]) / (2 * square(sigma[s])));
    real lp1 = bernoulli_lpmf(1 | psi) + bernoulli_lpmf(0 | p);
    real lp = lp1 - log_sum_exp(bernoulli_lpmf(0 | psi), lp1);

    N_total = N_total + bernoulli_rng(exp(lp));
  }
  dens = N_total / area;
}