Trivariate model with a binary margin

albert_chib_style_trivariate.stan

data {
  int N;
  int N_neg;
  int N_pos;
  int P;
  matrix[N, P] X;
  matrix[N, 3] Y_raw; // first column is 0s and 1s
  int pos_neg[N_neg];
  int pos_pos[N_pos];
}
parameters {
  vector[3] alpha;
  vector<lower = 0>[2] tau_raw;
  cholesky_factor_corr[3] L_Omega;
  matrix[3, P] Beta;
  vector<lower = 0>[N_pos] xi_pos;
  vector<upper = 0>[N_neg] xi_neg;
}
transformed parameters {
  vector[3] tau = append_row(rep_vector(1.0, 1), tau_raw);
  vector[N] xi;
  for(i in 1:N_neg) {
    xi[pos_neg[i]] = xi_neg[i];
  }
  for(i in 1:N_pos) {
    xi[pos_pos[i]] = xi_pos[i];
  }
}
model {
  matrix[N, 3] Y = append_col(xi, Y_raw[,2:3]);
  matrix[N, 3] z = (Y - rep_matrix(alpha',N) - X * Beta') * inverse(diag_matrix(tau) * L_Omega);
  
  // priors
  alpha ~ normal(0, 2);
  to_vector(Beta) ~ normal(0, 1);
  tau_raw ~ cauchy(0, 1);
  L_Omega ~ lkj_corr_cholesky(3.0);
  // likelihood
  to_vector(z) ~ normal(0, 1);
}