Using non-approximate GPs in forecasting

[SamuelBrand1]

Problem

The GP options described here give different available forms for Hilbert space approximate GPs (HSGPs) for Rt. This is really neat, and I assume helps sampling. However, given that this relies on sine basis functions I'd want a different option for out-of-sample (i.e. forecasting) to avoid possible extrapolation problems.

Possible solution

Given that EpiNow2 samples:

1. the length scale and variance parameters of the available HS approximate GPs
2. Their latent process trajectory

It is possible to sample their forward trajectories from their "true" GP form and avoid extrapolation problems associated with basis functions.

Once possible way to do this would be to:

• For each posterior sample get sampled HSGP trajectory $X = (X_t)_{t =1, \dots, T}$ at the infection time points, a length scale $\ell$ and variance $\sigma^2$.
  • Generate forward samples of the GP $\tilde{X} = (X_t)_{ t = T+1, T+2, ...}$ by sampling:

$$\tilde{X} \sim f(\tilde{X} | X, \ell, \sigma^2)$$

Which any GP package will do.

• Given forward sample $\tilde{X}$, then continue the Renewal model simulation into future using all relevant posterior sample parameters.

This should let EpiNow2 use posterior samples from the HSGP approximation to sample future trajectories from the "full" GP. This might be an improvement (would need checking) and seems well motivated.

[sbfnk]

This is a really nice idea - the only downside I can see is that it requires maintaining an R simulation engine in parallel to the stan model.

[SamuelBrand1]

Yes, thats the annoying bit from the engineering point of view.

Added context is that I'm doing something similar atm here https://github.com/cdcent/cfa-mech-experiment (which is not open but I know other CDC people sometimes come on here).

What I'm doing is time differential modelling (i.e. ODEs or RODEs) but would work similarly/identically in discrete time albeit in julia

Basic outline:

Forward sampling

function multistrain_seir!(du, u, p, t, W)
# stuff in here with noise process value at time t = W
end

This function signature can be passed to RODE solvers https://docs.sciml.ai/DiffEqDocs/stable/tutorials/rode_example/ where you specify the noise process $W_t$ parameters (e.g. volatility), in my case I'm looking at OU process/RW/brownian bridge. The RODE solver then handles simulating samples.

In-sample inference

I wrap the above inside a function signature like

function multistrain_seir!(du, u, p, t)
W = hilbert_approximation(t,p)
multistrain_seir!(du, u, p, t, W)
end

Where hilbert_approximation brings in the handy Hilbert space approximation which is useful for inference because you can diff through the ODE solver (the sine basis functions are smooth) to do gradient based inference on the parameters.

Stan implementation

I'm not great at stan... But I use a dispatch mechanism to let the code know which approach I want to do. In stan you could have a boolean to select which mode you want? Have the forward sampling version in the generated quantities?