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 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(loo)
library(tidybayes)
options(pillar.neg = FALSE, pillar.subtle=FALSE, pillar.sigfig=2)
library(tidyr)
library(dplyr)
library(ggplot2)
library(ggrepel)
library(latex2exp)
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
```

## 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)")
```