Processing math: 100%

Gaussian Process implementation details

1 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.

2 Definition

The one-dimensional Gaussian Processes (GPt) we use can be written as

GP(μ(t),k(t,t))

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

μ(t)0k(t,t)=k(|tt|)=k(Δt)

with the following choices available for the kernel k

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

k(Δt)=α2(1+3Δtρ)exp(3Δtρ)

with ρ>0 and α>0 the length scale and magnitude, respectively, of the kernel. Note that here and later we use a slightly different definition of α compared to Riutort-Mayol et al.[1], where this is defined as our α2.

2.2 Squared exponential kernel

k(Δt)=α2exp(12(Δt2)l2)

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

k(Δt)=α2exp(Δt2ρ2)

2.4 Matérn 5/2 covariance kernel

k(Δt)=α(1+5Δtρ+53(Δtl)2)exp(5Δtρ)

3 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.

GP(0,k(Δt))mj=1(Sk(λj))12ϕj(t)βj

with m the number of basis functions to use in the approximation, which we calculate from the number of time points tGP to which the Gaussian Process is being applied (rounded up to give an integer value), as is recommended[1].

m=btGP

and values of λj given by

λj=(jπ2L)2

where L is a positive number termed boundary condition, and βj are regression weights with standard normal prior

βjNormal(0,1)

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

Sk(x)=α22πΓ(ν+1/2)(2ν)νΓ(ν)ρ2ν(2νρ2+ω2)(ν+12)

For ν=3/2 (the default in EpiNow2) this simplifies to

Sk(ω)=α24(3/ρ)3((3/ρ)2+ω2)2=(2α(3/ρ)32(3/ρ)2+ω2)2

For ν=1/2 it is

Sk(ω)=α22ρ(1/ρ2+ω2)

and for ν=5/2 it is

Sk(ω)=α23(5/ρ)52((5/ρ)2+ω2)3

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

Sk(ω)=α22πρexp(12ρ2ω2)

The functions ϕj(x) are the eigenfunctions of the Laplace operator,

ϕj(t)=1Lsin(λj(t+L))

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

t=t12tGP12tGP

Relevant default priors are

αHalfNormal(0,0.01)ρLogNormal(μρ,σρ)

with ρ additionally constrained with an upper bound of 60 and μρ and σρ calculated using a mean of 21 and standard deviation of 7.

Furthermore, by default we set.

b=0.2L=1.5

These values as well as the prior distributions of relevant parameters can all be changed by the user.

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