Locally stationary GPs

[seabbs]

Extend the current GP functionality to be locally stationary by using a patch based approach of stationary Gps with the GP parameters being related with hierarchical distributions. The main issue with this approach is what to do at the boundary between patches. In the first instance I think it makes sense do nothing and see how it goes but if/when that fails I think extending patches to overlap and using a sigmoid or similar to transition from one patch to another might be a good way to go.

As a feature this would support better fitting to retrospective data when the lengthscale has changed over time. It might also be a useful forecasting model in some formulations if you tweak the hierarchical priors to have long lengthscales (as this makes the prior model for future patches that they have long lenghtscales (i.e slower change over time with less variance).

[seabbs]

I plan to work on this in the next few months

[seabbs]

@SamuelBrand1 did you have any thoughts as to the best approximate thing we can do at patch boundaries here or anything else we could do that 1. We can implement 2. Might perform better.

[SamuelBrand1]

Hey @seabbs,

My off the top thoughts are:

• I don't see a reason not to at least reason on the whole GP being made up of compositions of stationary kernels via the "change point" kernel (e.g. implemented here). I think in practise this is a sigmoid interpolation between the patches.
• The approximation deployed in https://link.springer.com/article/10.1007/s11222-022-10167-2 (and implemented in EpiNow2) could be thought of as the errors around a patch trend line (this accounts for the Dirichlet zero boundary conditions). From f2f/off-line chat you know I think this has a maths representation as a Gaussian bridge but I haven't ever got round to locking the maths down on that precisely.

In terms of Maths, what I'm thinking (pre apols for any errors, I'm being quick):

• The Change point kernel for two patches is:

$$K(x, x') = K_1(x, x') \cdot \left(1 - \sigma(x)\right)\left(1 - \sigma(x')\right) + K_2(x, x') \cdot \sigma(x)\sigma(x')$$

where $\sigma(x)$ is the logistic function:

$$\sigma(x) = \frac{1}{1 + e^{-s(x - x_0)}}$$

So we have added s (steepness) to the parameters to fit and maybe x_0 (location) (although maybe we choose this).

• The actual $K_1(x,x')$ and $K_2(x,x')$ can be:

$$K_1(x,x') = \beta^2 \langle x, x' \rangle + \tilde{K}_1.$$

Where $\tilde{K}_1$ is the approx kernel from the practical Hilbert space paper and the $\beta^2 \langle x, x' \rangle$ term is from having a linear slope with grad $\beta$. Similarly $K_2$.

• We put the kernel approxs directly into the change point form and then we should have everything?

[SamuelBrand1]

IMO, something like this should be pretty principled. The "true" form is definitely a GP, the computation approx form is well described by the linked paper and has been implemented in EpiNow2 as a "one patch".

[SamuelBrand1]

Btw, other sigmoids work too e.g. tanh