I don’t know of a Bayesian NMF implementation in Stan.
For PSIS-LOO and bad Pareto k, I typically avoid the conditional likelihood and use marginal instead. If you wanted to use conditional likelihood, I think the recommendation would be to set up your own cross-validation. Some more details are here (especially see the FAQ).
For parameterizing the multivariate normal via mean and covariance matrix, I think the code below will do it. Mu and Sigma are the multivariate normal parameters; xbar and S are the observed mean and covariance matrix; N is number of observations; p is the dimension of the multivariate normal:
out = -.5 * N * ( sum( inverse_spd(Sigma) .* (S + (xbar - Mu) * (xbar - Mu)') ) + log_determinant(Sigma) + p * log(2 * pi()) );