Complicated variance-covariance matrix not positive definite

Modeling
, ,
1

I am modeling a set of quantiles from a Gaussian mixture distribution, and it involves calculating the quantile function of a Gaussian mixture as well as populating a variance-covariance matrix with values that are a function of the density and quantile function. I’m using an algebraic solver to calculate the quantile function. I frequently, but not always, get the error that the VCOV matrix is not positive definite. This is surprising because all the entries should be positive. I try adding .000001 in case the numerical solutions give answers very close to 0, but it hasn’t fixed anything. Are there any suggestions on how to fix the issue? The model code is included below as well as how data might be generated.

functions{
      
  vector GMInv_CDF(vector Q, data real p, vector pi, vector mu, 
                   vector sigma, int N) {
    vector[1] z;
    real cdf;
    // this encodes the cdf from which we want to sample
    cdf = pi[1]*exp(normal_lcdf(Q | mu[1], sigma[1])) +
          pi[2]*exp(normal_lcdf(Q | mu[2], sigma[2])) +
          pi[3]*exp(normal_lcdf(Q | mu[3], sigma[3])) //+
          //pi[4]*exp(normal_lcdf(Q | mu[4], sigma[4]))
          ;
    z[1] = cdf - p;
    return z;
  }
  
  
  real GM_PDF(vector mu, vector sigma, vector pi, real y) {
    
    real density;
    
    density = pi[1]*(exp(normal_lpdf(y | mu[1], sigma[1]))) +
              pi[2]*(exp(normal_lpdf(y | mu[2], sigma[2]))) +
              pi[3]*(exp(normal_lpdf(y | mu[3], sigma[3]))) +
              //pi[4]*exp(normal_lpdf(y | mu[4], sigma[4])) +
              .00001;
              
    // density = log_sum_exp({log(pi[1]) + normal_lpdf(y | mu[1], sigma[1]),
    //                        log(pi[2]) + normal_lpdf(y | mu[2], sigma[2]),
    //                        log(pi[3]) + normal_lpdf(y | mu[3], sigma[3]),
    //                        log(pi[4]) + normal_lpdf(y | mu[4], sigma[4]),
    //                        log(pi[5]) + normal_lpdf(y | mu[5], sigma[5])});
    
    return density;
  }
  
  
}

data {
  int<lower=0> N;
  vector[N] Q;
  vector<lower=0, upper=1>[N] p;
  int<lower=1> n_components;
  real m;
  real<lower=0> c; 
  real<lower=0> sv; // sd prior sd
  real<lower=0> nv; // n prior sd
}

transformed data{
  real scaling_step = 1e-5;
  real rel_tol = 1e-6;
  real f_tol = 1;
  int max_steps = 3000;
  vector[1] y_guess = [0.5]';
}

parameters {
  ordered[n_components] mus;
  vector<lower=0>[n_components] sigmas;
  real<lower=0> n;
  vector<lower=0,upper=1>[n_components] pi;

}

transformed parameters {
  vector[N] Qi;
  cov_matrix[N] q_var;
  vector[n_components] pit;
  

  pit = softmax(pi);
  
  for (i in 1:N) {
    // // vector[1] Q_guess = [Q[i]]';
    Qi[i] = solve_powell_tol(GMInv_CDF, y_guess,
                             rel_tol, f_tol, max_steps,
                             p[i], pit, mus, sigmas, N)[1];
                             
    // Qi[i] = solve_newton_tol(GMInv_CDF, y_guess,
    //                          scaling_step, f_tol, max_steps,
    //                          p[i], pit, mus, sigmas, N)[1];
  }
  
  for (i in 1:N) {
    for (j in 1:i) {
      
      q_var[i, j] = ((fmin(p[i],p[j])*(1 - fmax(p[i],p[j]))) / 
                    (n*GM_PDF(mus, sigmas, pit, Qi[i])*
                     GM_PDF(mus, sigmas, pit, Qi[j]) + .000001)) +
                     .000001;
      
      q_var[j, i] = q_var[i, j];
      
    }
  }
  
  // q_var_cholesky = cholesky_decompose(q_var);

}

model {
  pi ~ normal(4,5);
  mus ~ normal(4, 5);
  sigmas ~ normal(0,7);
  n ~ normal(0, nv);

  Q ~ multi_normal(Qi, q_var);
  // Q ~ multi_normal_cholesky(Qi, q_var_cholesky);
}

I am treating the data I’m working with as though it came from a Gaussian mixture, but if the data was generated from a simple normal in R it could be as here


probs <- seq(0.05, .95, by = 0.05)
samp <-  rnorm(500)
quantiles <- quantile(samp, probs = probs)

dat <- list(
    N = length(quantiles),
    n = size,
    Q = quantiles,
    p = probs,
    n_components = 3,
    m = 2,
    c = 3,
    sv = 3,
    nv = 3000,
    pv = 3
  )
2

First, I would suggest adding a test to the code and printing when you find values you don’t expect. Such as

if (q_var[i, j] < 0) print("q_var[", i, ",", j, "] = ", q_var[i, j]);

This could be a problem with our solver’s tolerances.

Being symmetric with positive entries isn’t enough to ensure positive definiteness. You can, on the other hand, take a symmetric matrix and add to the diagonal until it’s positive definite. You can also get much fancier and do something like an eigendecomposition, then push the eigenvalues to where they’re positive, then put everything back together, but I wouldn’t recommend going that far for something like this.

Was there some other reason you have to believe the matrix you’re constructing will be positive definite?

P.S. I’d suggest pulling out reusable code, which can help avoid cut-and-paste errors. For example, you can define GMInv_CDF in terms of GM_PDF:

vector GMInv_CDF(vector Q, data real p, vector pi, vector mu, 
                   vector sigma, int N) {
  return GM_PDF(mu, sigma, pi, Q) - p;
}
3

Thanks for the tips! I’ll try them out and let you know how it goes.