functions { vector gp_pred_exp_quad_rng(vector[] x2, vector y, vector[] x1, real m, real l, real sig, real sigma) { int N1 = rows(y); int N2 = size(x2); vector[N2] f; { matrix[N1, N1] K = add_diag(gp_exp_quad_cov(x1, m, l) + gp_dot_prod_cov(x1, sig), sigma); matrix[N1, N1] L_K = cholesky_decompose(K); vector[N1] L_K_div_y = mdivide_left_tri_low(L_K, y); vector[N1] K_div_y = mdivide_right_tri_low(L_K_div_y', L_K)'; matrix[N1, N2] k_x1_x2 = gp_exp_quad_cov(x1, x2, m, l); f = (k_x1_x2' * K_div_y); } return f; } vector gp_pred_dot_prod_rng(vector[] x2, vector y, vector[] x1, real m, real l, real sig, real sigma) { int N1 = rows(y); int N2 = size(x2); vector[N2] f; { matrix[N1, N1] K = add_diag(gp_exp_quad_cov(x1, m, l) + gp_dot_prod_cov(x1, sig), sigma); matrix[N1, N1] L_K = cholesky_decompose(K); vector[N1] L_K_div_y = mdivide_left_tri_low(L_K, y); vector[N1] K_div_y = mdivide_right_tri_low(L_K_div_y', L_K)'; matrix[N1, N2] k_x1_x2 = gp_dot_prod_cov(x1, x2, sig); f = (k_x1_x2' * K_div_y); } return f; } vector gp_pred_rng(vector[] x2, vector y, vector[] x1, real m, real l, real sig, real sigma) { int N1 = rows(y); int N2 = size(x2); vector[N2] f; { matrix[N1, N1] K = add_diag(gp_exp_quad_cov(x1, m, l) + gp_dot_prod_cov(x1, sig), sigma); matrix[N1, N1] L_K = cholesky_decompose(K); vector[N1] L_K_div_y = mdivide_left_tri_low(L_K, y); vector[N1] K_div_y = mdivide_right_tri_low(L_K_div_y', L_K)'; matrix[N1, N2] k_x1_x2 = gp_exp_quad_cov(x1, x2, m, l) + gp_dot_prod_cov(x1, x2, sig); f = (k_x1_x2' * K_div_y); } return f; } } data { int N; int D; vector[D] x[N]; vector[N] y; int N_pred; vector[D] x_pred[N_pred]; } transformed data { vector[N] mu = rep_vector(0, N); } parameters { real l; real m; real sig0; real sigma; } transformed parameters { matrix[N, N] L_K; { matrix[N, N] K = gp_exp_quad_cov(x, m, l) + gp_dot_prod_cov(x, sig0); K = add_diag(K, sigma); L_K = cholesky_decompose(K); } } model { m ~ normal(0, 1); l ~ normal(0, 1); sig0 ~ normal(0, 1); sigma ~ normal(0, 1); y ~ multi_normal_cholesky(mu, L_K); } generated quantities { vector[N_pred] f1_pred = gp_pred_exp_quad_rng(x_pred, y, x, m, l, sig0, sigma); vector[N_pred] f2_pred = gp_pred_dot_prod_rng(x_pred, y, x, m, l, sig0, sigma); vector[N_pred] f3_pred = gp_pred_rng(x_pred, y, x, m, l, sig0, sigma); }