Prior for a covariance matrix

Modeling
, , ,
#1

Hi Everyone,

I am stuck in a fix.
I want to sample a covariance matrix \Sigma_e.
So, I am using Inv-Wishart distribution to sample \Sigma_e using scale matrix \Lambda and degree of freedom v. And then use this covariance matrix \Sigma_e as a parameter in Multivariate normal distribution.
Everything works fine if the dimension of the covariance matrix is small e.g., ten or less, but if the dimension of the matrix is large, e.g., 11 or more.
The model does not fit properly and shows fitting errors.
Can anyone help me understand what is happening and how could I fix it?

sigma_err   ~    inv_wishart(sigma_error_dof, diag_matrix(rep_vector(sigma_error_variance,N)));
count       ~    multi_normal(count_hat, sigma_err);

this is the error that i am getting:

Exception: inv_wishart_lpdf: LDLT_Factor of random variable is not positive definite.  last conditional variance is -4.44089e-16.
#2

How did you declare sigma_err? Try

cov_matrix[K] sigma_err;

If this does not help, we recommend using the LKJ distribution and it’s cholesky decomposition. See the Stan manual.
It’s much faster and numerically stable.

#3

This is how I am modelling the process,

model_code = """
    data {
        matrix[24, 6]         X[75];
        vector[75]            crime_count[24];
    }

    parameters {
        vector[75]            alpha;
        vector[6]             beta[75];
        cov_matrix[75]        sigma_err;
    }

    transformed parameters {
        vector[75]            crime_count_hat[24];
        for(j in 1:75) {
            for(t in 1:24) {
                crime_count_hat[t, j] = alpha[j] + dot_product(X[j, t], beta[j]);
            }
        }
    }
    
    model {
        alpha       ~    normal(0, 100);
        beta        ~    multi_normal(rep_vector(0, 6), diag_matrix(rep_vector(50,6)));
        sigma_err   ~    inv_wishart(76, diag_matrix(rep_vector(0.1, 75)));
        crime_count ~    multi_normal(crime_count_hat, sigma_err);
    }
    """
#4

I have tried passing different values to the scale matrix, instead of 0.1, for e.g. 100, 1000

But in each case, it gives same error.

#5

I know it’s not ideal. If you want to keep the inverse wishart, you may use the barlett decomposition of the wishart distribution and since it’s triangular the calc. of it’s inverse is cheap. There is a chapter in Stan’s manual too.

#6

Could you provide me with a link. It seem I am not able to find it.

If possible , can you help me with a example code for bartlett decomposition of the Wishart distribution?

#7

Pages 280ff, especially 282-283.

#8

Thank you.

#9

Is there any way we can put some kind of try-catch construct around inv-wishart sampling, so that when samples are not positive definite it just resamples until it finds one?

#10

I don’t think it’s the way to go.
You can provide your own initialization, which helps with high-dimensional LKJ sampling problems, which is the way recommended to sample a covariance matrix.
What you can try is to sample from a wishart distribution apply the cholesky decomposition and invert it.