<parameter> is 0, but must be positive! when trying to sample from prior of gaussian process

Modeling
1

Hello Guys,

I’m trying to model a time series using a Gaussian process, akin to do a time series decomposition similar to the birthday model, but I get stuck when trying to sample from the prior. I read the chapter about the birthday model in the Bayesian Workflow book and started out with the code from there, but then decided to not use the Hilbert space basis functions in my first iteration, because they looked a bit intimidating to me and I wanted to get going quickly.

Here is my code:

data {
  int<lower=1> N;
  vector[N] x;
}
transformed data {
  vector[N] mu = rep_vector(0, N);
}
parameters {
  real<lower=0,upper=N> rho;
  real<lower=0> alpha;
  vector[N] y_sim;
}
model {
  matrix[N, N] L_K;
  rho ~ std_normal();
  alpha ~ std_normal();

  matrix[N, N] K = gp_exp_quad_cov(to_array_1d(x), 1.0, 1.0) + diag_matrix(rep_vector(1e-10, N));
  L_K = cholesky_decompose(K);
  y_sim ~ multi_normal_cholesky(rep_vector(0, N), L_K);
}

The values in x are just the integers from 1 to 582. However, despite taking very long, I get

Chain [1] Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain [1] Exception: gp_exp_quad_cov: sigma is 0, but must be positive! (in ‘sim.stan’, line 23, column 2 to column 99)
Chain [1] If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain [1] but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

basically all the time. When I set rho and alpha to concrete 1.0, the informational message disappears and the code runs fine (but still takes very long.)

I must be doing something fundamentally wrong. Any help is appreciated!

2

Sorry this didn’t get answered sooner. For some reason it’s not popping up in our orphaned topics list.

The code doesn’t have 23 lines. Are you sure this is the right code? The code you provide here has a fixed 1.0 for sigma in gp_cov_quad function.

3

Hi Bob,

thanks for your answer. I’m sorry, I was a bit sloppy when I wrote my post. I tried to run my program again yesterday, still getting the same messages. This is the exact code I used:

data {
  int<lower=1> N;      // number of observations
  vector[N] x;         // univariate covariate
}
transformed data {
  vector[N] mu = rep_vector(0, N);
}
parameters {
  real<lower=0,upper=N> rho;
  real<lower=0> alpha;
  real<lower=0> sigma;
  vector[N] y;
}
model {
  matrix[N, N] L_K;
  rho ~ std_normal();
  alpha ~ std_normal();
  sigma ~ std_normal();

  matrix[N, N] K = gp_exp_quad_cov(to_array_1d(x), alpha, rho) + diag_matrix(rep_vector(1e-10, N));
  L_K = cholesky_decompose(K);
  y ~ multi_normal_cholesky(rep_vector(0, N), L_K);
}

Now that I actually use the parameters alpha and rho in the call to gp_exp_quad_cov I get messages telling me that both sigma and length_scale are 0 but must be positive.