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);
}
```