Constrained multinormal

Modeling
3

Thanks @bgoodri, that works nicely. Based on your paper, I added the case with both lower and upper bounds

functions {
  vector[] make_stuff(vector mu, matrix L, vector lb, vector ub, vector u) {
    int K = rows(mu); vector[K] d; vector[K] z; vector[K] out[2];
    for (k in 1:K) {
      int km1 = k - 1;
      
      if (lb[k] == negative_infinity() && ub[k] == positive_infinity()) { // no bounds
        z[k] = inv_Phi(u[k]);
        d[k] = 1;
      } else {
        real v; 
        
        if (lb[k] != negative_infinity() && ub[k] != positive_infinity()) { // lower & upper bounds
          real z_star_lb = (lb[k] - (mu[k] + ((k > 1) ? L[k,1:km1] * head(z, km1) : 0))) / L[k,k];
          real z_star_ub = (ub[k] - (mu[k] + ((k > 1) ? L[k,1:km1] * head(z, km1) : 0))) / L[k,k];
          real u_star_lb = Phi(z_star_lb);
          real u_star_ub = Phi(z_star_ub);
          v = u_star_lb + (u_star_ub - u_star_lb) * u[k];
          d[k] = u_star_ub - u_star_lb;          
        } else if (ub[k] != positive_infinity()) { // upper bound
          real z_star = (ub[k] - (mu[k] + ((k > 1) ? L[k,1:km1] * head(z, km1) : 0))) / L[k,k];
          real u_star = Phi(z_star);
          v = u_star * u[k];
          d[k] = u_star;
        } else { // lower bound
          real z_star = (lb[k] - (mu[k] + ((k > 1) ? L[k,1:km1] * head(z, km1) : 0))) / L[k,k];
          real u_star = Phi(z_star);
          d[k] = 1 - u_star;
          v = u_star + d[k] * u[k];
        }
        z[k] = inv_Phi(v);
      }
    }
    out[1] = z;
    out[2] = d;
    return out;
  }
}
data {
  int<lower=2> K;             // number of dimensions
  vector[K] mu;
  vector[K] sigma;
  vector[K] lb;              // lower bounds, -Inf for no bound
  vector[K] ub;              // upper bounds,  Inf for no bound
  real rho;
}

transformed data {
  matrix[K, K] Omega;
  matrix[K, K] Sigma;
  cholesky_factor_cov[K,K] L;
  for (m in 1:K) for (n in 1:K) Omega[m, n] = (m==n) ? 1.0 : rho;
  Sigma = quad_form_diag(Omega, sigma);
  L = cholesky_decompose(Sigma);
}

parameters {
  vector<lower=0,upper=1>[K] u;
}
model {
  target += log(make_stuff(mu, L, lb, ub, u)[2]); // Jacobian adjustments
  // implicit: u ~ uniform(0,1)
}
generated quantities {
  vector[K] y = mu + L * make_stuff(mu, L, lb, ub, u)[1];
}