# Fitting normal-inverse-wishart

**URL:** <https://discourse.julialang.org/t/fitting-normal-inverse-wishart/73273>\
**Category:** Probabilistic Programming\
**Created:** [December 17, 2021, 5:12pm UTC](https://discourse.julialang.org/t/fitting-normal-inverse-wishart/73273 "2021-12-17T17:12:45Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Paul\_McVay](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paul_mcvay/32/4307_2.png) [@Paul\_McVay](https://discourse.julialang.org/u/Paul_McVay)\
**Post date:** [December 17, 2021, 5:12pm UTC](https://discourse.julialang.org/t/fitting-normal-inverse-wishart/73273/1 "2021-12-17T17:12:45Z")

</div>

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.

```julia
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])
         ]

```

---

<div class="post-metadata">

**Author:** ![opera\_malenky](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/opera_malenky/32/8213_2.png) [@opera\_malenky](https://discourse.julialang.org/u/opera_malenky)\
**Post date:** [December 17, 2021, 5:49pm UTC](https://discourse.julialang.org/t/fitting-normal-inverse-wishart/73273/2 "2021-12-17T17:49:26Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![trahflow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trahflow/32/30585_2.png) [@trahflow](https://discourse.julialang.org/u/trahflow)\
**Post date:** [December 17, 2021, 6:09pm UTC](https://discourse.julialang.org/t/fitting-normal-inverse-wishart/73273/3 "2021-12-17T18:09:10Z")

</div>

> [@Paul\_McVay](#):
>
> ```julia
> 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])
> ]
> 
> ```

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

---

<div class="post-metadata">

**Author:** ![Paul\_McVay](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paul_mcvay/32/4307_2.png) [@Paul\_McVay](https://discourse.julialang.org/u/Paul_McVay)\
**Post date:** [December 17, 2021, 6:18pm UTC](https://discourse.julialang.org/t/fitting-normal-inverse-wishart/73273/4 "2021-12-17T18:18:41Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jso](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jso/32/209698_2.png) [@jso](https://discourse.julialang.org/u/jso)\
**Post date:** [June 4, 2024, 10:06am UTC](https://discourse.julialang.org/t/fitting-normal-inverse-wishart/73273/5 "2024-06-04T10:06:38Z")

</div>

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](https://arxiv.org/abs/2405.16088). 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.