Hi,
I’ve got a model that uses the multi_normal Stan function as below to model a set of data with the variance covariance matrix V. I’m trying to generate and store the pointwise log-likelihood in the generated quantities block as recommended by the loo package for later model comparison. My issue is the multi_normal_lpdf takes in a vectors and a covariance matrix and outputs a single real, and I am unsure how to get individual pointwise log-likelihoods for the observations using this approach.
I’m probably missing something obvious, and thanks for any assistance.
MG
model{
...
mu = X*beta;
Y ~ multi_normal(mu ,V);
}
generated quantities{
...
log_lik = multi_normal_lpdf(Y | mu, V);}
}
As far as I know, there is not a function to obtain the pointwise log-likelihood in this way. I believe you need a loop in generated quantities to compute the log-likelihood separately for each observation in Y. Under this approach, it is faster to obtain the Cholesky decomposition of V once, and then use multi_normal_cholesky_lpdf()each time through the loop.
Thanks! But multi_normal_cholesky_lpdf() returns a real log-likelihood value, but takes a vector - how would the loop compute the log-likelihood separately for each value when all Y values are used in each iteration?
Best,
MG
It might help to see how you declare Y in the data block. For N observations and multivariate normal of dimension p, I typically do
vector[p] Y[N];
You can then do the loop in generated quantities:
for (i in 1:N) {
log_lik[i] = multi_normal_cholesky(Y[i] | mu, V_chol);
}
My data declaration above might be “illegal” using the new Stan array declarations, but I’m not immediately sure because I typically use rstan. You might need something like vector Y[N, p];
instead (would be good to verify that, though).
I think there’s a confusion here. When the outcome is multivariate, I think @mgrab is expecting to be able to write down a point-wise log likelihood such that each element of the multivariate outcome makes its own point-wise contribution. But it doesn’t work this way. The multivariate outcome isn’t factorizable in this way, because the likelihood associated with the ith element of the outcome is not conditionally (on the parameters) independent of the values taken by all the other elements of the outcome.
Note however, that in a model where you have a multivariate normal observation-level random effect, and then some other conditionally independent response distribution, now you can again talk about the point-wise likelihood per element. However, you will not likely be able to use PSIS-loo in this situation, because you will have a random effect parameter per observation, which will induce too much flexibility in the model for the PSIS approximation to the leave-one-out posterior to be reliable.
Thanks @edm and @jsocolar, appreciate the help. Question - ok, given PSIS-loo is not appropriate here, what is an equivalent way of quantifying predictive accuracy given a multivariate outcome?
I also found the link below, and checked out the original article by Bürkner et al. (2020), which refers to “LOO-CV for multivariate normal models” - is this a way forward?
Thanks,
MG