This is a collection of models that I frequent if I’m doing any modeling with Gaussian processes. Contains logistic regression models with Gaussian process priors, spatial models, survival models, regression models with separate length scales, etc. There’s been some changes to the language, so you need to replace the name of the covariance functions to what’s currently in the language. This was in an attempt to transition the Bayesian modeling examples in BDA3 to Stan from GPstuff, the MATLAB GP library that’s maintained by Aalto university. I’ve used parts of this code for some modeling on medical data and when I worked a contract for the Tampa Bay Rays. No warranty or guaranteed support, but I’m happy to answer any questions at my leisure. Hopefully people can use this code for any modeling. Some of it is commented out for debugging, etc. Some of the survival models are lost, but it’s fairly trivial to replace the likelihood function with a different one.
All of these models use my contributions to the Stan math library, part of the C++ backend to Stan.
Regression Models
data {
int<lower=1> N;
int<lower=1> D;
vector[D] x[N];
vector[N] y;
}
transformed data {
vector[N] mu;
mu = rep_vector(0, N);
}
parameters {
real<lower=0> magnitude;
real<lower=0> 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);
}data {
int<lower=1> N;
int<lower=1> D;
vector[D] x[N];
vector[N] y;
}
transformed data {
vector[N] mu;
mu = rep_vector(0, N);
}
parameters {
real<lower=0> magnitude;
real<lower=0> 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);
}data {
int<lower=1> N;
real x[N];
real<lower=0> length_scale;
real<lower=0> magnitude;
real<lower=0> 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);
}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<lower=1> N;
int<lower=1> D;
vector[D] x[N];
vector[N] y;
int<lower=1> N_pred;
vector[D] x_pred[N_pred];
}
transformed data {
vector[N] mu;
mu = rep_vector(0, N);
}
parameters {
real<lower=0> magnitude;
real<lower=0> length_scale[D];
real<lower=0> sig;
real<lower=0> 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);
}
data { int<lower=1> N; real x[N]; real<lower=0> length_scale; real<lower=0> magnitude; real<lower=0> 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);}data {
int<lower=1> N;
real x[N];
real<lower=0> length_scale;
real<lower=0> magnitude;
real<lower=0> 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);
}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<lower=1> N;
int<lower=1> D;
vector[D] x[N];
vector[N] y;
int<lower=1> N_pred;
vector[D] x_pred[N_pred];
}
transformed data {
vector[N] mu;
mu = rep_vector(0, N);
}
parameters {
real<lower=0> magnitude;
real<lower=0> length_scale[D];
real<lower=0> sig_0;
real<lower=0> 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
}Time Series Models
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<lower=1> N;
int<lower=1> 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<lower=1e-13, upper=400> length_scale_1;
real<lower=1e-13> sigma_1;
// f_2(t), squared exponential
real<lower=1e-13> length_scale_2;
real<lower=1e-13> sigma_2;
// f_3(t), weekly quasi-periodic
real<lower=1e-13> length_scale_3_1;
real<lower=1e-13> sigma_3_1;
real<lower=1e-13> length_scale_3_2;
// f_4(t), yearly smooth seasonal
real<lower=1e-13> length_scale_4_1;
real<lower=1e-13> sigma_4_1;
real<lower=1e-13> length_scale_4_2;
// f_5(t), indicators for weekends and special days
real<lower=1e-13> sigma_5_1;
real<lower=1e-13> sigma_5_2;
// noise
real<lower=1e-13> 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);
}
}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<lower=1> N;
int<lower=1> 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<lower=0> length_scale_1;
real<lower=0> magnitude_1;
// f_2(t), squared exponential
real<lower=0> length_scale_2;
real<lower=0> magnitude_2;
// f_3(t), weekly quasi-periodic
real<lower=0> length_scale_3_1;
real<lower=0> magnitude_3_1;
real<lower=0> length_scale_3_2;
// f_4(t), yearly smooth seasonal
real<lower=0> length_scale_4_1;
real<lower=0> magnitude_4_1;
real<lower=0> length_scale_4_2;
// f_5(t), indicators for weekends and special days
real<lower=0> magnitude_5_1;
real<lower=0> magnitude_5_2;
real<lower=0> 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);
}
}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<lower=1> N; int<lower=1> D; int<lower=1> N_pred; int<lower=1> n_events; int<lower=0> 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<lower=0> magnitude; real<lower=0> length_scale[D]; real<lower=0> sig_0; real<lower=0> 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);}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<lower=1> N;
real xn[N];
vector[N] yn;
}
transformed data {
vector[N] mu = rep_vector(0, N);
}
parameters {
// kernel 1, smooth non-periodic component
real<lower=0> length_scale_1;
real<lower=0> magnitude_1;
// kernel 2, periodic component
real<lower=0> length_scale_2_1;
real<lower=0> magnitude_2;
real<lower=0> length_scale_2_2;
}
transformed parameters {
real<lower=0> 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);
}
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<lower=1> N;
real xn[N];
vector[N] yn;
}
transformed data {
vector[N] mu = rep_vector(0, N);
}
parameters {
// kernel 1, smooth non-periodic component
real<lower=0> length_scale_1;
real<lower=0> magnitude_1;
// kernel 2, periodic component
real<lower=0> length_scale_2_1;
real<lower=0> magnitude_2;
real<lower=0> length_scale_2_2;
}
transformed parameters {
real<lower=0> 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);
}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<lower=1> N;
real xn[N];
vector[N] yn;
}
transformed data {
vector[N] mu = rep_vector(0, N);
}
parameters {
// kernel 1, smooth non-periodic component
real<lower=1e-13> length_scale_1;
real<lower=1e-13> magnitude_1;
// kernel 2, periodic component
real<lower=1e-13> length_scale_2_1;
real<lower=1e-13> magnitude_2;
real<lower=1e-13> length_scale_2_2;
// student t priors
real<lower=0> stud_t_sigma;
real<lower=0> stud_t_nu;
// model noise
real<lower=1e-13> 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);
}Survival Models
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<lower=1> N;
int<lower=1> D;
int<lower=1> N_pred;
int<lower=1> n_events;
int<lower=0> 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<lower=0> magnitude;
real<lower=0> length_scale[D];
real<lower=0> sig_0;
real<lower=0> 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);
}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<lower=1> N;
int<lower=1> D;
int<lower=1> N_pred;
int<lower=1> n_events;
int<lower=0> 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<lower=0> magnitude;
real<lower=0> length_scale[D];
real<lower=0> sig_0;
real<lower=0> 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);
}
likely llc 
Leave a comment