// gp_regression.stan

functions {
    real mean_fn(array[] real t, array[] real x, vector ky, real a, real r){
      return dot_product(gp_exp_quad_cov(t, x, a, r)[1],ky);
    }
}
data {
    int<lower=0> n; // number of observations
    int<lower=0> m; // number of grid points
    vector[n] y; // response variables
    array[n] real x; // vector of covariates
    real<lower=0> a;
    real<lower=0> b;
}
transformed data {
    vector[n] mu = rep_vector(0, n);
    array[m] real grid; // vector of grid points
    for (i in 1:m)
    	grid[i] = min(x)+(max(x)-min(x))*(i-1)/(m-1);
}
parameters {
    real<lower=0> alpha; //kernel amplitude
    real<lower=0> rho; //kernel lengthscale
}
transformed parameters {
    matrix[n,n] Sigma = gp_exp_quad_cov(x, alpha, rho);
    for (i in 1:n)
    	Sigma[i,i] += 1;
}
model {
    y ~ multi_student_t(2*a, mu, b/a * Sigma);
}
generated quantities {
    vector[n] Ky = inverse_spd(Sigma) * y;
    vector[m] fn_vals;
    for (i in 1:m)
       fn_vals[i] = mean_fn(segment(grid,i,1),x,Ky,alpha,rho);
}