Marginalisation over correlation matrices with the non-centred parameterisation of the multivariate normal

Modeling
1

Another covariance matrix marginalisation question from me, I’m afraid. I have n covariance matrices which have been converted to correlation matrices, inputted as data. I want to marginalise over them to get means (assuming each correlation matrix has an implicit uniform prior). In the centered parameterisation, this is easy, and various people on this board have been very helpful in showing me how to do it. The prior for the means, is (with “mean” and “scale” being parameters):

lp = rep_vector(log_unif, n);
for (i in 1:n) {
    lp[i] += multi_normal_cholesky_lpdf(mean | 0, diag_pre_multiply(scale, cholesky_of_correlation_matrix[i]));
}
target += log_sum_exp(lp);
###prior on scale###

I can’t work out how to convert this to the non-centred version though. If you were not marginalising (i.e. you had a single known correlation matrix), for a 0 mean MVN it would be:

mean = scale * (cholesky_of_correlation_matrix * mean_tilde);
mean_tilde ~ normal(0,1);
###prior on scale###

The logical way to marginalise would be to do the above for each matrix, but that ends up with n sets of means, which is clearly not what’s required. If anyone could point me in the direction of something explaining how to do this, I would greatly appreciate it!

2

I think this might kind of defeat the purpose of the reparameterization. For a given ‘mean’, each correlation matrix induces a different transformation from mean_tilde, so you’d need to have a different mean_tilde for each correlation matrix in order to wind up with the same mean. I think the most logical way to do that would be to sample the means and then back transform them to mean_tilde… which is essentially what the original, centered parameterization does within multi_normal_cholesky_lpdf.

3

I suppose it could potentially be helpful to come up with some kind of ‘central’ correlation matrix that relates a single set of raw parameters to a single set of means. Not sure if I would expect it to behave nicely but could be worth a try.

transformed_data {
matrix new_choleskies[n];
for (i in 1:n) {
new_choleskies[i] = mdivide_left_tri_low(cholesky_of_correlation_matrix[i]), cholesky_of_average_correlation_matrix);
}
}

model {
mean = scale * (cholesky_of_average_correlation_matrix * mean_tilde);
for (i in 1:n) {
lp[i] += std_normal_lpdf(new_choleskies[i] * mean_tilde);
}
target += log_sum_exp(lp);
}