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<lower = 0, upper = B>[N_ind + N_z] x2 = append_row(X, x_new);
  vector<lower = 0, upper = 1>[N_ind + N_z] p = exp(-(x2 .* x2)
                                                    / (2 * square(sigma)));
}
model {
  X ~ uniform(0, B);
  x_new ~ uniform(0, B);
  for (n in 1:N_ind) {
    target += bernoulli_lpmf(1 | psi) + bernoulli_lpmf(1 | p[n]);
  }
  for (n in (N_ind + 1):(N_ind + N_z)) {
    target += log_sum_exp(bernoulli_lpmf(0 | psi),
                          bernoulli_lpmf(1 | psi)
                          + bernoulli_lpmf(0 | p[n]));
  }
}
generated quantities {
  int<lower = 0, upper = N_ind + N_z> N = N_ind;
  real D;
  for (n in (N_ind + 1):(N_ind + N_z)) {
    real zp = psi * (1 - p[n]) / (psi * (1 - p[n]) + (1 - psi));
    N = N + bernoulli_rng(zp);
  }
  D = N / 60.0;
}