GP interpolation with brms

I want to use the R package brms to fit a Gaussian process interpolating the points. I want to use Bayesian techniques for hyperparameter selection, inference, and uncertainty propagation. Since I am doing interpolation, I do not want there to be a “residual error” term being fitted, or if there is (maybe for numerical reasons) it needs to be very small.

Below is the closest I got:

library(brms)
testdat <- cbind(x = seq(0, 1, length = 5), y = c(0.1, 1.2, 2, 0.8, -0.4))
testfit <- brm(y ~ gp(x), data = testdat, chains = 2, prior = set_prior("normal(0, 0.000000000001)", class = "sigma"))
plot(conditional_effects(testfit, ndraws = 200, spaghetti = TRUE))

I tried replacing the Normal prior with constant(0) but that resulted in problems, so I opted for a Normal prior with a very small standard deviation. This code kills my R session when I attempt plot (so the problem likely lies in conditional_effects.

(Realize that this is a toy example and that I plan on doing GP interpolation over a 2D plane, and would like to be able to interpolate in arbitrary dimensions in the end.)

So how can I do GP interpolation with brms?

(I also asked this question on CrossValidated.)

How about constant(0.0001) ?

You could also consider creating the Stan code with brms function make_stancode() and then edit that code for your purposes.

Do you plan to have other components than GP is your model? If you need just the GP interpolation, it might be better to use some other software that supports “noiseless” observations in GPs. gplite might be useful.

(A few years ago we used GP interpolation succesfully up to 50 dimensions with “noiseless” electron density functional calculations of the energy and atomic forces, but we used Matlab based GPstuff for the development, and the production version is in C tailored for the specific application)

Thanks for the answer. It seems like it gets to the issue, then there’s a ton of complaints about divergent chains and tree depth and ESS and R-hat and so on. I can’t seem to figure it out, and for my purposes I don’t think it’s worth it. But thank you for the answer.

Ah, brms is using the latent variable presentation of GP, and you definitely want to use the marginalized version for this use case. See how multivariate_normal is used in Stan manual and Gaussian process demonstration with Stan

1 Like