<jupyter_start><jupyter_text>Quantify approximation error and calibration of a PN MethodThis notebook explains how to use the functions in probnum-evaluation to quantify the approximation error as well as the uncertainty calibration of the output of a probabilistic numerical method.<jupyter_code>from probnum import quad, diffeq
import numpy as np

np.random.seed(1)
from probnumeval import multivariate, timeseries<jupyter_output><empty_output><jupyter_text>One-dimensional exampleThe first example shall be Bayesian quadrature, which computes a probabilistic numerical solution of (e.g.) the integral$$F = \int_0^\pi \sin(x) \,\text{d}x.$$The true solution of this integral is $F = 2.0$. Begin by setting the problem up.<jupyter_code>input_dim = 1
domain = (0, np.pi)

# Define the integrand
def f(x):
    return np.sin(x)


# True solution of \int_0^\pi sin(x) dx is 2.0
sol = 2.0<jupyter_output><empty_output><jupyter_text>The solution of this integration problem is computed with `probnum.quad` (https://probnum.readthedocs.io/en/latest/public_api/quad.html). Let us only use few function evaluations, so the solution is sufficiently inaccurate and the covariance of the output will be $\gg 0$.<jupyter_code>F, _ = quad.bayesquad(fun=f, input_dim=input_dim, domain=domain, nevals=5)<jupyter_output><empty_output><jupyter_text>In `ProbNum-Evaluation`, finite-dimensional outputs (such as $F$) are quantified with `probnumeval.multivariate`. Among other things, this module contains error measures such as root mean-square errors (RMSEs), mean absolute errors, or their relative counterparts. For instance, the RMSE and relative RMSE are computed as follows.<jupyter_code>rmse = multivariate.rmse(F.mean, reference_solution=sol)
relative_rmse = multivariate.relative_rmse(F.mean, reference_solution=sol)
print(rmse, relative_rmse)<jupyter_output>0.06127865361488327 0.030639326807441636<jupyter_text>We can also check how well the covariance quantifies the approximation error -- for instance, with the averaged normalised estimation error squared (ANEES). A value close to 1 is good.<jupyter_code>anees = multivariate.anees(F, reference_solution=sol)
print(anees)<jupyter_output>0.12591126045245504<jupyter_text>Time-series exampleThe same functionality can be used for a more complex posterior.The next example is a probabilistic numerical ODE solution of $$\dot y(t) = y(t) * (1 - y(t)), \quad y(0) = 0.5$$for $t\in (0.0, 5.0)$. This logistic ODE has the analytic solution$$y(t) = \frac{y_0 e^t}{1 + y_0 (e^t - 1)}.$$Let us begin by setting up the problem.<jupyter_code>def f(t, x):
    return x * (1 - x)


y0 = np.array([0.5])


t0 = 0.0
tmax = 5.0


def sol(t):
    return (y0 * np.exp(t)) / (1 + y0 * (np.exp(t) - 1))<jupyter_output><empty_output><jupyter_text>In `ProbNum`, ODEs are solved with `probnum.diffeq`. The outputs are usually `TimeSeriesPosterior`, i.e. functions that carry the states of the solution at a grid, and can also be called at any input.To evaluate the posterior, let us include the zeroth point of the domain (which is error-free by construction, and has covariance zero). If we did not do this, computation of calibration indices may be subject to numerical instability.<jupyter_code>pnsol = diffeq.probsolve_ivp(f, y0=y0, t0=t0, tmax=tmax)

approx_sol = pnsol.states[1:]
true_sol = sol(pnsol.locations[1:])[:, None]<jupyter_output><empty_output><jupyter_text>On the given set of evaluations of the solution, we can compute any error measure. This time, we choose the mean absolute error (MAE) and its relative counterpart, the relative MAE.<jupyter_code>mae = multivariate.mae(approx_sol.mean, true_sol)
relative_mae = multivariate.relative_mae(approx_sol.mean, true_sol)
print(mae, relative_mae)<jupyter_output>0.0006241640303712251 0.0006614154047129657<jupyter_text>Since covariances of ODE solutions can be ill-conditioned, we make sure that all the covariance calculations use pseudo-inverses.This step is totally optional.<jupyter_code>from probnumeval import config

config.set_covariance_inversion_parameters(
    strategy="pinv", symmetrize=True, damping=0.0
)<jupyter_output><empty_output><jupyter_text>We can also quantify the uncertainty calibration, for instance by calling the ANEES again.This time however, we opt for the non-credibility index and the inclination index.The difference to the ANEES is that NCI and II compare the calibration of the posterior to a "reference calibration", and thus respond differently to optimism and pessimism.The difference between NCI and II is that NCI is an absolute number (the lower the better), whereas the II indicates over- and underconfidence. For both values, the closer to zero the better.<jupyter_code>nci = multivariate.non_credibility_index(approx_sol, true_sol)
ii = multivariate.inclination_index(approx_sol, true_sol)
print(nci, ii)<jupyter_output>12.803778627298394 11.306190199210022<jupyter_text>All of the above could have been obtained by applying the `timeseries` module to the (callable) ODE posteriors. While the `multivariate` module expects (lists of) random variables, the `timeseries` acts directly on the ODE solution.<jupyter_code># Reshape the reference solution to match it to the ODE solution
def sol_reshaped(*args, **kwargs):
    return sol(*args, **kwargs)[:, None]


# Extract the mean function from the PN solution.
def pnsol_mean(*args, **kwargs):
    return pnsol(*args, **kwargs).mean


mae = timeseries.mae(pnsol_mean, sol_reshaped, locations=pnsol.locations[1:])
relative_mae = timeseries.relative_mae(
    pnsol_mean, sol_reshaped, locations=pnsol.locations[1:]
)

nci = timeseries.non_credibility_index(
    pnsol, sol_reshaped, locations=pnsol.locations[1:]
)
ii = timeseries.inclination_index(pnsol, sol_reshaped, locations=pnsol.locations[1:])
print(mae, relative_mae, nci, ii)<jupyter_output>0.0006241640303712251 0.0006614154047129657 12.803778627298394 11.306190199210022
{'fraction_non_alphanumeric': 0.09137055837563451, 'fraction_numerical': 0.04297800338409476, 'mean_word_length': 5.995266272189349, 'pattern_counts': {'":': 0, '<': 40, '<?xml version=': 0, '>': 40, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '97220', 'license': 'MIT', 'n_tokens_mistral': 2092, 'n_tokens_neox': 1851, 'n_words': 622, 'path': 'docs/source/tutorials/measure_error_and_calibration.ipynb', 'repo_name': 'anukaal/probnum-evaluation', 'size': '5910'}