Stan model for occupancy model with covariates, translated from Section 10.4 of "Advanced Hierarchical Modeling in Ecology, vol. 1" by Kéry and Royle.

occ_cov2.stan

data {
  int<lower = 1> M;
  int<lower = 1> J;
  int<lower = 0, upper = 1> Y[M, J];
  vector[M] VegHt;
  matrix[M, J] Wind;
  vector[M] XvegHt;
  vector[M] Xwind;
}
transformed data {
  int present[M];
  int C;
  for (m in 1:M) {
    present[m] = sum(Y[m, ]) > 0 ? 1: 0;
  }
  C = sum(present);
}
parameters {
  real mean_psi;
  real mean_p;
  real alpha1;
  real beta1;
}
transformed parameters {
 real alpha0 = logit(mean_p);
 real beta0 = logit(mean_psi);
 vector[M] logit_psi = logit(mean_psi) + beta1 * VegHt;
 matrix[M, J] logit_p = logit(mean_p) + alpha1 * Wind;
}
model {
  for (m in 1:M) {
    if (present[m] == 0) {
      real lp[2];
      lp[1] = bernoulli_logit_lpmf(0 | logit_psi[m]);
      lp[2] = bernoulli_logit_lpmf(1 | logit_psi[m])
              + bernoulli_logit_lpmf(0 | logit_p[m, ]);
      target += log_sum_exp(lp);
    } else {
      target += bernoulli_logit_lpmf(1 | logit_psi[m])
                + bernoulli_logit_lpmf(Y[m, ] | logit_p[m, ]);
    }
  }
  alpha1 ~ normal(0, 10);
  beta1 ~ normal(0, 10);
}
generated quantities {
  int<lower = C, upper = M> N_occ = 0;
  real<lower = 0, upper = 1> psi_fs;
  vector[100] psi_pred = inv_logit(beta0 + beta1 * XvegHt);
  vector[100] p_pred = inv_logit(alpha0 + alpha1 * Xwind);
  for (m in 1:M) {
    if (present[m] == 1) {
      N_occ = N_occ + 1;
    } else {
      // Bern(1 | psi) Bern(0 | p) /
      //     (Bern(1 | psi) Bern(0 | p) + Bern(0 | psi))
      real psi = inv_logit(logit_psi[m]);
      vector[J] p = inv_logit(logit_p[m, ]');
      real log_p = bernoulli_lpmf(1 | psi) + bernoulli_lpmf(0 | p)
                   - log_sum_exp(bernoulli_lpmf(1 | psi) + bernoulli_lpmf(0 | p),
                                 bernoulli_lpmf(0 | psi));
      N_occ = N_occ + bernoulli_rng(exp(log_p));
    }
  }
  psi_fs = (N_occ * 1.0) / M;
}