Demonstration of covariance matrix and basis function implementation of Gaussian process model in Stan.

The basics of the covariance matrix approach is based on the Chapter 10 of Stan User’s Guide, Version 2.26 by Stan Development Team (2021). https://mc-stan.org/docs/stan-users-guide/

The basics of the Hilbert space basis function approximation is based on Riutort-Mayol, Bürkner, Andersen, Solin, and Vehtari (2020). Practical Hilbert space approximate Bayesian Gaussian processes for probabilistic programming. https://arxiv.org/abs/2004.11408

1 Motorcycle accident acceleration data

Data are measurements of head acceleration in a simulated motorcycle accident, used to test crash helmets.

Data are modelled first with normal distribution having Gaussian process prior on mean: \[ y \sim \mbox{normal}(f(x), \sigma)\\ f \sim GP(0, K_1)\\ \sigma \sim \mbox{normal}^{+}(0, 1), \] and then with normal distribution having Gaussian process prior on mean and log standard deviation: \[ y \sim \mbox{normal}(f(x), \exp(g(x))\\ f \sim GP(0, K_1)\\ g \sim GP(0, K_2). \]

Load packages

library(cmdstanr) 
library(posterior)
library(tidybayes)
options(pillar.neg = FALSE, pillar.subtle=FALSE, pillar.sigfig=2)
library(tidyr) 
library(dplyr) 
library(ggplot2)
library(ggrepel)
library(bayesplot)
theme_set(bayesplot::theme_default(base_family = "sans", base_size=16))
set1 <- RColorBrewer::brewer.pal(7, "Set1")
SEED <- 48927 # set random seed for reproducibility

1.1 Load and plot data

Load data

data(mcycle, package="MASS")
head(mcycle)
  times accel
1   2.4   0.0
2   2.6  -1.3
3   3.2  -2.7
4   3.6   0.0
5   4.0  -2.7
6   6.2  -2.7

Plot data

mcycle %>%
  ggplot(aes(x=times,y=accel))+
  geom_point()+
  labs(x="Time (ms)", y="Acceleration (g)")