Skip to contents

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 (GPt\mathrm{GP}_t) we use can be written as

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

where μ(t)\mu(t) and k(t,t)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 kk

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

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

with l>0l>0 and α>0\alpha > 0 the length scale and magnitude, respectively, of the kernel.

Squared exponential kernel

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

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

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

Matérn 5/2 covariance kernel

k(Δt)=α(1+5Δtl+53(Δtl)2)exp(5Δtl)\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[1], centered around mean zero.

𝒢𝒫(0,k(Δt))j=1m(Sk(λj))12ϕj(t)βj\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 mm the number of basis functions to use in the approximation, which we calculate from the number of time points tGPt_\mathrm{GP} to which the Gaussian Process is being applied (rounded up to give an integer value), as is recommended[1].

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

and values of λj\lambda_j given by

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

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

βj𝒩𝓇𝓂𝒶𝓁(0,1)\begin{equation} \beta_j \sim \mathcal{Normal}(0, 1) \end{equation}

The function Sk(x)S_k(x) is the spectral density relating to a particular covariance function kk. In the case of the Matérn 3/2 kernel (the default in EpiNow2) this is given by

Sk(x)=4α2(3ρ)3((3ρ)2+w2)2\begin{equation} S_k(x) = 4 \alpha^2 \left( \frac{\sqrt{3}}{\rho}\right)^3 \left(\left( \frac{\sqrt{3}}{\rho} \right)^2 + w^2 \right)^{-2} \end{equation}

and in the case of a squared exponential kernel by

Sk(x)=α22πρexp(12ρ2w2)\begin{equation} S_k(x) = \alpha^2 \sqrt{2\pi} \rho \exp \left( -\frac{1}{2} \rho^2 w^2 \right) \end{equation}

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

ϕj(t)=1Lsin(λj(t*+L))\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,

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

Relevant priors are

α𝒩𝓇𝓂𝒶𝓁(0,σα)ρ𝒩𝓇𝓂𝒶𝓁(μρ,σρ)\begin{align} \alpha &\sim \mathcal{Normal}(0, \sigma_{\alpha}) \\ \rho &\sim \mathcal{LogNormal} (\mu_\rho, \sigma_\rho)\\ \end{align}

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

b=0.2L=1.5mρ=21sρ=7ρmin=0ρmax=60σα=0.05\begin{align} b &= 0.2 \\ L &= 1.5 \\ m_\rho &= 21 \\ s_\rho &= 7 \\ \rho_\mathrm{min} &= 0\\ \rho_\mathrm{max} &= 60\\ \sigma_\alpha &= 0.05\\ \end{align}

References

1. Riutort-Mayol, G., Bürkner, P.-C., Andersen, M. R., Solin, A., & Vehtari, A. (2020). Practical hilbert space approximate bayesian gaussian processes for probabilistic programming. https://arxiv.org/abs/2004.11408