Fitting normal-inverse-wishart

Specific Domains Probabilistic programming
1

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.

data = [ ([0, 0], [1 2; 2 1]),
         ([0.1, 0.1], [1.1 2.1; 2.1 1.1]),
         ([0.05, 0.05], [1.05 2.05; 2.05 1.05]),
         ([0, 0], [1 2; 2 1]),
         ([0, 0], [1 2; 2 1])
         ]
2

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.)

1 Like
3

In your data here, your covariance matrices aren’t positive semidefinite though.

4

Sorry, I just wanted to show the structure of the data and made up values

1 Like
5

On the off chance this is still needed, you can fit parameters of the NIW distribution using the procedure described here [2405.16088] Estimating the normal-inverse-Wishart distribution. I have a python/JAX implementation and I will probably get around to a Julia implementation at some point, but I could do it sooner if it is useful for anyone.