data { int N; int D; vector[D] x[N]; vector[N] y; } transformed data { vector[N] mu; real delta = 1e-9; mu = rep_vector(0, N); } parameters { real alpha; vector[D] rho; real sigma; } transformed parameters { matrix[N, N] L_K; { // K = gp_exp_quad_cov(x, magnitude, length_scale); // K = add_diag(K, square(sigma)); // K = cov_exp_quad(X, magnitude, length_scale); // K = K + diag_matrix(rep_vector(square(sigma), N)); matrix[N, N] K; real sq_alpha = square(alpha); for (i in 1:(N-1)) { K[i, i] = sq_alpha + delta; for (j in (i + 1):N) { K[i, j] = sq_alpha * exp(-0.5 * dot_self((x[i] - x[j]) ./ rho)); K[j, i] = K[i, j]; } } K[N, N] = sq_alpha + delta; L_K = cholesky_decompose(K); } } model { alpha ~ normal(0, 3); rho ~ inv_gamma(5, 5); sigma ~ normal(0, 1); y ~ multi_normal_cholesky(mu, L_K); }