Lkj_corr_lpdf: Correlation matrix is not positive definite

General
7

You can try using the onion method to construct the matrix. This works for sizes much larger than 30.

functions {
  array[] vector create_shapes(int K, real eta) {
    real alpha = eta + (K - 2) / 2.0;
    array[2] vector[K - 1] shape;
    
    shape[1, 1] = alpha;
    shape[2, 1] = alpha;
    for (k in 2 : (K - 1)) {
      alpha -= 0.5;
      shape[1, k] = k / 2.0;
      shape[2, k] = alpha;
    }
    
    return shape;
  }
  
  matrix lkj_onion(int K, row_vector l, vector R2, data array[] vector shape) {
    matrix[K, K] L = rep_matrix(0, K, K); // cholesky_factor corr matrix
    {
      int start = 1;
      int end = 2;
      
      L[1, 1] = 1.0;
      L[2, 1] = 2.0 * R2[1] - 1.0;
      L[2, 2] = sqrt(1.0 - square(L[2, 1]));
      for (k in 2 : (K - 1)) {
        int kp1 = k + 1;
        row_vector[k] l_row = segment(l, start, k);
        real scale = sqrt(R2[k] / dot_self(l_row));
        L[kp1, 1 : k] = l_row[1 : k] * scale;
        L[kp1, kp1] = sqrt(1.0 - R2[k]);
        start = end + 1;
        end = start + k - 1;
      }
    }
    return L;
  }
}
data {
  int<lower=0> K; // dim of corr mat
  real<lower=0> eta; // lkj param
}
transformed data {
  array[2] vector[K - 1] shapes = create_shapes(K, eta);
}
parameters {
  row_vector[choose(K, 2) - 1] l; // do NOT init with 0 for all elements
  vector<lower=0, upper=1>[K - 1] R2; // first element is not really a R^2 but is on (0,1)  
}
transformed parameters {
  cholesky_factor_corr[K] L = lkj_onion(K, l, R2, shapes);
}
model {
  target += std_normal_lpdf(l);
  target += beta_lpdf(R2 | shapes[1], shapes[2]);
}
generated quantities {
  matrix[K, K] Sigma = multiply_lower_tri_self_transpose(L);
}