Ensure correlation matrix is positive semi-definite?

Modeling
6

You can find a near PSD correlation matrix to the given priors. PSD matrices are really rare compared to all square matrices, especially so as the dimension increases.

You need to somehow get PSD by construction. Stan’s corr_matrix declaration in the parameter block will construct a valid PSD matrix. By using this I can put your priors on the lower elements of this matrix to find a corr matrix near the given priors. In the 3 x 3 case you mention, it’s likely that the priors will be close to the corr matrix entries. However, there’s still a large region of values that are not permissible. Visualized as a floating samosa in 3d:
image
(from my blog post Correlation Cats Part I – The Stan Blog )

Meaning, that there are some priors which will be non-permissible for PSD and this non-permissibility increases as the dimensions of the correlation matrix increase. You can mitigate this somewhat by choosing a sigma on the soft normal constraint to be larger (ie set sigma > 0.01).

functions {
  vector lower_elements(matrix M, int tri_size){
    int n = rows(M);
    int counter = 1;
    vector[tri_size] out;
    for (i in 2:n){
      for (j in 1:(i - 1)) {
        out[counter] = M[i, j];
        counter += 1;
      }
    }
    return out;
  }
}
data { 
  int<lower=1> Q; //  dimension of the v-cov matrices
  int<lower=1> I; //  no. of subjects
} 
transformed data {
  int<lower=0> off_diag_Q = (Q * (Q - 1))/2;
}
parameters { 
  vector<lower=0, upper=1>[off_diag_Q] r_raw[I]; //values to draw from Beta distribution
  real<lower=0> hyper_alpha; // hyperparameter on the Beta prior
  real<lower=0> hyper_beta; //hyperparameter on the Beta prior
  corr_matrix[Q] R[I];
}
model {
  //hyperparameters on the S_is 
  hyper_alpha ~ exponential(1);
  hyper_beta ~ exponential(1);

  for(i in 1:I) {
    r_raw[i] ~ beta(hyper_alpha, hyper_beta);
    lower_elements(R[i], off_diag_Q) ~ normal(2 * r_raw[i] - 1, 0.01);
  }
}
8 Likes