Error with inverse Wishart prior

[dunajevas]

i'm trying to fit basic Bayesian mixture model for bivariate normal and uniform distributions. Basically when I increase degrees of freedom parameter (wishart_df) I get this warning (mostly during warm up stage, but if I increase degrees of freedom it becomes quite frequent):

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Error in function stan::prob::inv_wishart_log(d): LDLT_Factor of random variableunderlying matrix is not positive definite. LDLT_Factor of random variablelast conditional variance is -2.84217e-14.
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

I saw this error in https://github.com/stan-dev/stan/issues/662, but my variance matrix is correcly specified as cov_matrix type. Is there anything that might be wrong with my model? Here is my code

library(rstan) 
library(MASS)

mix_model <- 
  '
  data {
    int<lower = 1> N; // number of data points
    matrix[N, 2] xy; // observations
    cov_matrix[2] sigma_for_prior; //sigma for wishart prior
    int<lower = 1> wishart_df;
  }
  parameters {
    simplex[2] theta; // mixing proportions
    cov_matrix[2] sigma;
    vector<lower = 0, upper = 1>[2] mu;
  }
  model {
    real ps[2]; // temp for log component densities
    sigma ~ inv_wishart(wishart_df, sigma_for_prior);
    for (n in 1:N) {
      ps[1] <- log(theta[1])
        + multi_normal_log(xy[n], mu, sigma);
      ps[2] <- log(theta[2])
        + uniform_log(xy[n, 1], 0, 1)
        + uniform_log(xy[n, 2], 0, 1);
      increment_log_prob(log_sum_exp(ps));
    }
  }'

#(variance of unform distribution on [0,1/8])
variance <- 1/12 * (1/8)^2
correlation <- matrix(c(1, 0.3, 0.3, 1), ncol=2)
sigma_gen <- variance*correlation

xy <- rbind(mvrnorm(n = 50, mu = c(0.5, 0.5), Sigma=sigma_gen),
            matrix(runif(100, 0, 1), ncol=2))

plot(xy)
mix_data <- list(N = 100,
                   xy = xy,
                   sigma_for_prior = sigma_gen,
                   wishart_df = 2)
fit <- stan(model_code = mix_model, data = mix_data,
             iter = 1000, chains = 4)
fit
#check covariance matrix
fitext <- extract(fit)
sigma_gen
colMeans(fitext$sigma[,,1])
colMeans(fitext$sigma[,,2])

[bob-carpenter]

Are you running RStan 2.4?

If you're running into the problem mostly during warmup, you should be fine. Stan has to figure out the right scale to make jumps and before it does that, the arithmetic instability causes these errors in poorly conditioned models (and that's just about anything in high dimensions [and 50 counts as high here] trying to factor a covariance matrix [which multi_normal and multi_normal_prec do]). I wouldn't expect so many errors in 2 dimensions, though I don't see anything wrong with your model.

From a user point of view, it seriously reduces the usability of Stan. The problem is that using covariance matrices and factoring them (as required for the covariance matrix parameterization of multinormal) is very imprecise arithmetically.

I thought we'd relaxed the error checking in 2.4, but apparently not enough. Ben may have some suggestions on how to rewrite your models.

I'm closing this issue and referring you to #256.

[bgoodri]

For a 2x2 covariance matrix, it less likely to run into numerical errors if you specify two standard deviations and a correlation in the parameters block and construct the matrix in the model block. But then in order to use an inverse Wishart prior, you would have to do some stuff with the Jacobian. There is an example of this for the Wishart (not inverse Wishart at)

https://raw.githubusercontent.com/stan-dev/stan/develop/src/models/basic_distributions/wishart2x2.stan

But that is probably overkill; just ignore the informational message.