LDLT_Factor of random variable is not positive definite

Modeling
,
2

Unfortunately, large covariance matrices have initialization problems. I believe if you try this with N < 30 then it will work. Once you get over a dimension of ~30 then you will start to see this error.

You can use the Using the onion method in occupancy model hierarchical model - #2 by spinkney.

This will give you a cholesky factor that you can multiply by the transpose of itself to get a corr mat. Then multiply by the vector of variances to get a cov mat.

Or you can use the wishart_cholesky directly

functions {
matrix onion_cholesky_lp(
		       int N,
                       real eta,
		       vector r_2,
		       row_vector l_omega) {       
  {
    matrix[n_pred, n_pred] L_Omega = rep_matrix(0, N, N);
    real alpha = eta + (N - 2.0) / 2.0;
    vector[N - 1] shape1;
    vector[N - 1] shape2;
 
     shape1[1] = alpha;
     shape2[1] = alpha;

    int start = 1;
    int end = 2;

    L_Omega[1,1] = 1.0;
    L_Omega[2,1] = 2.0 * r_2[1] - 1.0;
    L_Omega[2,2] = sqrt(1.0 - square(L_Omega[2,1]));
    for (k in 2:N - 1) {
      int kp1 = k + 1;
      row_vector[k] l_row = segment(l_omega, start, k);
      real scale = sqrt(r_2[k] / dot_self(l_row));
      L_Omega[kp1, 1:k] = l_row[1:k] * scale;
      L_Omega[kp1, kp1] = sqrt(1.0 - r_2[k]);
      
     alpha -= 0.5;
    shape1[k] = k / 2.0;
    shape2[k] = alpha;

      start = end + 1;
      end = start + k - 1;
    }

    target += beta_lpdf(r_2, shape1, shape2);

    return L_Omega;
  }
}
}
data {
  int<lower=1> N;
  array[N] vector[2] index;
  vector[N] y_obs;
  real<lower=0> eta;
}


parameters{
  real<lower=0> sigma;
  real<lower=0> length_scale;
 row_vector[choose(N, 2) - 1]  l_omega_p;
vector<lower = 0,upper = 1>[N - 1] R2;
 // cov_matrix[N] KW;
}

transformed parameters{
  matrix[N, N] L_KW = onion_cholesky(N, eta, R2, l_omega_p);
  matrix[N, N] L_K = cholesky_decompose(gp_exponential_cov(index, sigma, length_scale) + diag_matrix(rep_vector(1e-9, N)));
}

model{
  sigma ~ gamma(1, 1);
  length_scale ~ gamma(1, 1);
  L_KW ~ wishart_cholesky(N, L_K);
  
  y_obs ~ multi_normal_cholesky(rep_vector(0, N), L_KW);
}