K = gp_exp_quad_cov(X, magnitude, length_scale);
K = add_diag(K, square(eta));
with
K = cov_exp_quad(X, magnitude, length_scale) + diag_matrix(rep_vector(square(eta), N));
Now I get a
Error in cpp_object_initializer(.self, .refClassDef, ...) :
could not find function "cpp_object_initializer"
failed to create the sampler; sampling not done
What version of Stan/math are you using? If you’re on an older version, before 2.18 I believe, you can replace gp_exp_quad_cov with cov_exp_quad, as we’re making the names of the covariance functions more unified. I think the gp_exp_quad_cov is in 2.18.
As to your second question, try the initialization of K in two seperate lines, like so:
K = cov_exp_quad(x, magnitude, length_scale);
K = K + diag_matrix(rep_vector(square(eta), N);
And I think this should solve the compilation issue. With that said, many of the kernels (matern 3/2, matern 5/2, exponential) are not yet apart of stan math, and you can find the implementations here: https://github.com/drezap/math/tree/kernels-beta, if you want to play around with some other kernels (I don’t think ARD works yet, though on this branch, the matern5/2 got merged a few weeks ago).
I’ve been slacking at getting these merged, apologies.
Please let me know if that solution works.
I’d also recommend checking out Betancourt’s case studies, for suggestions as to how to scale the hyper prior distributions, as I’ve neglected this in my cases study.
Hi @anon79882417.
First of all, thank you for your document. It pairs very well with the workflow presented by @betanalpha.
This said, I changed the transformed parameters block following your advice:
transformed parameters {
matrix[N, N] L_K;
{
matrix[N, N] K;
K = cov_exp_quad(X, magnitude, length_scale);
K = K + diag_matrix(rep_vector(square(eta), N));
L_K = cholesky_decompose(K);
}
}
but the error remains:
recompiling to avoid crashing R session
Error in cpp_object_initializer(.self, .refClassDef, ...) :
could not find function "cpp_object_initializer"
failed to create the sampler; sampling not done
Can you please run in command Stan to make sure it’s not an issue with Rstan? I’ve never seen this error before. If I have a syntax error in Rstan, it usually just crashes. command Stan will catch it with a Lint error.
Could you also post the full Stan code? not just the transformed parameters block? This might help.
short long story: I am trying to reproduce the two-dimensional krigging on the mesue dataset, so in this case D=2 and Y corresposds to the zinc varable, that I log-transformed (see transformed data block).
cmdstan (v. 2.18.0) gives me no error: just one warning:
The interpretation of the hyperparameters wasn’t clear to me initially, either. If you want help interpreting them, I mention it (once, I think) in my case study, Betancourt has some info about it, but the main reference for GPs is GPML (Gaussian Processes for Machine Learning) which I think you can find with a google search. Highly recommended you check it out. Please put some thought into the scale of the length-scale hyper parameter, at minimum. There’s also a Trangucci et al paper about it.
On that note - I would look into the posterior predictive distribution. Using this, you can generate predictions for any range of the input data you desire. So if you want to do prediction for spatial data (with lat/long) an approach could be to set x_pred to be an evenly spaced grid along each dimension. I’ve also seen spatial analysis that use census tracts. Generating the posterior predictive is straight forward multivariate gaussian algebra.
I’m finding GPs more useful for prediction or as prior distributions where we need more flexibility, or if I need to use a basis function approach to data analysis. Honestly, it’s a lot easier to interpret linear regression coefficients than the hyper parameters of a GP regression. Also, I think for spatial data people use the Matern5/2 covariance function, not as smooth as the exp quad one. It all depends on the application though. This just got merged into stan/math, and it has ARD “built-in”.
Stan could be less stable for higher dimensional data (i.e. autodiff consumes lots of memory with large covariance matrices, and the posterior induced by a GP prior can be nasty, so HMC sometimes doesn’t behave nicely if the model isn’t well specified or is has a complex geometry), but let us know how it goes!