Hello Stan users,
I am trying to re-purpose some of Rob Trangucci’s excellent code on Gaussian processes but I am getting some errors and I am not sure what to do https://github.com/rtrangucci/multi_output_gps
Rejecting initial value:
Error evaluating the log probability at the initial value.
Here is the modified code. I originally commented out some of the program thinking that was the issue but I am still getting the same error. The big change that I am working on is that I need a Poisson model where I can observe covariate level data at the observation level. The end result being that I would have observation level effects that are independent of the GP grid.
functions {
matrix gp_exp_quad_chol(real[] x, real alpha, real len, real jitter) {
int dim_x = size(x);
matrix[dim_x, dim_x] L_K_x;
matrix[dim_x, dim_x] K_x = cov_exp_quad(x, alpha, len);
for (n in 1:dim_x)
K_x[n,n] = K_x[n,n] + jitter;
L_K_x = cholesky_decompose(K_x);
return L_K_x;
}
}
data {
// Model is for outcomes that have been
// observed on a 2 dimensional grid
int<lower=1> N; // Number of observations
int<lower=1> K; // Number of predictors
int<lower=1> dim_x_1; // size of grid dimension 1
int<lower=1> dim_x_2; // size of grid dimension 2
real x_1[dim_x_1]; // locations - 1st dimension
real x_2[dim_x_2]; // locations - 2nd dimension
int<lower=0> y[N]; // Outcome
//matrix[N, K] X; // predictor matrix
int row_i[N]; // row position of each observation
int col_j[N]; // col position of each observation
}
parameters {
//real w0; //intercept
//vector[K] beta; //PH predictors
real<lower=0> len_scale_x_1; // length-scale - 1st dimension
real<lower=0> len_scale_x_2; // length-scale - 2nd dimension
real<lower=0> alpha; // Scale of outcomes
// Standardized latent GP
matrix[dim_x_1, dim_x_2] y_tilde;
}
model {
matrix[dim_x_1, dim_x_2] latent_gp;
{
matrix[dim_x_1, dim_x_1] L_K_x_1 = gp_exp_quad_chol(x_1, 1.0, len_scale_x_1, 1e-12);
matrix[dim_x_2, dim_x_2] L_K_x_2 = gp_exp_quad_chol(x_2, alpha, len_scale_x_2, 1e-12);
// latent_gp is matrix-normal with among-column covariance K_x_1
// among-row covariance K_x_2
latent_gp = L_K_x_1 * y_tilde * L_K_x_2';
}
// priors
len_scale_x_1 ~ gamma(8, 2);
len_scale_x_2 ~ gamma(8, 2);
alpha ~ normal(0, 1);
to_vector(y_tilde) ~ normal(0, 1);
// weakly informative prior for the intercept
//w0 ~ normal(0,10);
//Weakly informative prior for the other variables in data
//beta ~ normal(0,10); //w0 + X * beta +
// likelihood
y ~ poisson_log(to_vector(latent_gp[row_i, col_j]'));
}
Thanks for any help you can provide!
-Patrick