Those are some really nice implementations of truncated MVN functions.
However, I am still confused about what transformation/generation to do where. I do need Beta not just as a generated quantity, I need it as a parameter. I should have clarified that. Let me provide some mock code to illustrate the setup:
functions{
real myLLfun(... matrix Beta ...) //lower level log-likelihood function for reduce_sum
real multi_normal_cholesky_truncated_lpdf()
vector multi_normal_cholesky_truncated_rng()
}
data{
X, mu_0, mu_sig, //data, prior specs
S, // number of subjects
Q, // number of parameters
lb[Q,S],ub[Q,S],lb_ind[Q,S],ub_ind[Q,S]
}
parameters {
vector[Q] mu; // upper MVN: mean vector of ind betas
vector<lower=0>[Q] tau;
cholesky_factor_corr[Q] L_Omega;
matrix<lower=0,upper=1>[Q, S] Z; // individual z, unif(0,1)
}
transformed parameters {
//chol-of-covmat
matrix[Q,Q] L = diag_pre_multiply(tau, L_Omega);
//lower level para
matrix[Q, S] Beta;
{
for (r in 1:S) {
Beta[,r] = multi_normal_cholesky_truncated_rng(Z[,r], mu, L, lb[,r],ub[,r],lb_ind[,r],ub_ind[,r]);
}
}
}
model {
//2nd stage
tau ~ cauchy(0, 2.5);
L_Omega ~ lkj_corr_cholesky(2);
mu ~ normal(mu_0, mu_sig);
//1st stage
to_vector(Z) ~ uniform(); // implicit
for (r in 1:S) {
target+= multi_normal_cholesky_truncated_lpdf(Z[,r], mu, L, lb[,r],ub[,r],lb_ind[,r],ub_ind[,r])
}
//logll
target += reduce_sum(myLLfun, ..., Beta, X, ...);
}
The _rng function is really just a transformation going from unconstrained (uniform-distributed) space to constrained space (the MVN-trun distributed space). The _lpdf function is what I need for the logpost. Somehow I think I am missing something here. Or, rather, I am double counting something: I have Z and transformed Z in here, but I think the only thing to add in model{} would be the Jacobian.