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

Motorcycle

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

Data are modelled with normal distribution having Gaussian process prior on mean and log standard deviation: \[ y \sim \mbox{normal}(\mu(x), \exp(\eta(x))\\ \mu \sim GP(0, K_1)\\ \eta \sim GP(0, K_2) \] \(K_1\) and \(K_2\) are exponentiated quadratic covariance functions.

Load packages

library(cmdstanr) 
library(posterior)
options(pillar.neg = FALSE, pillar.subtle=FALSE, pillar.sigfig=2)
library(tidyr) 
library(dplyr) 
library(ggplot2)
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 reproducability

Motorcycle accident acceleration 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)")