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.)