How to use matrix normal distribution in Stan?

Modeling
, ,
#1

kind regards,

I have a response variable with the following dimensions:

data {
int N; // number of observations.
int K; // number of locations.
int R; // number of sensors at each location.

matrix[R, K] y[N];
}

I want to model y using a matrix normal distribution but don’t know how to specify it in Stan. The definition of the matrix normal distribution is the following:

p(\mathbf{X}\mid\mathbf{M}, \mathbf{U}, \mathbf{V}) = \frac{\exp\left( -\frac{1}{2} \, \mathrm{tr}\left[ \mathbf{V}^{-1} (\mathbf{X} - \mathbf{M})^{T} \mathbf{U}^{-1} (\mathbf{X} - \mathbf{M}) \right] \right)}{(2\pi)^{np/2} |\mathbf{V}|^{n/2} |\mathbf{U}|^{p/2}}

\mathbf{X} (p \times n)
\mathbf{M} (p \times n)
\mathbf{U} (p \times p)
\mathbf{V} (n \times n)

source

1 Like
#2

A few options.

  1. You can code it yourself.
  2. Use the relation between a matrix-variate normal and a multivariate normal,
\text{MatN}(X \mid M,U,V) = MVN(vec(X) \mid M,V⊗U).

But stan doesn’t have kronecker product either. So you’d have to code that up.

  1. Quick google search and I also found this attempt in Stan.

  2. Long term option. Wait for it to come to Stan. There is a matrix normal in Stan-math that has not been exposed to Stan. Maybe some of the devs can weigh in on when that will make its way in @bbbales2 @mcol

4 Likes
#3

If you compute the log density manually then you can increment the target variable in Stan to make it part of the log density Stan is sampling.

I feel like I’m not saying that clearly, but in code that is:

real matrix_normal_log_density = ...; // a bunch of calculations
target += matrix_normal_log_density;

Indeed it hasn’t, and this is kind of a difficult function to compute so we should probably do something. I’d like to have the non-precision version as well. I made a stanc3 issue.

3 Likes
#4

I was able to get model running using your suggestion, thanks.

2 Likes