What library (and method) is recommended for taking a sequence of mean-vectors and covariance matrices and fitting a normal-inverse-wishart distribution?

For example, my data is a vector of tuples where the first element is the tuple is an array (of mean values) and the second element is a covariance matrix. I would like to pass the data to a function and have it output the parameters of a normal-inverse-wishart distribution that fits it the best.

I would try Turing.jl (or Stan – there are Julia bindings), and consider a model that does something like the following:

\mu_i \sim N(\mu, \Sigma), for i = 1...k \Sigma_i \sim IW(\nu, \Psi), for i = 1...k \nu \sim Truncated(Normal(0, 10), p-1, Inf) \Psi \sim IW(p-1, I)

where p is the dimension of the elements in your vector, and I is a p-dimensional identity matrix.

In this model, you are estimating the means of the means (allowing for covariance), and estimating the parameters (\nu and \Psi) of the inverse Wishart that best fit the observed covariance matrices. And because we actually want to estimate those parameters, we also need hyperpriors for them.

I’ll warn you, however, I’ve had very bad luck getting inverse Wishart distributions to work in Turing for anything with more than a very small number of dimensions.

I don’t know if this will work, but it’s where I would probably start if I were doing it myself.

(Edit: I didn’t write it down, but you also need priors for \mu and \Sigma, of course.)