Currently, I am trying to fit a lognormal model for one of the parameters in the Stan program. For one-dimension case, the easy way is just using the built-in Stan function as:

\beta \sim lognormal(\mu, \sigma)

However, currently my parameter \beta is a vector of length 3, where \beta = (\beta_1, \beta_2, \beta_3) and I would like to jointly model the three components as multivariate lognormal distribution with some mean vector and covariance matrix (maybe treat them as hyperparameter in the model). May I ask that does Stan provide such built-in functions for multivariate lognormal distribution?

Usually when we include multivariate random effects in a model, we are modeling covariances between multiple vectors, not multiple elements of a single vector. If \beta_1, \beta_2, and \beta_3 are all elements of a single random-effect vector, then there is just the variance and no additional covariance. If, on the other hand, \beta_1, \beta_2, and \beta_3 all represent single-element random effect vectors, then how do you expect to be able to estimate the underlying means, variances, and covariances? If the means, variances, and covariances are tightly constrained a priori then this could be a useful approach, but thatâ€™s not what you seem to be driving at when you write

This implies that beta has a multi_lognormal(mu, Sigma) distribution, where Sigma = L_Sigma * L_Sigma'. With this code, you can use beta just as you would had it been given a multivariate lognormal directly.

Multivariate logistic normal works the same way, only with inv_logit() in place of exp().

[edit: added he covariance matrix as data and its cholesky factor as a transformâ€”you could also take this to be a parameter, for which you might consider the LKJ prior on a cholesky factor of correlation matrices scaled by a vector of scales with their own independent prior]