diff --git a/posterior_database/models/stan/birthday2_6_28.stan b/posterior_database/models/stan/birthday2_6_28.stan new file mode 100755 index 00000000..ef467be2 --- /dev/null +++ b/posterior_database/models/stan/birthday2_6_28.stan @@ -0,0 +1,273 @@ +functions { + vector gp_pred_exp_quad_rng(real[] x2, + vector y, real[] x1, + matrix k, + real magnitude, real length_scale, real sigma, // squared exp params + model noise + real jitter) { + // x2: test data + // x1: training data + // magnitude: magnitude of squared exponential kernel + // length_scale: length_scale of squared exponential kernel + // sigma: model noise + // jitter: for numerical stability for inverse, etc. + + int N1 = rows(y); + int N2 = size(x2); + vector[N2] f; + { + matrix[N1, N1] K = k;// + diag_matrix(rep_vector(square(sigma), N1)); + 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 = cov_exp_quad(x1, x2, magnitude, length_scale); + vector[N2] f_mu = (k_x1_x2' * K_div_y); + + matrix[N1, N2] v_pred = mdivide_left_tri_low(L_K, k_x1_x2); + matrix[N2, N2] cov_f = cov_exp_quad(x2, magnitude, length_scale) - v_pred' * v_pred + + diag_matrix(rep_vector(jitter, N2)); + vector[N2] temp; + matrix[N2, N2] cov_f_diag; + for (n2 in 1:N2) + temp[n2] = cov_f[n2, n2]; + cov_f_diag = diag_matrix(temp); + + f = multi_normal_rng(f_mu, cov_f_diag); + + } + return f; + } +/////////////////////////// +// periodic + squared exponential prediction, +// for kernels 3 and 4 in the sum + vector gp_pred_periodic_exp_quad_rng(real[] x2, + vector y, real[] x1, + matrix k, + real magnitude_1, real length_scale_1, + real magnitude_2, real length_scale_2, real period, + real sigma, real jitter) { + // x2: test data + // x1: training data + // magnitude_1: magnitude of squared exponential kernel + // length_scale_1 length scale of squared exponential kernel + // magnitude_2: magnitude of periodic kernel + // length_scale_2: length_scale of periodic kernel + // period: perio of periodic kernel + // sigma: model noise + // jitter: for numerical stability for inverse, etc. + + int N1 = rows(y); + int N2 = size(x2); + vector[N2] f; + { + matrix[N1, N1] K = k; // + diag_matrix(rep_vector(square(sigma), N1)); + 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_periodic_cov(x1, x2, magnitude_2, length_scale_2, period) .* + cov_exp_quad(x1, x2, magnitude_1, length_scale_1); + vector[N2] f_mu = (k_x1_x2' * K_div_y); + + matrix[N1, N2] v_pred = mdivide_left_tri_low(L_K, k_x1_x2); + matrix[N2, N2] cov_f = gp_periodic_cov(x2, magnitude_2, length_scale_2, period) .* + cov_exp_quad(x2, magnitude_1, length_scale_1) - v_pred' * v_pred // K_t - (K_tn) ( K_n + I * sigma^2)^-1 K_nt) + sigma^2 I + + diag_matrix(rep_vector(jitter, N2)); + vector[N2] temp; + matrix[N2, N2] cov_f_diag; + for (n2 in 1:N2) + temp[n2] = cov_f[n2, n2]; + cov_f_diag = diag_matrix(temp); + + f = multi_normal_rng(f_mu, cov_f_diag); + + } + return f; + } +/////////////////////////// + vector gp_pred_dot_prod_rng(real[] x2, + vector y, real[] x1, + matrix k, + real intercept, real sigma, // squared exp params + model noise + real jitter) { + // x2: test data + // x1: training data + // intercept: intercept term in dot product kernel, sigma0 in GPML + // sigma: model noise + // jitter: for numerical stability for inverse, etc. + + int N1 = rows(y); + int N2 = size(x2); + vector[N2] f; + { + matrix[N1, N1] K = k;// + diag_matrix(rep_vector(square(sigma), N1)); + 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, intercept); + vector[N2] f_mu = (k_x1_x2' * K_div_y); + + matrix[N1, N2] v_pred = mdivide_left_tri_low(L_K, k_x1_x2); + matrix[N2, N2] cov_f = gp_dot_prod_cov(x2, intercept) - v_pred' * v_pred + + diag_matrix(rep_vector(jitter, N2)); +// f = multi_normal_rng(f_mu, cov_f); + vector[N2] temp; + matrix[N2, N2] cov_f_diag; + for (n2 in 1:N2) + temp[n2] = cov_f[n2, n2]; + cov_f_diag = diag_matrix(temp); + + f = multi_normal_rng(f_mu, cov_f_diag); + } + return f; + } +} +/////////////////////////// +data { + int N; + int N_pred; + int M; + real x[N]; + vector[N] y; + + // for the weekends and special days + vector[M] I_ss[N]; // N x M matrix + vector[M] I_ws[N]; // "" "" + real x_pred[N_pred]; +} +transformed data { + vector[N] mu = rep_vector(0, N); +} +//////////////////////////// +parameters { + // f_1(t), squared exponential + real length_scale_1; + real sigma_1; + + // f_2(t), squared exponential + real length_scale_2; + real sigma_2; + + // f_3(t), weekly quasi-periodic + real length_scale_3_1; + real sigma_3_1; + real length_scale_3_2; + + // f_4(t), yearly smooth seasonal + real length_scale_4_1; + real sigma_4_1; + real length_scale_4_2; + + // f_5(t), indicators for weekends and special days + real sigma_5_1; + real sigma_5_2; + + // noise + real sigma; +} +transformed parameters { + // TODO change this: + matrix[N, N] k = cov_exp_quad(x, sigma_1, length_scale_1) + + cov_exp_quad(x, sigma_2, length_scale_2) + + gp_periodic_cov(x, sigma_3_1, length_scale_3_1, 7) .* + cov_exp_quad(x, 1.0, length_scale_3_2) + + gp_periodic_cov(x, sigma_4_1, length_scale_4_1, 365.25) .* + cov_exp_quad(x, 1.0, length_scale_4_2) + + gp_dot_prod_cov(I_ss, sigma_5_1) + + gp_dot_prod_cov(I_ws, sigma_5_2); + real sq_sigma = square(sigma); + for (n in 1:N) + k[n, n] = k[n, n] + 1e-12; +} + +model { + matrix[N, N] L_k; + L_k = cholesky_decompose(k); + + sigma ~ normal(0, 1); + + // smooth non-periodic component + length_scale_1 ~ lognormal(log(365), 1); + sigma_1 ~ normal(0, 1); +// sigma_1 ~ normal(.7, 4); // and also set inits as the MAP for this joint distribution? + // making sigma priors more diffuse tends to help, + // whereas making them more strict tends to make the model not initialize + + // faster changing non-periodic component + length_scale_2 ~ lognormal(log(10), 1); + sigma_2 ~ normal(0, 1); +// sigma_2 ~ normal(.4, 4); + + // 7 day period + length_scale_3_1 ~ lognormal(log(2), sqrt(2)); + sigma_3_1 ~ normal(0, 1); // ends up at 20 // decay?? + length_scale_3_2 ~ lognormal(log(20), 1); + + // 365.25 day period + length_scale_4_1 ~ lognormal(log(100), sqrt(2)); + sigma_4_1 ~ normal(0, 1); + length_scale_4_2 ~ lognormal(log(1000), sqrt(2)); + + sigma_5_1 ~ normal(0, 1); + sigma_5_2 ~ normal(0, 1); + + y ~ multi_normal_cholesky(mu, L_k); +} + +/////////////////////////// +generated quantities { + + // f1 predictive + vector[N_pred] f1_pred; + vector[N_pred] y1_pred; + // f2 predictive + vector[N_pred] f2_pred; + vector[N_pred] y2_pred; + // f3 predictive + vector[N_pred] f3_pred; + vector[N_pred] y3_pred; + // f4 predictive + vector[N_pred] f4_pred; + vector[N_pred] y4_pred; + // f5_1 predictive + vector[N_pred] f5_1_pred; + vector[N_pred] y5_1_pred; + // f5_2 predictive + vector[N_pred] f5_2_pred; + vector[N_pred] y5_2_pred; + + + // // f1 predictive + // f1_pred = gp_pred_exp_quad_rng(x_pred, y, x, k, sigma_1, length_scale_1, sigma, 1e-6); + + // // f2 predictive + // f2_pred = gp_pred_exp_quad_rng(x_pred, y, x, k, sigma_2, length_scale_2, sigma, 1e-6); + + // // f3 predictive + // f3_pred = gp_pred_periodic_exp_quad_rng(x_pred, y, x, k, + // sigma_3_1, length_scale_3_1, + // 1.0, length_scale_3_2, 7, + // sigma, 1e-6); + // // f4 predictive + // f4_pred = gp_pred_periodic_exp_quad_rng(x_pred, y, x, k, + // sigma_4_1, length_scale_4_1, + // 1.0, length_scale_4_2, 365.25, + // sigma, 1e-6); + // // f5_1 predictive + // f5_1_pred = gp_pred_dot_prod_rng(x_pred, y, x, k, sigma_5_1, sigma, 1e-6); + + // // f5_2 predictive + // f5_2_pred = gp_pred_dot_prod_rng(x_pred, y, x, k, sigma_5_2, sigma, 1e-6); + + + for (n in 1:N_pred) + { +// y1_pred[n] = normal_rng(f1_pred[n], sigma); +// y2_pred[n] = normal_rng(f2_pred[n], sigma); + // y3_pred[n] = normal_rng(f3_pred[n], sigma); + // y4_pred[n] = normal_rng(f3_pred[n], sigma); + // y5_1_pred[n] = normal_rng(f5_1_pred[n], sigma); + // y5_2_pred[n] = normal_rng(f5_2_pred[n], sigma); + } +} \ No newline at end of file diff --git a/posterior_database/models/stan/birthday_demo1.stan b/posterior_database/models/stan/birthday_demo1.stan new file mode 100755 index 00000000..7ae08152 --- /dev/null +++ b/posterior_database/models/stan/birthday_demo1.stan @@ -0,0 +1,123 @@ +functions { + vector gp_pred_exp_quad_rng(real[] x2, + vector y, real[] x1, + real magnitude_1, real length_scale_1, + real magnitude_2, real length_scale_2_1, + real length_scale_2_2, real period, + real sigma) { + int N1 = rows(y); + int N2 = size(x2); + vector[N2] f; + { + matrix[N1, N1] K = add_diag(gp_exp_quad_cov(x1, magnitude_1, length_scale_1) + + gp_periodic_cov(x1, magnitude_2, length_scale_2_1, period) .* + gp_exp_quad_cov(x1, 1.0, length_scale_2_2), 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, magnitude_1, length_scale_1); + f = (k_x1_x2' * K_div_y); + } + return f; + } + vector gp_pred_periodic_exp_quad_rng(real[] x2, + vector y, real[] x1, + real magnitude_1, real length_scale_1, + real magnitude_2, real length_scale_2_1, + real length_scale_2_2, real period, + real sigma) { + int N1 = rows(y); + int N2 = size(x2); + vector[N2] f; + { + matrix[N1, N1] K = add_diag(gp_exp_quad_cov(x1, magnitude_1, length_scale_1) + + gp_periodic_cov(x1, magnitude_2, length_scale_2_1, period) .* + gp_exp_quad_cov(x1, 1.0, length_scale_2_2), 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_periodic_cov(x1, x2, magnitude_2, length_scale_2_1, period) .* + gp_exp_quad_cov(x1, x2, 1.0, length_scale_2_2); + f = (k_x1_x2' * K_div_y); + } + return f; + } + vector gp_pred_full_rng(real[] x2, + vector y, real[] x1, + real magnitude_1, real length_scale_1, + real magnitude_2, real length_scale_2_1, + real length_scale_2_2, real period, + real sigma) { + int N1 = rows(y); + int N2 = size(x2); + vector[N2] f; + { + matrix[N1, N1] K = add_diag(gp_exp_quad_cov(x1, magnitude_1, length_scale_1) + + gp_periodic_cov(x1, magnitude_2, length_scale_2_1, period) .* + gp_exp_quad_cov(x1, 1.0, length_scale_2_2), 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_periodic_cov(x1, x2, magnitude_2, length_scale_2_1, period) .* + gp_exp_quad_cov(x1, x2, 1.0, length_scale_2_2) + + gp_exp_quad_cov(x1, 1.0, length_scale_2_2); + f = (k_x1_x2' * K_div_y); + } + return f; + } +} +data { + int N; + real xn[N]; + vector[N] yn; +} +transformed data { + vector[N] mu = rep_vector(0, N); +} +parameters { + // kernel 1, smooth non-periodic component + real length_scale_1; + real magnitude_1; + + // kernel 2, periodic component + real length_scale_2_1; + real magnitude_2; + real length_scale_2_2; + real period; + + real sigma; +} +transformed parameters { + matrix[N, N] L_K; + { + matrix[N, N] K = gp_exp_quad_cov(xn, magnitude_1, length_scale_1) + + gp_periodic_cov(xn, magnitude_2, length_scale_2_1, period) .* + gp_exp_quad_cov(xn, 1.0, length_scale_2_2); + K = add_diag(K, sigma); + L_K = cholesky_decompose(K); + } +} +model { + // kernel 1 priors + length_scale_1 ~ lognormal(log(10), log(2)); + magnitude_1 ~ lognormal(log(1), log(2)); + + // kernel 2 priors + length_scale_2_1 ~ lognormal(log(2), log(2)); + magnitude_2 ~ lognormal(log(.05), log(2)); + length_scale_2_2 ~ lognormal(log(20), log(2)); + period ~ normal(7, .01); + + sigma ~ normal(0, 1); + + yn ~ multi_normal_cholesky(mu, L_K); +} +generated quantities { + vector[N] f1_pred = gp_pred_exp_quad_rng(xn, yn, xn, magnitude_1, length_scale_1, + magnitude_2, length_scale_2_1, + length_scale_2_2, period, sigma); + vector[N] f2_pred = gp_pred_periodic_exp_quad_rng(xn, yn, xn, magnitude_1, length_scale_1, + magnitude_2, length_scale_2_1, + length_scale_2_2, period, sigma); +} + diff --git a/posterior_database/models/stan/birthday_demo1_2.stan b/posterior_database/models/stan/birthday_demo1_2.stan new file mode 100755 index 00000000..f3de9d4e --- /dev/null +++ b/posterior_database/models/stan/birthday_demo1_2.stan @@ -0,0 +1,162 @@ +functions { + vector gp_pred_exp_quad_rng(real[] x2, + vector y, real[] x1, + matrix k, + real magnitude, real length_scale, real sigma, // squared exp params + model noise + real jitter) { + // x2: test data + // x1: training data + // magnitude: magnitude of squared exponential kernel + // length_scale: length_scale of squared exponential kernel + // sigma: model noise + // jitter: for numerical stability for inverse, etc. + + int N1 = rows(y); + int N2 = size(x2); + vector[N2] f; + { + matrix[N1, N1] K = k;// + diag_matrix(rep_vector(square(sigma), N1)); + 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 = cov_exp_quad(x1, x2, magnitude, length_scale); + vector[N2] f_mu = (k_x1_x2' * K_div_y); + + matrix[N1, N2] v_pred = mdivide_left_tri_low(L_K, k_x1_x2); + matrix[N2, N2] cov_f = cov_exp_quad(x2, magnitude, length_scale) - v_pred' * v_pred + + diag_matrix(rep_vector(jitter, N2)); + vector[N2] temp; + matrix[N2, N2] cov_f_diag; + for (n2 in 1:N2) + temp[n2] = cov_f[n2, n2]; + cov_f_diag = diag_matrix(temp); + +// f = multi_normal_rng(f_mu, cov_f_diag); + f = f_mu; + } + return f; + } +/////////////////////////// +// periodic + squared exponential prediction, +// for kernels 3 and 4 in the sum + vector gp_pred_periodic_exp_quad_rng(real[] x2, + vector y, real[] x1, + matrix k, + real magnitude_1, real length_scale_1, + real magnitude_2, real length_scale_2, real period, + real sigma, real jitter) { + // x2: test data + // x1: training data + // magnitude_1: magnitude of squared exponential kernel + // length_scale_1 length scale of squared exponential kernel + // magnitude_2: magnitude of periodic kernel + // length_scale_2: length_scale of periodic kernel + // period: perio of periodic kernel + // sigma: model noise + // jitter: for numerical stability for inverse, etc. + + int N1 = rows(y); + int N2 = size(x2); + vector[N2] f; + { + matrix[N1, N1] K = k; // + diag_matrix(rep_vector(square(sigma), N1)); + 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_periodic_cov(x1, x2, magnitude_2, length_scale_2, period) .* + cov_exp_quad(x1, x2, magnitude_1, length_scale_1); + vector[N2] f_mu = (k_x1_x2' * K_div_y); + + matrix[N1, N2] v_pred = mdivide_left_tri_low(L_K, k_x1_x2); + matrix[N2, N2] cov_f = gp_periodic_cov(x2, magnitude_2, length_scale_2, period) .* + cov_exp_quad(x2, magnitude_1, length_scale_1) - v_pred' * v_pred // K_t - (K_tn) ( K_n + I * sigma^2)^-1 K_nt) + sigma^2 I + + diag_matrix(rep_vector(jitter, N2)); + vector[N2] temp; + matrix[N2, N2] cov_f_diag; + for (n2 in 1:N2) + temp[n2] = cov_f[n2, n2]; + cov_f_diag = diag_matrix(temp); + +// f = multi_normal_rng(f_mu, cov_f_diag); + f = f_mu; + } + return f; + } +} +data { + int N; + real xn[N]; + vector[N] yn; +} +transformed data { + vector[N] mu = rep_vector(0, N); +} +parameters { + // kernel 1, smooth non-periodic component + real length_scale_1; + real magnitude_1; + + // kernel 2, periodic component + real length_scale_2_1; + real magnitude_2; + real length_scale_2_2; + + // student t priors + real stud_t_sigma; + real stud_t_nu; + + // model noise + real sigma; + vector[N] eta; +} +transformed parameters { + vector[N] f; + matrix[N, N] L_K; + matrix[N, N] K = cov_exp_quad(xn, magnitude_1, length_scale_1) + + gp_periodic_cov(xn, magnitude_2, length_scale_2_1, 7.0) .* + cov_exp_quad(xn, 1.0, length_scale_2_2); + + real sq_sigma = square(sigma); + for (n in 1:N) { + K[n, n] = K[n, n] + sq_sigma; + } + L_K = cholesky_decompose(K); + f = L_K * eta; +} +model { + // kernel 1 priors + length_scale_1 ~ lognormal(log(10), log(2)); // may be second parameter is wrong?? + magnitude_1 ~ lognormal(log(1), log(2)); + + // kernel 2 priors + length_scale_2_1 ~ lognormal(log(2), log(2)); + magnitude_2 ~ lognormal(log(.5), log(2)); + length_scale_2_2 ~ lognormal(log(20), log(2)); + + // student_t likelihood priors + // stud_t_sigma ~ lognormal(log(.05), log(2)); + // stud_t_nu ~ lo + + eta ~ normal(0, 1); + + target += student_t_lpdf(yn | stud_t_nu, mu, stud_t_sigma); +} +generated quantities { + // f1 predictive + vector[N] f1_pred; + vector[N] y1_pred; + // f2 predictive + vector[N] f2_pred; + vector[N] y2_pred; + + // f1_predictive + f1_pred = gp_pred_exp_quad_rng(xn, yn, xn, K, magnitude_1, length_scale_1, sigma, 1e-6); + + // f2 predictive + f2_pred = gp_pred_periodic_exp_quad_rng(xn, yn, xn, K, + magnitude_2, length_scale_2_1, + 1.0, length_scale_2_2, 7, + sigma, 1e-6); +} \ No newline at end of file diff --git a/posterior_database/models/stan/birthday_demo1_mpi.stan b/posterior_database/models/stan/birthday_demo1_mpi.stan new file mode 100755 index 00000000..3dc35baf --- /dev/null +++ b/posterior_database/models/stan/birthday_demo1_mpi.stan @@ -0,0 +1,54 @@ +functions { + matrix gp(real[] x, vector params) { + int N = size(x); + matrix[N, N] K = gp_exp_quad_cov(x, params[1], params[2]) + + + gp_periodic_cov(x, params[3], params[4], 7.0) .* + gp_exp_quad_cov(x, 1.0, params[5]); + for (n in 1:N) K[n, n] = K[n, n] + 1e-6; + return K; + } +} +data { + int N; + real xn[N]; + vector[N] yn; +} +transformed data { + vector[N] mu = rep_vector(0, N); +} +parameters { + // kernel 1, smooth non-periodic component + real length_scale_1; + real magnitude_1; + + // kernel 2, periodic component + real length_scale_2_1; + real magnitude_2; + real length_scale_2_2; +} +transformed parameters { + real params[5]; + matrix[N, N] L_K; + params[1] = length_scale_1; + params[2] = magnitude_1; + params[3] = length_scale_2_1; + params[4] = magnitude_2; + params[5] = length_scale_2_2; + { + // L_K = cholesky_decompose(K); + } +} +model { + + // kernel 1 priors + length_scale_1 ~ lognormal(log(10), log(2)); // may be second parameter is wrong?? + magnitude_1 ~ lognormal(log(1), log(2)); + + // kernel 2 priors + length_scale_2_1 ~ lognormal(log(2), log(2)); + magnitude_2 ~ lognormal(log(.5), log(2)); + length_scale_2_2 ~ lognormal(log(20), log(2)); + + // yn ~ multi_normal_cholesky(mu, L_K); +} diff --git a/posterior_database/models/stan/birthday_demo1_no_noise.stan b/posterior_database/models/stan/birthday_demo1_no_noise.stan new file mode 100755 index 00000000..5d8cc73d --- /dev/null +++ b/posterior_database/models/stan/birthday_demo1_no_noise.stan @@ -0,0 +1,97 @@ +functions { + vector gp_pred_exp_quad_rng(real[] x2, + vector y, real[] x1, + real magnitude_1, real length_scale_1, + real magnitude_2, real length_scale_2_1, + real length_scale_2_2) { + // real sigma) { + int N1 = rows(y); + int N2 = size(x2); + vector[N2] f; + { + matrix[N1, N1] K = add_diag(gp_exp_quad_cov(x1, magnitude_1, length_scale_1) + + gp_periodic_cov(x1, magnitude_2, length_scale_2_1, 7.0) .* + gp_exp_quad_cov(x1, 1.0, length_scale_2_2), 1e-6); + 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, magnitude_1, length_scale_1); + f = (k_x1_x2' * K_div_y); + } + return f; + } + vector gp_pred_periodic_exp_quad_rng(real[] x2, + vector y, real[] x1, + real magnitude_1, real length_scale_1, + real magnitude_2, real length_scale_2_1, + real length_scale_2_2) { + // real sigma) { + int N1 = rows(y); + int N2 = size(x2); + vector[N2] f; + { + /* matrix[N1, N1] K = add_diag(gp_exp_quad_cov(x1, magnitude_1, length_scale_1) + */ + /* gp_periodic_cov(x1, magnitude_2, length_scale_2_1, 7.0) .* */ + /* gp_exp_quad_cov(x1, 1.0, length_scale_2_2), sigma); */ + matrix[N1, N1] K = add_diag(gp_exp_quad_cov(x1, magnitude_1, length_scale_1) + + gp_periodic_cov(x1, magnitude_2, length_scale_2_1, 7.0) .* + gp_exp_quad_cov(x1, 1.0, length_scale_2_2), 1e-6); + 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_periodic_cov(x1, x2, magnitude_2, length_scale_2_1, 7.0) .* + gp_exp_quad_cov(x1, x2, 1.0, length_scale_2_2); + f = (k_x1_x2' * K_div_y); + } + return f; + } +} +data { + int N; + real xn[N]; + vector[N] yn; +} +transformed data { + vector[N] mu = rep_vector(0, N); +} +parameters { + // kernel 1, smooth non-periodic component + real length_scale_1; + real magnitude_1; + + // kernel 2, periodic component + real length_scale_2_1; + real magnitude_2; + real length_scale_2_2; +} +transformed parameters { + matrix[N, N] L_K; + { + matrix[N, N] K = gp_exp_quad_cov(xn, magnitude_1, length_scale_1) + + gp_periodic_cov(xn, magnitude_2, length_scale_2_1, 7.0) .* + gp_exp_quad_cov(xn, 1.0, length_scale_2_2); + K = add_diag(K, 1e-6); + L_K = cholesky_decompose(K); + } +} +model { + // kernel 1 priors + length_scale_1 ~ lognormal(log(10), log(2)); + magnitude_1 ~ lognormal(log(1), log(2)); + + // kernel 2 priors + length_scale_2_1 ~ lognormal(log(2), log(2)); + magnitude_2 ~ lognormal(log(.5), log(2)); + length_scale_2_2 ~ lognormal(log(20), log(2)); + + yn ~ multi_normal_cholesky(mu, L_K); +} +generated quantities { + vector[N] f1_pred = gp_pred_exp_quad_rng(xn, yn, xn, magnitude_1, length_scale_1, + magnitude_2, length_scale_2_1, + length_scale_2_2); + vector[N] f2_pred = gp_pred_periodic_exp_quad_rng(xn, yn, xn, magnitude_1, length_scale_1, + magnitude_2, length_scale_2_1, + length_scale_2_2); +} + diff --git a/posterior_database/models/stan/birthday_demo2.stan b/posterior_database/models/stan/birthday_demo2.stan new file mode 100755 index 00000000..db14cf0f --- /dev/null +++ b/posterior_database/models/stan/birthday_demo2.stan @@ -0,0 +1,369 @@ +functions { + vector gp_pred_exp_quad_rng(real[] x2, + vector y, real[] x1, + vector[] I_s, vector[] I_ss, vector[] I_ws, + real mag, real len, + + real magnitude_1, real length_scale_1, + real magnitude_2, real length_scale_2, + + real magnitude_3, real length_scale_3_1, + real length_scale_3_2, + + real magnitude_4, real length_scale_4_1, + real length_scale_4_2, + + real magnitude_5_1, real magnitude_5_2, real magnitude_5_3, + real jitter) { + int N1 = rows(y); + int N2 = size(x2); + vector[N2] f; + { + matrix[N1, N1] K = cov_exp_quad(x1, magnitude_1, length_scale_1) + + cov_exp_quad(x1, magnitude_2, length_scale_2) + + gp_periodic_cov(x1, magnitude_3, length_scale_3_1, 7) .* + cov_exp_quad(x1, 1.0, length_scale_3_2) + + gp_periodic_cov(x1, magnitude_4, length_scale_4_1, 365.25) .* + cov_exp_quad(x1, 1.0, length_scale_4_2) + + gp_dot_prod_cov(I_s, magnitude_5_1) + + gp_dot_prod_cov(I_ss, magnitude_5_2) + + gp_dot_prod_cov(I_ws, magnitude_5_3) + diag_matrix(rep_vector(jitter, N2)); + 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 = cov_exp_quad(x1, x2, mag, len); + vector[N2] f_mu = (k_x1_x2' * K_div_y); + matrix[N1, N2] v_pred = mdivide_left_tri_low(L_K, k_x1_x2); + matrix[N2, N2] cov_f = cov_exp_quad(x2, mag, len) - v_pred' * v_pred + + diag_matrix(rep_vector(jitter, N2)); + /* vector[N2] temp; */ + /* matrix[N2, N2] cov_f_diag; */ + /* for (n2 in 1:N2) */ + /* temp[n2] = cov_f[n2, n2]; */ + /* cov_f_diag = diag_matrix(temp); */ + /* f = multi_normal_rng(f_mu, cov_f_diag); */ + f = f_mu; + } + return f; + } +/////////////////////////// +// periodic + squared exponential prediction, +// for kernels 3 and 4 in the sum + vector gp_pred_periodic_exp_quad_rng(real[] x2, + vector y, real[] x1, + vector[] I_s, vector[] I_ss, vector[] I_ws, + + real mag_2, real len_2, + real mag_3, real len_3_1, real len_3_2, + real period, + + real magnitude_1, real length_scale_1, + real magnitude_2, real length_scale_2, + + real magnitude_3, real length_scale_3_1, + real length_scale_3_2, + + real magnitude_4, real length_scale_4_1, + real length_scale_4_2, + + real magnitude_5_1, real magnitude_5_2, real magnitude_5_3, + + real jitter) { + int N1 = rows(y); + int N2 = size(x2); + vector[N2] f; + { + matrix[N1, N1] K = cov_exp_quad(x1, magnitude_1, length_scale_1) + + cov_exp_quad(x1, magnitude_2, length_scale_2) + + gp_periodic_cov(x1, magnitude_3, length_scale_3_1, 7) .* + cov_exp_quad(x1, 1.0, length_scale_3_2) + + gp_periodic_cov(x1, magnitude_4, length_scale_4_1, 365.25) .* + cov_exp_quad(x1, 1.0, length_scale_4_2) + + gp_dot_prod_cov(I_s, magnitude_5_1) + + gp_dot_prod_cov(I_ss, magnitude_5_2) + + gp_dot_prod_cov(I_ws, magnitude_5_3) + + diag_matrix(rep_vector(jitter, N2)); + 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_periodic_cov(x1, x2, mag_3, len_3_1, period) .* + cov_exp_quad(x1, x2, mag_2, len_2); + vector[N2] f_mu = (k_x1_x2' * K_div_y); + + matrix[N1, N2] v_pred = mdivide_left_tri_low(L_K, k_x1_x2); + matrix[N2, N2] cov_f = gp_periodic_cov(x2, mag_3, len_3_1, period) .* + cov_exp_quad(x2, mag_2, len_2) - v_pred' * v_pred + + diag_matrix(rep_vector(jitter, N2)); + + /* vector[N2] temp; */ + /* matrix[N2, N2] cov_f_diag; */ + /* f_mu = (k_x1_x2' * K_div_y); */ + /* for (n2 in 1:N2) */ + /* temp[n2] = cov_f[n2, n2]; */ + /* cov_f_diag = diag_matrix(temp); */ + + /* f = multi_normal_rng(f_mu, cov_f_diag); */ + f = f_mu; + } + return f; + } +/////////////////////////// + vector gp_pred_dot_prod_rng(real[] x2, + vector y, real[] x1, + vector[] I_s, vector[] I_ss, vector[] I_ws, + + real sig, + + real magnitude_1, real length_scale_1, + real magnitude_2, real length_scale_2, + + real magnitude_3, real length_scale_3_1, + real length_scale_3_2, + + real magnitude_4, real length_scale_4_1, + real length_scale_4_2, + + real magnitude_5_1, real magnitude_5_2, real magnitude_5_3, + real jitter) { + int N1 = rows(y); + int N2 = size(x2); + vector[N2] f; + { + matrix[N1, N1] K = cov_exp_quad(x1, magnitude_1, length_scale_1) + + cov_exp_quad(x1, magnitude_2, length_scale_2) + + gp_periodic_cov(x1, magnitude_3, length_scale_3_1, 7) .* + cov_exp_quad(x1, 1.0, length_scale_3_2) + + gp_periodic_cov(x1, magnitude_4, length_scale_4_1, 365.25) .* + cov_exp_quad(x1, 1.0, length_scale_4_2) + + gp_dot_prod_cov(I_s, magnitude_5_1) + + gp_dot_prod_cov(I_ss, magnitude_5_2) + + gp_dot_prod_cov(I_ws, magnitude_5_3) + diag_matrix(rep_vector(jitter, N2)); + 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); + vector[N2] f_mu = (k_x1_x2' * K_div_y); + + matrix[N1, N2] v_pred = mdivide_left_tri_low(L_K, k_x1_x2); + matrix[N2, N2] cov_f = gp_dot_prod_cov(x2, sig) - v_pred' * v_pred + + diag_matrix(rep_vector(jitter, N2)); + /* vector[N2] temp; */ + /* matrix[N2, N2] cov_f_diag; */ + /* for (n2 in 1:N2) */ + /* temp[n2] = cov_f[n2, n2]; */ + /* cov_f_diag = diag_matrix(temp); */ + + // f = multi_normal_rng(f_mu, cov_f_diag); + f = f_mu; + } + return f; + } +} +/////////////////////////// +data { + int N; + int N_pred; + int M; + real x[N]; + vector[N] y; + + // for the weekends and special days + vector[M] I_s[N]; + vector[M] I_ss[N]; + vector[M] I_ws[N]; + real x_pred[N_pred]; +} +transformed data { + vector[N] mu = rep_vector(0, N); +} +parameters { + // f_1(t), squared exponential + real length_scale_1; + real magnitude_1; + + // f_2(t), squared exponential + real length_scale_2; + real magnitude_2; + + // f_3(t), weekly quasi-periodic + real length_scale_3_1; + real magnitude_3_1; + real length_scale_3_2; + + // f_4(t), yearly smooth seasonal + real length_scale_4_1; + real magnitude_4_1; + real length_scale_4_2; + + // f_5(t), indicators for weekends and special days + real magnitude_5_1; + real magnitude_5_2; + real magnitude_5_3; +} +transformed parameters { + matrix[N, N] L_k; + { + matrix[N, N] k = cov_exp_quad(x, magnitude_1, length_scale_1) + + cov_exp_quad(x, magnitude_2, length_scale_2) + + gp_periodic_cov(x, magnitude_3_1, length_scale_3_1, 7) .* + cov_exp_quad(x, 1.0, length_scale_3_2) + + gp_periodic_cov(x, magnitude_4_1, length_scale_4_1, 365.25) .* + cov_exp_quad(x, 1.0, length_scale_4_2) + + gp_dot_prod_cov(I_s, magnitude_5_1) + + gp_dot_prod_cov(I_ss, magnitude_5_2) + + gp_dot_prod_cov(I_ws, magnitude_5_3); + for (n in 1:N) + k[n, n] = k[n, n] + 1e-12; + L_k = cholesky_decompose(k); + } +} +model { +// matrix[N, N] L_k; +// L_k = cholesky_decompose(k); + + // smooth non-periodic component + length_scale_1 ~ lognormal(log(365), 1); + magnitude_1 ~ normal(0, 1); + + // faster changing non-periodic component + length_scale_2 ~ lognormal(log(10), 1); + magnitude_2 ~ normal(0, 1); + + // 7 day period + length_scale_3_1 ~ lognormal(log(2), sqrt(2)); + magnitude_3_1 ~ normal(0, 1); + length_scale_3_2 ~ lognormal(log(20), 1); + + // 365.25 day period + length_scale_4_1 ~ lognormal(log(100), sqrt(2)); + magnitude_4_1 ~ normal(0, 1); + length_scale_4_2 ~ lognormal(log(1000), sqrt(2)); + + magnitude_5_1 ~ normal(0, 1); + magnitude_5_2 ~ normal(0, 1); + magnitude_5_3 ~ normal(0, 1); + + y ~ multi_normal_cholesky(mu, L_k); +} +/////////////////////////// +generated quantities { + + // f1 predictive + vector[N_pred] f1_pred; + vector[N_pred] y1_pred; + // f2 predictive + vector[N_pred] f2_pred; + vector[N_pred] y2_pred; + // f3 predictive + vector[N_pred] f3_pred; + vector[N_pred] y3_pred; + // f4 predictive + vector[N_pred] f4_pred; + vector[N_pred] y4_pred; + // f5_1 predictive + vector[N_pred] f5_pred; + vector[N_pred] y5_pred; + + // f1 predictive + f1_pred = gp_pred_exp_quad_rng(x_pred, y, x, + I_s, I_ss, I_ws, + + magnitude_1, length_scale_1, + + 1.0, length_scale_1, + magnitude_2, length_scale_2, + + magnitude_3_1, length_scale_3_1, + length_scale_3_2, + + magnitude_4_1, length_scale_4_1, + length_scale_4_2, + + magnitude_5_1, magnitude_5_2, magnitude_5_3, + + 1e-6); + + // f2 predictive + f2_pred = gp_pred_exp_quad_rng(x_pred, y, x, + I_s, I_ss, I_ws, + + magnitude_2, length_scale_2, + + magnitude_1, length_scale_1, + magnitude_2, length_scale_2, + + magnitude_3_1, length_scale_3_1, + length_scale_3_2, + + magnitude_4_1, length_scale_4_1, + length_scale_4_2, + + magnitude_5_1, magnitude_5_2, magnitude_5_3, + + 1e-6); + + // f3 predictive + f3_pred = gp_pred_periodic_exp_quad_rng(x_pred, y, x, + I_s, I_ss, I_ws, + + magnitude_2, length_scale_2, + 1.0, length_scale_3_1, length_scale_3_2, + 7.0, + + magnitude_1, length_scale_1, + magnitude_2, length_scale_2, + + magnitude_3_1, length_scale_3_1, + length_scale_3_2, + + magnitude_4_1, length_scale_4_1, + length_scale_4_2, + + magnitude_5_1, magnitude_5_2, magnitude_5_3, + 1e-6); + // f4 predictive + f4_pred = gp_pred_periodic_exp_quad_rng(x_pred, y, x, + I_s, I_ss, I_ws, + + magnitude_2, length_scale_2, + 1.0, length_scale_3_1, length_scale_3_2, + 365.25, + + magnitude_1, length_scale_1, + magnitude_2, length_scale_2, + + magnitude_3_1, length_scale_3_1, + length_scale_3_2, + + magnitude_4_1, length_scale_4_1, + length_scale_4_2, + + magnitude_5_1, magnitude_5_2, magnitude_5_3, + 1e-6); + // f5 predictive + f5_pred = gp_pred_dot_prod_rng(x_pred, y, x, + I_s, I_ss, I_ws, + + magnitude_5_1 + magnitude_5_2 + magnitude_5_3, + + magnitude_1, length_scale_1, + magnitude_2, length_scale_2, + + magnitude_3_1, length_scale_3_1, + length_scale_3_2, + + magnitude_4_1, length_scale_4_1, + length_scale_4_2, + + magnitude_5_1, magnitude_5_2, magnitude_5_3, + 1e-6); + + for (n in 1:N_pred) + { + y1_pred[n] = normal_rng(f1_pred[n], 1.0); + y2_pred[n] = normal_rng(f2_pred[n], 1.0); + y3_pred[n] = normal_rng(f3_pred[n], 1.0); + y4_pred[n] = normal_rng(f3_pred[n], 1.0); + y5_pred[n] = normal_rng(f5_pred[n], 1.0); + } +} diff --git a/posterior_database/models/stan/concrete.stan b/posterior_database/models/stan/concrete.stan new file mode 100755 index 00000000..dcc9dd1c --- /dev/null +++ b/posterior_database/models/stan/concrete.stan @@ -0,0 +1,26 @@ +data { + int N; + int D; + + vector[D] x[N]; + vector[N] y; +} +transformed data { + vector[N] mu; + mu = rep_vector(0, N); +} +parameters { + real magnitude; + real length_scale[D]; +} +model { + matrix[N, N] L_K; + { + matrix[N, N] K; + K = cov_exp_quad(x, magnitude, length_scale); + for (n in 1:N) + K[n, n] = K[n, n] + 1e-6; + L_K = cholesky_decompose(K); + } + y ~ multi_normal_cholesky(mu, L_K); +} \ No newline at end of file diff --git a/posterior_database/models/stan/gp_regression.stan b/posterior_database/models/stan/gp_regression.stan new file mode 100755 index 00000000..69176887 --- /dev/null +++ b/posterior_database/models/stan/gp_regression.stan @@ -0,0 +1,84 @@ +functions { + vector gp_pred_rng(real[] x_pred, + vector y1, vector[] x, + real magnitude, real length_scale, + real sigma) { + int N = rows(y1); + int N_pred = size(x_pred); + vector[N_pred] f2; + { + /* matrix[N, N] K = add_diag(gp_exp_quad_cov(x, magnitude, length_scale), */ + /* sigma); */ + /* matrix[N, N] K = add_diag(gp_matern52_cov(x, magnitude, length_scale), */ + /* sigma); */ + matrix[N, N] K = add_diag(gp_matern32_cov(x, magnitude, length_scale), + sigma); + /* matrix[N, N] K = add_diag(gp_exponential_cov(x, magnitude, length_scale), */ + /* sigma); */ + matrix[N, N] L_K = cholesky_decompose(K); + vector[N] L_K_div_y1 = mdivide_left_tri_low(L_K, y1); + vector[N] K_div_y1 = mdivide_right_tri_low(L_K_div_y1', L_K)'; + /* matrix[N, N_pred] k_x_x_pred = gp_exp_quad_cov(x, x_pred, magnitude, length_scale); */ + /* matrix[N, N_pred] k_x_x_pred = gp_matern52_cov(x, x_pred, magnitude, length_scale); */ + matrix[N, N_pred] k_x_x_pred = gp_matern32_cov(x, x_pred, magnitude, length_scale); + /* matrix[N, N_pred] k_x_x_pred = gp_exponential_cov(x, x_pred, magnitude, length_scale); */ + + f2 = (k_x_x_pred' * K_div_y1); + } + return f2; + } +} +data { + int Nobs; + int Ncens; + int D; + int N; + vector[D] x[N]; // N-D GP + // real x[N]; + vector[N] y; + + int N_pred; + // vector[D] x_pred[N_pred]; + real x_pred[N_pred]; +} +transformed data { + vector[N] mu; + mu = rep_vector(0, N); +} +parameters { + real magnitude; + real length_scale[D]; + real sig; + + real sigma; +} +transformed parameters { + matrix[N, N] L_K; + { + matrix[N, N] K = gp_exp_quad_cov(x, magnitude, length_scale); + /* matrix[N, N] K = gp_matern52_cov(x, magnitude, length_scale); */ + /* matrix[N, N] K = gp_matern32_cov(x, magnitude, length_scale); */ + /* matrix[N, N] K = gp_exponential_cov(x, magnitude, length_scale); */ + /* matrix[N, N] K = gp_periodic_cov(x, 1.0, 1.0, 1.0); */ + /* matrix[N, N] K2 = gp_periodic_cov(x, x, 1.0, 1.0, 1.0); */ + /* matrix[N, N] K = gp_dot_prod_cov(x, 1.0); */ + /* matrix[N, N] K2 = gp_dot_prod_cov(x, x, 1.0); */ + + K = add_diag(K, square(sigma)); + L_K = cholesky_decompose(K); + } +} +model { + magnitude ~ normal(0, 3); + length_scale ~ normal(1, 3); + sig ~ normal(0, 1); + + sigma ~ normal(0, 1); + + y ~ multi_normal_cholesky(mu, L_K); +} +generated quantities { + vector[N_pred] f_pred = gp_pred_rng(x_pred, y, x, magnitude, length_scale, sigma); + vector[N_pred] y_pred; + for (n in 1:N_pred) y_pred[n] = normal_rng(f_pred[n], sigma); +} diff --git a/posterior_database/models/stan/gp_regression_2.stan b/posterior_database/models/stan/gp_regression_2.stan new file mode 100755 index 00000000..b1a28002 --- /dev/null +++ b/posterior_database/models/stan/gp_regression_2.stan @@ -0,0 +1,70 @@ +functions { + vector gp_pred_rng(vector[] x_pred, + vector y1, vector[] x, + real magnitude, real[] length_scale, + real sig_0, + real sigma) { + int N = rows(y1); + int N_pred = size(x_pred); + vector[N_pred] f2; + { + matrix[N, N] K = add_diag(gp_dot_prod_cov(x, sig_0) + + gp_exp_quad_cov(x, magnitude, length_scale), + sigma); + matrix[N, N] L_K = cholesky_decompose(K); + vector[N] L_K_div_y1 = mdivide_left_tri_low(L_K, y1); + vector[N] K_div_y1 = mdivide_right_tri_low(L_K_div_y1', L_K)'; + matrix[N, N_pred] k_x_x_pred = gp_dot_prod_cov(x, x_pred, sig_0) + + gp_exp_quad_cov(x, x_pred, magnitude, length_scale); + f2 = (k_x_x_pred' * K_div_y1); + } + return f2; + } +} +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; + mu = rep_vector(0, N); +} +parameters { + real magnitude; + real length_scale[D]; + + real sig_0; + + real sigma; +} +transformed parameters { + matrix[N, N] L_K; + { + matrix[N, N] K = gp_dot_prod_cov(x, sig_0) + gp_exp_quad_cov(x, magnitude, length_scale); + K = add_diag(K, square(sigma)); + L_K = cholesky_decompose(K); + } +} +model { + magnitude ~ normal(0, 3); + length_scale ~ inv_gamma(5, 5); + sig_0 ~ normal(0, 1); + + sigma ~ normal(0, 1); + + y ~ multi_normal_cholesky(mu, L_K); +} +generated quantities { + vector[N_pred] f_pred = gp_pred_rng(x_pred, y, x, magnitude, length_scale, sig_0,sigma); + vector[N] f_pred_in = gp_pred_rng(x, y, x, magnitude, length_scale, sig_0, sigma); + vector[N_pred] y_pred; + vector[N] y_pred_in; + + for (n in 1:N_pred) y_pred[n] = normal_rng(f_pred[n], sigma); // out of sample predictions + for (n in 1:N) y_pred_in[n] = normal_rng(f_pred_in[n], sigma); // out of sample predictions +} diff --git a/posterior_database/models/stan/gp_regression_ard.stan b/posterior_database/models/stan/gp_regression_ard.stan new file mode 100755 index 00000000..f9df7027 --- /dev/null +++ b/posterior_database/models/stan/gp_regression_ard.stan @@ -0,0 +1,66 @@ +functions { + vector gp_pred_rng(vector[] x_pred, + vector y1, vector[] x, + real magnitude, real[] length_scale, + real sigma) { + int N = rows(y1); + int N_pred = size(x_pred); + vector[N_pred] f2; + { + matrix[N, N] K = add_diag(gp_exp_quad_cov(x, magnitude, length_scale), + sigma); + matrix[N, N] L_K = cholesky_decompose(K); + vector[N] L_K_div_y1 = mdivide_left_tri_low(L_K, y1); + vector[N] K_div_y1 = mdivide_right_tri_low(L_K_div_y1', L_K)'; + matrix[N, N_pred] k_x_x_pred = gp_exp_quad_cov(x, x_pred, magnitude, length_scale); + f2 = (k_x_x_pred' * K_div_y1); + } + return f2; + } +} +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; + mu = rep_vector(0, N); +} +parameters { + real magnitude; + real length_scale[D]; + real sig; + + real sigma; +} +transformed parameters { + matrix[N, N] L_K; + { + matrix[N, N] K = gp_exp_quad_cov(x, magnitude, length_scale); + K = add_diag(K, square(sigma)); + L_K = cholesky_decompose(K); + } +} +model { + magnitude ~ normal(0, 3); + length_scale ~ inv_gamma(5, 5); + sig ~ normal(0, 1); + + sigma ~ normal(0, 1); + + y ~ multi_normal_cholesky(mu, L_K); +} +generated quantities { + vector[N_pred] f_pred = gp_pred_rng(x_pred, y, x, magnitude, length_scale, sigma); + vector[N] f_pred_in = gp_pred_rng(x, y, x, magnitude, length_scale, sigma); + vector[N_pred] y_pred; + vector[N] y_pred_in; + + for (n in 1:N_pred) y_pred[n] = normal_rng(f_pred[n], sigma); + for (n in 1:N) y_pred_in[n] = normal_rng(f_pred_in[n], sigma); +} diff --git a/posterior_database/models/stan/gp_regression_mike.stan b/posterior_database/models/stan/gp_regression_mike.stan new file mode 100755 index 00000000..7afa4e6b --- /dev/null +++ b/posterior_database/models/stan/gp_regression_mike.stan @@ -0,0 +1,61 @@ +functions { + vector gp_pred_rng(vector[] x2, + vector y1, vector[] x1, + real alpha, real rho, real sigma, real delta) { + int N1 = rows(y1); + int N2 = size(x2); + vector[N2] f2; + { + matrix[N1, N1] K = cov_exp_quad(x1, alpha, rho) + + diag_matrix(rep_vector(square(sigma), N1)); + matrix[N1, N1] L_K = cholesky_decompose(K); + + vector[N1] L_K_div_y1 = mdivide_left_tri_low(L_K, y1); + vector[N1] K_div_y1 = mdivide_right_tri_low(L_K_div_y1', L_K)'; + matrix[N1, N2] k_x1_x2 = cov_exp_quad(x1, x2, alpha, rho); + vector[N2] f2_mu = (k_x1_x2' * K_div_y1); + matrix[N1, N2] v_pred = mdivide_left_tri_low(L_K, k_x1_x2); + matrix[N2, N2] cov_f2 = cov_exp_quad(x2, alpha, rho) - v_pred' * v_pred + + diag_matrix(rep_vector(delta, N2)); + f2 = multi_normal_rng(f2_mu, cov_f2); + } + return f2; + } +} + +data { + int N; + int D; + vector[D] x[N]; + vector[N] y; + + int N_pred; + vector[D] x_pred[N_pred]; +} + +parameters { + real rho; + real alpha; + real sigma; +} + +model { + matrix[N, N] cov = cov_exp_quad(x, alpha, rho) + + diag_matrix(rep_vector(square(sigma), N)); + matrix[N, N] L_cov = cholesky_decompose(cov); + + // P[rho < 2.0] = 0.01 + // P[rho > 10] = 0.01 + rho ~ inv_gamma(5, 5); + alpha ~ normal(0, 2); + sigma ~ normal(0, 1); + + y ~ multi_normal_cholesky(rep_vector(0, N), L_cov); +} + +generated quantities { + vector[N_pred] f_pred = gp_pred_rng(x_pred, y, x, alpha, rho, sigma, 1e-10); + vector[N_pred] y_pred; + for (n in 1:N_pred) + y_pred[n] = normal_rng(f_pred[n], sigma); +} diff --git a/posterior_database/models/stan/linear_regression.stan b/posterior_database/models/stan/linear_regression.stan new file mode 100755 index 00000000..0bcab2f4 --- /dev/null +++ b/posterior_database/models/stan/linear_regression.stan @@ -0,0 +1,25 @@ +data { + int N; + int D; + matrix[N, D] x; + vector[N] y; + + /* int N_pred; */ + /* matrix[N_pred, D] x_pred; */ +} +parameters { + real alpha; + vector[D] beta; + real sigma; +} +model { + alpha ~ normal(0, 1); // intercept prior + beta ~ normal(0, 1); // regression coef priors + sigma ~ normal(0, 1); // model noise + + y ~ normal(alpha + x * beta, sigma); // likelihood function +} +generated quantities { + /* vector[N_pred] y_pred; */ + /* for (n in 1:N_pred) y_pred[n] = normal_rng(alpha + x_pred[n] * beta, sigma); */ +} diff --git a/posterior_database/models/stan/loggaussian_survival_gp.stan b/posterior_database/models/stan/loggaussian_survival_gp.stan new file mode 100755 index 00000000..edfd8ad6 --- /dev/null +++ b/posterior_database/models/stan/loggaussian_survival_gp.stan @@ -0,0 +1,127 @@ +functions { + vector gp_pred_exp_quad_rng(vector[] x_pred, + vector y1, vector[] x, + real magnitude, real[] length_scale, real sigma) { + int N1 = rows(y1); + int N2 = size(x_pred); + vector[N2] f2; + { + matrix[N1, N1] K = gp_exp_quad_cov(x, magnitude, length_scale) + + diag_matrix(rep_vector(square(sigma), N1)); + matrix[N1, N1] L_K = cholesky_decompose(K); + + vector[N1] L_K_div_y1 = mdivide_left_tri_low(L_K, y1); + vector[N1] K_div_y1 = mdivide_right_tri_low(L_K_div_y1', L_K)'; + matrix[N1, N2] k_x_x_pred = gp_exp_quad_cov(x, x_pred, magnitude, length_scale); + vector[N2] f2_mu = (k_x_x_pred' * K_div_y1); + matrix[N1, N2] v_pred = mdivide_left_tri_low(L_K, k_x_x_pred); + matrix[N2, N2] cov_f2 = gp_exp_quad_cov(x_pred, magnitude, length_scale) - v_pred' * v_pred + + diag_matrix(rep_vector(1e-6, N2)); + f2 = multi_normal_rng(f2_mu, cov_f2); + } + return f2; + } +} +data { + int N; + int D; + int N_pred; + int n_events; + int n_censored; + vector[D] x[N]; + + + vector[D] x_male_age[N_pred]; + vector[D] x_female_age[N_pred]; + + /* vector[D] x_male_tdi[N_pred]; */ + /* vector[D] x_female_tdi[N_pred]; */ + + // vector[D] x_male_wbc[N_pred]; + // vector[D] x_female_wbc[N_pred]; + + // vector[D] x_female_wbc_tdi1[N_pred]; + // vector[D] x_female_wbc_tdi6[N_pred]; + + vector[N] y; // time for observation n + vector[N] event; + int event_idx[n_events]; + int censored_idx[n_censored]; +} + +transformed data { + vector[N] y_new = y / exp(mean(log(y))); +} +parameters { + real magnitude; + real length_scale[D]; + real sig_0; + real sigma; + vector[N] eta; +} +transformed parameters { + vector[N] f; +{ + matrix[N, N] L_K; + // matrix[N, N] K = gp_dot_prod_cov(x, sig_0) + gp_exp_quad_cov(x, magnitude, length_scale); + matrix[N, N] K = gp_exp_quad_cov(x, magnitude, length_scale); + K = add_diag(K, sigma); + L_K = cholesky_decompose(K); + f = L_K * eta; + } +} +model { + length_scale ~ gamma(2, 2); + magnitude ~ normal(0, 2); + sig_0 ~ normal(0, 2); + sigma ~ normal(0, 1); + eta ~ normal(0, 1); + + target += lognormal_lpdf(y[event_idx] | f[event_idx], sigma); + target += lognormal_lccdf(y[censored_idx] | f[censored_idx], sigma); +} +generated quantities { + vector[N_pred] f_male_age = gp_pred_exp_quad_rng(x_male_age, f, x, + magnitude, length_scale, sigma); + vector[N_pred] f_female_age = gp_pred_exp_quad_rng(x_female_age, f, x, + magnitude, length_scale, sigma); + + /* vector[N_pred] f_male_tdi = gp_pred_rn(x_male_tdi, f, x, */ + /* magnitude, length_scale, sigma); */ + + vector[N_pred] y_male_age; + vector[N_pred] y_female_age; + /* vector[N_pred] y_male_tdi; */ + for (n in 1:N_pred) { + y_male_age[n] = lognormal_rng(f_male_age[n], sigma); + y_female_age[n] = lognormal_rng(f_female_age[n], sigma); + /* y_male_tdi[n] = lognormal_rng(f_male_tdi[n], sigma); */ + } + + // vector[N_pred] f_male_tdi; + // vector[N_pred] y_male_tdi; + // vector[N_pred] f_female_tdi; + // vector[N_pred] y_female_tdi; + + // vector[N_pred] f_male_wbc; + // vector[N_pred] y_male_wbc; + // vector[N_pred] f_female_wbc; + // vector[N_pred] y_female_wbc; + + // vector[N_pred] f_female_wbc_tdi1; + // vector[N_pred] y_female_wbc_tdi1; + // vector[N_pred] f_female_wbc_tdi6; + // vector[N_pred] y_female_wbc_tdi6; + + // f_male_age = gp_pred_exp_quad_rng(x_male_age, f, x, K, 1.0, 1.0, sigma, 1e-12); + // f_female_age = gp_pred_exp_quad_rng(x_female_age, f, x, K, 1.0, 1.0, sigma, 1e-12); + + // f_male_tdi = gp_pred_exp_quad_rng(x_male_tdi, f, x, K, 1.0, 1.0, sigma, 1e-12); + // f_female_tdi = gp_pred_exp_quad_rng(x_female_tdi, y, x, K, 1.0, 1.0, sigma, 1e-12); + + // f_male_wbc = gp_pred_exp_quad_rng(x_male_wbc, f, x, K, 1.0, 1.0, sigma, 1.0, 1e-12); + // f_female_wbc = gp_pred_exp_quad_rng(x_female_wbc, f, x, K, 1.0.0, 1.0, sigma, 1e-12); + + // f_female_wbc_tdi1.0 = gp_pred_exp_quad_rng(x_female_wbc_tdi, 1.0, f, x, K, 1.0, 1.0, sigma, 1.0, 1e-12); + // f_female_wbc_tdi6 = gp_pred_exp_quad_rng(x_female_wbc_tdi6, f, x, K, 1.0, 1.0, sigma, 1e-12); +} diff --git a/posterior_database/models/stan/logistic_gp.stan b/posterior_database/models/stan/logistic_gp.stan new file mode 100755 index 00000000..557c1353 --- /dev/null +++ b/posterior_database/models/stan/logistic_gp.stan @@ -0,0 +1,63 @@ +functions { + vector gp_pred_rng(vector[] x_pred, + vector y1, vector[] x, + real magnitude, real[] length_scale, + real sig) { + int N = rows(y1); + int N_pred = size(x_pred); + vector[N_pred] f2; + { + matrix[N, N] K = gp_matern52_cov(x, magnitude, length_scale); + matrix[N, N] L_K = cholesky_decompose(K); + vector[N] L_K_div_y1 = mdivide_left_tri_low(L_K, y1); + vector[N] K_div_y1 = mdivide_right_tri_low(L_K_div_y1', L_K)'; + matrix[N, N_pred] k_x_x_pred = gp_matern52_cov(x, x_pred, magnitude, length_scale); + f2 = (k_x_x_pred' * K_div_y1); + } + return f2; + } +} +data { + int N; + int D; + vector[D] x[N]; + int y[N]; + + int N_pred; + vector[D] x_pred[N_pred]; +} +parameters { + real magnitude; + real length_scale; + real sig; + vector[N] eta; +} +transformed parameters { + vector[N] f; + { + matrix[N, N] K; + matrix[N, N] L_K; + K = gp_matern52_cov(x, magnitude, length_scale); + K = add_diag(K, 1e-12); + L_K = cholesky_decompose(K); + f = L_K * eta; + } +} +model { + magnitude ~ // what should the mangitude prior be? + length_scale ~ // what should the length scale prior be? + + sig ~ normal(0, 2); + + eta ~ normal(0, 1); + + y ~ bernoulli_logit(f); +} +generated quantities { + vector[N_pred] f_pred = gp_pred_rng(x_pred, f, x, magnitude, length_scale, sig); + int y_pred[N_pred]; + int y_pred_in[N]; + + for (n in 1:N) y_pred_in[n] = bernoulli_logit_rng(f[n]); // in sample prediction + for (n in 1:N_pred) y_pred[n] = bernoulli_logit_rng(f_pred[n]); // out of sample predictions +} diff --git a/posterior_database/models/stan/logistic_gp1.stan b/posterior_database/models/stan/logistic_gp1.stan new file mode 100755 index 00000000..dc452843 --- /dev/null +++ b/posterior_database/models/stan/logistic_gp1.stan @@ -0,0 +1,68 @@ +functions { + vector gp_pred_rng(vector[] x_pred, + vector y1, vector[] x, + real magnitude, real[] length_scale, + real sig_0, + real sigma) { + int N = rows(y1); + int N_pred = size(x_pred); + vector[N_pred] f2; + { + matrix[N, N] K = add_diag(gp_exp_quad_cov(x, magnitude, length_scale), sigma); + matrix[N, N] L_K = cholesky_decompose(K); + vector[N] L_K_div_y1 = mdivide_left_tri_low(L_K, y1); + vector[N] K_div_y1 = mdivide_right_tri_low(L_K_div_y1', L_K)'; + matrix[N, N_pred] k_x_x_pred = gp_exp_quad_cov(x, x_pred, magnitude, length_scale); + f2 = (k_x_x_pred' * K_div_y1); + } + return f2; + } +} +data { + int N; + int D; + vector[D] x[N]; + int y[N]; + + int N_pred; + vector[D] x_pred[N_pred]; +} +parameters { + real magnitude; + real length_scale[D]; + real sig_0; + vector[N] eta; + + real sigma; +} +transformed parameters { + vector[N] f; + { + matrix[N, N] K; + matrix[N, N] L_K; + K = gp_exp_quad_cov(x, magnitude, length_scale); + K = add_diag(K, sigma); + L_K = cholesky_decompose(K); + f = L_K * eta; + } +} +model { + magnitude ~ normal(0, 1); + length_scale ~ inv_gamma(5, 5); + + sig_0 ~ normal(0, 2); + + sigma ~ normal(0, 1); + + eta ~ normal(0, 1); + + y ~ bernoulli_logit(f); +} +generated quantities { + vector[N_pred] f_pred = gp_pred_rng(x_pred, f, x, magnitude, length_scale, sig_0, sigma); + int y_pred[N_pred]; + int y_pred_in[N]; + + for (n in 1:N) y_pred_in[n] = bernoulli_logit_rng(f[n]); + for (n in 1:N_pred) y_pred[n] = bernoulli_logit_rng(f_pred[n]); +} diff --git a/posterior_database/models/stan/simulate_exp_quad.stan b/posterior_database/models/stan/simulate_exp_quad.stan new file mode 100755 index 00000000..3caa7071 --- /dev/null +++ b/posterior_database/models/stan/simulate_exp_quad.stan @@ -0,0 +1,21 @@ +data { + int N; + real x[N]; + + real length_scale; + real magnitude; + real sigma; +} + +transformed data { + matrix[N, N] cov = add_diag(gp_exp_quad_cov(x, magnitude, length_scale), 1e-10); + matrix[N, N] L_cov = cholesky_decompose(cov); +} +parameters {} +model {} +generated quantities { + vector[N] f = multi_normal_cholesky_rng(rep_vector(0, N), L_cov); + vector[N] y; + for (n in 1:N) + y[n] = normal_rng(f[n], sigma); +} diff --git a/posterior_database/models/stan/simulate_exponential.stan b/posterior_database/models/stan/simulate_exponential.stan new file mode 100755 index 00000000..d0acd441 --- /dev/null +++ b/posterior_database/models/stan/simulate_exponential.stan @@ -0,0 +1,21 @@ +data { + int N; + real x[N]; + + real length_scale; + real magnitude; + real sigma; +} + +transformed data { + matrix[N, N] cov = add_diag(gp_exponetial_cov(x, magnitude, length_scale), 1e-10); + matrix[N, N] L_cov = cholesky_decompose(cov); +} +parameters {} +model {} +generated quantities { + vector[N] f = multi_normal_cholesky_rng(rep_vector(0, N), L_cov); + vector[N] y; + for (n in 1:N) + y[n] = normal_rng(f[n], sigma); +} diff --git a/posterior_database/models/stan/simulate_matern52.stan b/posterior_database/models/stan/simulate_matern52.stan new file mode 100755 index 00000000..4a34b6b7 --- /dev/null +++ b/posterior_database/models/stan/simulate_matern52.stan @@ -0,0 +1,21 @@ +data { + int N; + real x[N]; + + real length_scale; + real magnitude; + real sigma; +} + +transformed data { + matrix[N, N] cov = add_diag(gp_matern52_cov(x, magnitude, length_scale), 1e-10); + matrix[N, N] L_cov = cholesky_decompose(cov); +} +parameters {} +model {} +generated quantities { + vector[N] f = multi_normal_cholesky_rng(rep_vector(0, N), L_cov); + vector[N] y; + for (n in 1:N) + y[n] = normal_rng(f[n], sigma); +} diff --git a/posterior_database/models/stan/testing.stan b/posterior_database/models/stan/testing.stan new file mode 100755 index 00000000..8e58b105 --- /dev/null +++ b/posterior_database/models/stan/testing.stan @@ -0,0 +1,23 @@ +data { + int N; + int D; + + vector[D] x[N]; + vector[N] y; +} +transformed data { + vector[N] mu; + mu = rep_vector(0, N); +} +parameters { + real magnitude; + real length_scale[D]; +} +model { + matrix[N, N] L_K; + { + matrix[N, N] K; + K = cov_exp_quad(x[,1:6], magnitude, length_scale); + } + // y ~ multi_normal_cholesky(mu, L_K); +}