How meaningful are Gaussian ODE filter uncertainties?¶

We investigate whether the uncertaintes returned by Gaussian ODE Filters are meaningful as actual error estimates—or even better?

# Make inline plots vector graphics instead of raster graphics
%matplotlib inline
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('pdf', 'svg')

# Plotting
import matplotlib.pyplot as plt'../../probnum.mplstyle')
import numpy as np

import probnum.diffeq as diffeq
from probnum.diffeq import IVP
from probnum.diffeq import lotkavolterra
from probnum.diffeq import logistic, probsolve_ivp


We start with the Lotka-Volterra equations.


def f(t, y):
    y1, y2 = y
    return np.array([0.5 * y1 - 0.05 * y1 * y2, -0.5 * y2 + 0.05 * y1 * y2])

def df(t, y):
    y1, y2 = y
    return np.array([[0.5 - 0.05 * y2, -0.05 * y1], [0.05 * y2, -0.5 + 0.05 * y1]])

t0 = 0.
tmax = 20.
y0 = np.array([20, 20])

The EK0 ODE filter is quite fast and flexible, but does not yield accurate uncertainty estimates. We see below that the uncertainty increases monotonously, independent on the actual peaks and valleys of the ODE solution.

It is not surprising that the EK0 solution is agnostic of the trajectory. In fact, the covariance of the EK0 filter is independent of the data and as such, we cannot expect it to return useful uncertainty.

sol = probsolve_ivp(f, t0, tmax, y0, df=df, step=0.1, adaptive=False)
means, stds = sol.states.mean, sol.states.std

fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, sharex=True)
ax1.plot(sol.locations, means[:, 0], label="x1")
ax1.plot(sol.locations, means[:, 1], label="x2")
ax2.fill_between(sol.locations, stds[:, 0], alpha=0.25, label="x1-unc.")
ax2.fill_between(sol.locations, stds[:, 1], alpha=0.25, label="x2-unc.")
fig.suptitle("EK0 Solution")

Notice that the uncertainties are aware of the peaks and valleys. They even know whether the peak is flat (rhs of blue, lhs of orange; smaller ones of the blue peaks) or steep. On top of that, they increase over time.

For both we can also sample from the solution. Let us compute a low-res solution (so the samples actually look different from each other).

Beware, because large numbers of samples can take a long time to compute.

sol = probsolve_ivp(f, t0, tmax, y0, df=df, step=0.5, algo_order=1, adaptive=False)
num_samples = 20
locations = np.arange(0., 22.1, 0.1)
samples = sol.sample(t = locations, size=num_samples)

solution = sol(locations)
means = solution.mean
stds = solution.std
fig, (ax1) = plt.subplots(nrows=1, ncols=1, sharex=True)
for sample in samples:
    ax1.plot(locations, sample[:, 0], ":", color="C0", linewidth="1", alpha=1)
    ax1.plot(locations, sample[:, 1], ":", color="C1", linewidth="1", alpha=1)
ax1.plot(locations, means[:, 0], label="x1", color="k")
ax1.plot(locations, means[:, 1], label="x2", color="k")
    means[:, 0] - 2 * stds[:, 0],
    means[:, 0] + 2 * stds[:, 0],
    means[:, 1] - 2 * stds[:, 1],
    means[:, 1] + 2 * stds[:, 1],

ax1.plot(sol.locations, sol.states.mean, "o")
fig.suptitle(f"$N={num_samples}$ Samples From an EK0 Solution")
[ ]: