---
title: "Gaussian Process implementation details"
output:
  rmarkdown::html_vignette:
    toc: true
    number_sections: true
bibliography: library.bib
csl: https://raw.githubusercontent.com/citation-style-language/styles/master/apa-numeric-superscript-brackets.csl
vignette: >
  %\VignetteIndexEntry{Gaussian Process implementation details}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

# Overview

We make use of Gaussian Processes in several places in `EpiNow2`. For example, the default model for `estimate_infections()` uses a Gaussian Process to model the 1st order difference on the log scale of the reproduction number. This vignette describes the implementation details of the approximate Gaussian Process used in `EpiNow2`.

# Definition

The single dimension Gaussian Processes ($\mathrm{GP}_t$) we use can be written as

\begin{equation}
\mathcal{GP}(\mu(t), k(t, t'))
\end{equation}

where $\mu(t)$ and $k(t,t')$ are the mean and covariance functions, respectively.
In our case as set out above, we have

\begin{equation}
\mu(t) \equiv 0 \\
k(t,t') = k(|t - t'|) = k(\Delta t)
\end{equation}

with the following choices available for the kernel $k$

## Matérn 3/2 covariance kernel (the default)

\begin{equation}
k(\Delta t) = \alpha^2 \left( 1 + \frac{\sqrt{3} \Delta t}{l} \right) \exp \left( - \frac{\sqrt{3} \Delta t}{l}\right)
\end{equation}

with $l>0$ and $\alpha > 0$ the length scale and magnitude, respectively, of the kernel.
Note that here and later we use a slightly different definition of $\alpha$ compared to Riutort-Mayol et al. [@approxGP], where this is defined as our $\alpha^2$.

## Squared exponential kernel

\begin{equation}
k(\Delta t) = \alpha^2 \exp \left( - \frac{1}{2} \frac{(\Delta t^2)}{l^2} \right)
\end{equation}

## Ornstein-Uhlenbeck (Matérn 1/2) kernel

\begin{equation}
k(\Delta t) = \alpha^2 \exp{\left( - \frac{\Delta t}{2 l^2}  \right)}
\end{equation}

## Matérn 5/2 covariance kernel

\begin{equation}
k(\Delta t) = \alpha \left( 1 + \frac{\sqrt{5} \Delta t}{l} + \frac{5}{3} \left(\frac{\Delta t}{l} \right)^2 \right) \exp \left( - \frac{\sqrt{5} \Delta t}{l}\right)
\end{equation}

# Hilbert space approximation

In order to make our models computationally tractable, we approximate the Gaussian Process using a Hilbert space approximation to the Gaussian Process [@approxGP], centered around mean zero.

\begin{equation}
\mathcal{GP}(0, k(\Delta t)) \approx \sum_{j=1}^m \left(S_k(\sqrt{\lambda_j}) \right)^\frac{1}{2} \phi_j(t) \beta_j
\end{equation}

with $m$ the number of basis functions to use in the approximation, which we calculate from the number of time points $t_\mathrm{GP}$ to which the Gaussian Process is being applied (rounded up to give an integer value), as is recommended [@approxGP].

\begin{equation}
m = b t_\mathrm{GP}
\end{equation}

and values of $\lambda_j$ given by

\begin{equation}
\lambda_j = \left( \frac{j \pi}{2 L} \right)^2
\end{equation}

where $L$ is a positive number termed boundary condition, and $\beta_{j}$ are regression weights with standard normal prior

\begin{equation}
\beta_j \sim \mathcal{Normal}(0, 1)
\end{equation}

The function $S_k(x)$ is the spectral density relating to a particular covariance function $k$.
In the case of the Matérn kernel of order $\nu$ this is given by

\begin{equation}
S_k(x) = \alpha^2 \frac{2 \sqrt{\pi} \Gamma(\nu + 1/2) (2\nu)^\nu}{\Gamma(\nu) \rho^{2 \nu}} \left( \frac{2 \nu}{\rho^2} + \omega^2 \right)^{-\left( \nu + \frac{1}{2}\right )}
\end{equation}

For $\nu = 3 / 2$ (the default in `EpiNow2`) this simplifies to

\begin{equation}
    S_k(\omega) =
    \alpha^2 \frac{4 \left(\sqrt{3} / \rho \right)^3}{\left(\left(\sqrt{3} / \rho\right)^2 + \omega^2\right)^2} = 
    \left(\frac{2 \alpha \left(\sqrt{3} / \rho \right)^{\frac{3}{2}}}{\left(\sqrt{3} / \rho\right)^2 + \omega^2} \right)^2
\end{equation}

For $\nu = 1 / 2$ it is

\begin{equation}
S_k(\omega) = \alpha^2 \frac{2}{\rho \left(1 / \rho^2 + \omega^2\right)}
\end{equation}

and for $\nu = 5 / 2$ it is

\begin{equation}
S_k(\omega) =
    \alpha^2 \frac{3 \left(\sqrt{5} / \rho \right)^5}{2 \left(\left(\sqrt{5} / \rho\right)^2 + \omega^2\right)^3}
\end{equation}

In the case of a squared exponential the spectral kernel density is given by

\begin{equation}
S_k(\omega) = \alpha^2 \sqrt{2\pi} \rho \exp \left( -\frac{1}{2} \rho^2 \omega^2 \right)
\end{equation}

The functions $\phi_{j}(x)$ are the eigenfunctions of the Laplace operator,

\begin{equation}
\phi_j(t) = \frac{1}{\sqrt{L}} \sin\left(\sqrt{\lambda_j} (t^* + L)\right)
\end{equation}

with time rescaled linearly to be between -1 and 1,

\begin{equation}
t^* = \frac{t - \frac{1}{2}t_\mathrm{GP}}{\frac{1}{2}t_\mathrm{GP}}
\end{equation}

Relevant priors are

\begin{align}
\alpha &\sim \mathcal{Normal}(\mu_\alpha, \sigma_{\alpha}) \\
\rho   &\sim \mathcal{LogNormal} (\mu_\rho, \sigma_\rho)\\
\end{align}

with $\rho$ additionally constrained to be between $\rho_\mathrm{min}$ and $\rho_\mathrm{max}$, $\mu_{\rho}$ and $\sigma_\rho$ calculated from given mean $m_{\rho}$ and standard deviation $s_\rho$, and default values (all of which can be changed by the user):

\begin{align}
b &= 0.2 \\
L &= 1.5 \\
m_\rho &= 21 \\
s_\rho &= 7 \\
\rho_\mathrm{min} &= 0\\
\rho_\mathrm{max} &= 60\\
\mu_\alpha &= 0\\
\sigma_\alpha &= 0.01
\end{align}

# References