How to increase machine_precision

We would like to use an inverse Wishart prior distribution for a covariance matrix 40x40 , we get the error:

Ccov is not positive definite.

We suspect it is a machine precision problem, is it possible to increase machine precision? Since it looks like the maximum precision allowed is e(-16)

1e-16 is 64-bit floating point (Wikipedia has a nice table). It seems Eigen support long double.

Let’s consider a very simple example, it doesn’t work if P is greater than 30.

data {
int<lower=1> P; // number of variables
int<lower=1> N; // number of observations
matrix[N,P] x; // observations
vector[P] mu; // expected value
}

parameters{

cov_matrix[P] Ccov;
}

transformed parameters {

cov_matrix[P] identitaP;

for (i in 1:P){
for (j in 1:P)
{ if (i==j) { identitaP[i,j]=1;}
else identitaP[i,j]=0;
}
}
}

model {
for (j in 1:N)
x[j]~ multi_normal(mu,Ccov);
Ccov~inv_wishart(P+2,identitaP);
}

I got this message
Exception: multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite. last conditional variance is 0. (in ‘model10b4245857c9_programma_wishart’ at line 28)

Exception: multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite. last conditional variance is -2.22045e-016. (in ‘model10b4245857c9_programma_wishart’ at line 28)

[1] "Error in sampler$call_sampler(args_list[[i]]) : “
[2] " Exception: multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite. last conditional variance is -2.22045e-016. (in ‘model10b4245857c9_programma_wishart’ at line 28)”
[1] “error occurred during calling the sampler; sampling not done”

That isn’t very simple (for a 64bit computer) because P is greater than 30. That is one reason why people tend to work with Cholesky factors of covariance matrices in Stan rather than the covariance matrices themselves. We don’t have a density function for a Cholesky factor of a covariance matrix that is distributed inverse Wishart, but that is mostly because the inverse Wishart is not a good distribution for modeling anyway.

If you really want to assume that a covariance matrix is distributed inverse Wishart, it is better to work with the Cholesky factor of a precision matrix, which can be decomposed

and then construct the covariance matrix in the generated quantities block.

But it is better to just abandon the inverse Wishart distribution entirely and decompose a covariance matrix into a standard deviations and (the Cholesky factor of) a correlation matrix, which is discussed in the Stan User Manual and hundreds of threads.

1 Like