Email

Previously I looked at how to simulate Gaussian processes in R, following the methods in Rasmussen and Williams. But now that Andrew Gelman et al. (of Bayesian Data Analysis, and Data Analysis Using Regression and Multilevel/Hierarchical Models fame) have released their new Bayesian tool Stan, I wanted to redo the earlier analysis using the related Rstan package.

Turns out it’s pretty easy. The project team have provided abundant examples and so it was just a case of tweaking these demos to reproduce the earlier examples. I’ve created a Gist with all of the necessary code and here are the resulting plots, just to prove that it works.

Example of functions from a Gaussian process calculated with RStan

Example of Gaussian process trained on noise-free data, calculated with RStan

Example of Gaussian process trained on noisy data, calculated with RStan

Beware however that the Stan implementations are much slower than the manual method, as shown in the timings below, run for the basic unconstrained process.
## Manual method user system elapsed 0.38 0.00 0.39 ## RStan method user system elapsed 28.59 0.05 59.30

Of course if you already know the parameters of your covariance function, then using Stan is a bit like using a sledgehammer to crack a walnut. But if you have a very complex model that requires fitting, then Stan is ideal as it provides user-friendly syntax and powerful sampling for Bayesian inference.