Skip to content

Uncertainty quantification

Intervals for identified physics (pinnrl.uq): a Gaussian (Laplace) posterior from the surrogate's Fisher information, the region of equivalence sampled by an exploratory particle swarm, ensembles summarized with a dispersion diagnostic, loss weights from the noise model, and the calibration metrics that judge them. See the theory pages on the Bayesian view and on equivalence regions.

laplace_posterior 

laplace_posterior(
    pde,
    model,
    parameters,
    sigma,
    x=None,
    t=None,
    **sensitivity_kwargs,
)

Gaussian posterior over parameters at the trained estimate.

The covariance is the inverse Fisher information of the surrogate's sensitivities at the observation points (pde.observation_data by default, or x, t). This is the Laplace approximation restricted to the parameters with the generalized Gauss–Newton curvature; it coincides with the Cramér–Rao bound at the estimate.

Source code in pinnrl/uq/laplace.py 
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
def laplace_posterior(
    pde: Any,
    model: torch.nn.Module,
    parameters: Sequence[str],
    sigma: float | torch.Tensor,
    x: torch.Tensor | None = None,
    t: torch.Tensor | None = None,
    **sensitivity_kwargs: Any,
) -> LaplacePosterior:
    """Gaussian posterior over ``parameters`` at the trained estimate.

    The covariance is the inverse Fisher information of the surrogate's
    sensitivities at the observation points (``pde.observation_data`` by
    default, or ``x, t``). This is the Laplace approximation restricted to the
    parameters with the generalized Gauss–Newton curvature; it coincides with
    the Cramér–Rao bound at the estimate.
    """
    if x is None:
        obs = pde.observation_data
        if obs is None:
            raise ValueError(
                "laplace_posterior needs observation points: pde.observation_data or x, t"
            )
        x, t = obs["x"], obs["t"]
    sens = surrogate_sensitivities(pde, model, parameters, x, t, **sensitivity_kwargs)
    summary = fisher_information(sens.jacobian, sigma, parameters)
    values = pde.get_trainable_parameter_values()
    mean = torch.tensor(
        [values.get(n, _scalar(pde.get_parameter(n))) for n in parameters], dtype=torch.float64
    )
    return LaplacePosterior(
        parameters=tuple(parameters),
        mean=mean,
        covariance=summary.covariance,
        std=summary.std,
        correlation=summary.correlation,
    )

LaplacePosterior dataclass 

LaplacePosterior(
    parameters, mean, covariance, std, correlation
)

A Gaussian posterior over identified parameters.

Attributes:

Name Type Description
parameters tuple[str, ...]

Names.
mean Tensor

The point estimate (the trained values).
covariance Tensor

(p, p).
std Tensor

Marginal standard deviations.
correlation Tensor

(p, p).

interval 

interval(level=0.9)

Central interval per parameter, (2, p).

Source code in pinnrl/uq/laplace.py 
63
64
65
66
67
68
def interval(self, level: float = 0.9) -> torch.Tensor:
    """Central interval per parameter, ``(2, p)``."""
    from pinnrl.uq.metrics import _normal_quantile

    z = _normal_quantile(0.5 + level / 2)
    return torch.stack([self.mean - z * self.std, self.mean + z * self.std])

sample 

sample(n, seed=0)

n draws, (n, p).

Source code in pinnrl/uq/laplace.py 
70
71
72
73
74
75
def sample(self, n: int, seed: int = 0) -> torch.Tensor:
    """``n`` draws, ``(n, p)``."""
    gen = torch.Generator().manual_seed(seed)
    L = torch.linalg.cholesky(self.covariance)
    z = torch.randn(n, self.mean.numel(), generator=gen, dtype=self.covariance.dtype)
    return self.mean.double() + z @ L.T

likelihood_weights 

likelihood_weights(
    n_data,
    sigma_data,
    n_residual,
    sigma_residual=1.0,
    *,
    boundary=10.0,
    initial=10.0,
)

Loss weights that make the mean-squared loss a negative log-likelihood up to scale.

The loss uses means, so a term with n points and noise sigma contributes n / (2 sigma^2) times its mean square to -log L. With the residual weight fixed at 1, the data weight is (n_data / sigma_data^2) / (n_residual / sigma_residual^2), where sigma_residual is the model-error scale admitted for the equation.

Source code in pinnrl/uq/laplace.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
def likelihood_weights(
    n_data: int,
    sigma_data: float,
    n_residual: int,
    sigma_residual: float = 1.0,
    *,
    boundary: float = 10.0,
    initial: float = 10.0,
) -> dict[str, float]:
    """Loss weights that make the mean-squared loss a negative log-likelihood up to scale.

    The loss uses means, so a term with ``n`` points and noise ``sigma``
    contributes ``n / (2 sigma^2)`` times its mean square to ``-log L``. With
    the residual weight fixed at 1, the data weight is
    ``(n_data / sigma_data^2) / (n_residual / sigma_residual^2)``, where
    ``sigma_residual`` is the model-error scale admitted for the equation.
    """
    data = (n_data / sigma_data**2) / (n_residual / sigma_residual**2)
    return {"residual": 1.0, "boundary": boundary, "initial": initial, "data": float(data)}

sample_equivalence_region 

sample_equivalence_region(
    misfit,
    bounds,
    tolerance=None,
    *,
    sigma=None,
    level=0.9,
    n_particles=40,
    n_iterations=100,
    seed=0,
    exploration=0.9,
    stop_on_collapse=0.02,
)

Sample {theta : misfit(theta) <= tolerance} with a cloud particle swarm.

Parameters:

Name Type Description Default
misfit Callable[[ndarray], ndarray]

Vectorized misfit, (n, p) -> (n,). For Gaussian noise use the sum of squared residuals over sigma^2 and a tolerance of about N (1 + epsilon); see :func:misfit_tolerance.

 required
bounds Sequence[tuple[float, float]]

(low, high) per parameter: the prior range, which the swarm explores and particles are reflected into.

 required
tolerance float | None

Accept positions with misfit at or below it. When None, the likelihood-ratio region is used instead: positions within sigma^2 chi2_p(level) of the best misfit found, which for Gaussian noise is the asymptotic level confidence region and is directly comparable with the Fisher ellipse at that level.

 None
sigma float | None

Noise standard deviation, needed when tolerance is None and the misfit is a plain sum of squares (use 1 if the misfit is already divided by sigma^2).

 None
level float

Confidence level of the likelihood-ratio region.

 0.9
n_particles int

Swarm size.

 40
n_iterations int

Iterations; every evaluation is a candidate sample.

 100
seed int

Seed of the initial positions and of the per-particle parameters.

 0
exploration float

Fraction of the way to the second-order stability border at which each particle's acceleration is drawn (close to 1 is the exploratory setting the papers recommend).

 0.9
stop_on_collapse float

Stop when the dispersion falls below this fraction of its initial value.

 0.02
Source code in pinnrl/uq/swarm.py 
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
def sample_equivalence_region(
    misfit: Callable[[np.ndarray], np.ndarray],
    bounds: Sequence[tuple[float, float]],
    tolerance: float | None = None,
    *,
    sigma: float | None = None,
    level: float = 0.9,
    n_particles: int = 40,
    n_iterations: int = 100,
    seed: int = 0,
    exploration: float = 0.9,
    stop_on_collapse: float = 0.02,
) -> SwarmSample:
    """Sample ``{theta : misfit(theta) <= tolerance}`` with a cloud particle swarm.

    Args:
        misfit: Vectorized misfit, ``(n, p) -> (n,)``. For Gaussian noise use
            the sum of squared residuals over ``sigma^2`` and a tolerance of
            about ``N (1 + epsilon)``; see :func:`misfit_tolerance`.
        bounds: ``(low, high)`` per parameter: the prior range, which the
            swarm explores and particles are reflected into.
        tolerance: Accept positions with misfit at or below it. When ``None``,
            the likelihood-ratio region is used instead: positions within
            ``sigma^2 chi2_p(level)`` of the best misfit found, which for
            Gaussian noise is the asymptotic ``level`` confidence region and
            is directly comparable with the Fisher ellipse at that level.
        sigma: Noise standard deviation, needed when ``tolerance`` is ``None``
            and the misfit is a plain sum of squares (use ``1`` if the misfit
            is already divided by ``sigma^2``).
        level: Confidence level of the likelihood-ratio region.
        n_particles: Swarm size.
        n_iterations: Iterations; every evaluation is a candidate sample.
        seed: Seed of the initial positions and of the per-particle parameters.
        exploration: Fraction of the way to the second-order stability border
            at which each particle's acceleration is drawn (close to 1 is the
            exploratory setting the papers recommend).
        stop_on_collapse: Stop when the dispersion falls below this fraction
            of its initial value.
    """
    rng = np.random.default_rng(seed)
    lo = np.asarray([b[0] for b in bounds], dtype=float)
    hi = np.asarray([b[1] for b in bounds], dtype=float)
    p = lo.size
    x = lo + (hi - lo) * rng.random((n_particles, p))
    v = np.zeros_like(x)
    # Cloud parameters: inertia in [0.5, 0.95], total acceleration a fraction of
    # the way to the variance-stability border, time step in [0.8, 1.2].
    omega = rng.uniform(0.5, 0.95, n_particles)
    phi_max = 12 * (1 - omega**2) / (7 - 5 * omega)
    phi = exploration * phi_max * rng.uniform(0.8, 1.0, n_particles)
    dt = rng.uniform(0.8, 1.2, n_particles)

    f = np.asarray(misfit(x), dtype=float)
    local_best, local_f = x.copy(), f.copy()
    g_idx = int(np.argmin(f))
    global_best, global_f = x[g_idx].copy(), float(f[g_idx])
    positions, misfits = [x.copy()], [f.copy()]
    centre = np.median(x, axis=0)
    d0 = float(np.median(np.linalg.norm(x - centre, axis=1))) or 1.0
    dispersion = [1.0]
    for _ in range(1, n_iterations):
        r1, r2 = rng.random((n_particles, p)), rng.random((n_particles, p))
        a_g = phi[:, None] * r1  # each particle's accelerations sum to phi on average
        a_l = phi[:, None] * r2
        # Generalized PSO with time step dt: damping (1 - omega) dt, stiffness phi dt^2.
        v = (1 - (1 - omega[:, None]) * dt[:, None]) * v + dt[:, None] * (
            a_g * (global_best - x) + a_l * (local_best - x)
        )
        x = x + dt[:, None] * v
        # Reflect into the bounds.
        below, above = x < lo, x > hi
        x = np.where(below, 2 * lo - x, x)
        x = np.where(above, 2 * hi - x, x)
        x = np.clip(x, lo, hi)
        v = np.where(below | above, -v, v)
        f = np.asarray(misfit(x), dtype=float)
        improved = f < local_f
        local_best[improved], local_f[improved] = x[improved], f[improved]
        i = int(np.argmin(local_f))
        if local_f[i] < global_f:
            global_best, global_f = local_best[i].copy(), float(local_f[i])
        positions.append(x.copy())
        misfits.append(f.copy())
        centre = np.median(x, axis=0)
        dispersion.append(float(np.median(np.linalg.norm(x - centre, axis=1))) / d0)
        if dispersion[-1] < stop_on_collapse:
            break
    P = np.concatenate(positions)
    M = np.concatenate(misfits)
    if tolerance is None:
        if sigma is None:
            raise ValueError("give a tolerance, or sigma for the likelihood-ratio region")
        tolerance = global_f + sigma**2 * chi2_quantile(p, level)
    return SwarmSample(
        positions=P,
        misfits=M,
        accepted=M <= tolerance,
        best=global_best,
        best_misfit=global_f,
        dispersion=np.asarray(dispersion),
        tolerance=tolerance,
    )

SwarmSample dataclass 

SwarmSample(
    positions,
    misfits,
    accepted,
    best,
    best_misfit,
    dispersion,
    tolerance,
)

The outcome of a swarm run.

Attributes:

Name Type Description
positions ndarray

Every evaluated position, (n_iter * n_particles, p).
misfits ndarray

Their misfits.
accepted ndarray

Mask of the positions with misfit <= tol.
best ndarray

The lowest-misfit position.
best_misfit float

Its misfit.
dispersion ndarray

Per iteration, the median distance of the particles to their centre of gravity relative to iteration 0 (1.0 at the start; a value near 0 means the swarm collapsed).
tolerance float

The tolerance used.

samples property 

samples

The accepted positions, the sample of the region of equivalence.

interval 

interval(level=0.9)

Central empirical interval per parameter from the accepted samples, (2, p).

Source code in pinnrl/uq/swarm.py 
55
56
57
58
def interval(self, level: float = 0.9) -> np.ndarray:
    """Central empirical interval per parameter from the accepted samples, ``(2, p)``."""
    s = self.samples
    return np.percentile(s, [50 * (1 - level), 50 * (1 + level)], axis=0)

ranges 

ranges()

Min and max per parameter over the accepted samples, (2, p).

Source code in pinnrl/uq/swarm.py 
60
61
62
63
def ranges(self) -> np.ndarray:
    """Min and max per parameter over the accepted samples, ``(2, p)``."""
    s = self.samples
    return np.stack([s.min(axis=0), s.max(axis=0)])

misfit_tolerance 

misfit_tolerance(n_observations, sigma, epsilon=0.1)

Tolerance on the squared misfit sum (y - u)^2 for Gaussian noise.

The expected squared misfit of the true parameters is N sigma^2; the region of equivalence is taken as the set within (1 + epsilon) of it.

Source code in pinnrl/uq/swarm.py 
197
198
199
200
201
202
203
def misfit_tolerance(n_observations: int, sigma: float, epsilon: float = 0.1) -> float:
    """Tolerance on the squared misfit ``sum (y - u)^2`` for Gaussian noise.

    The expected squared misfit of the true parameters is ``N sigma^2``; the
    region of equivalence is taken as the set within ``(1 + epsilon)`` of it.
    """
    return (1 + epsilon) * n_observations * sigma**2

second_order_stable 

second_order_stable(omega, phi)

Whether a PSO parameter pair is inside the second-order (variance) stability region.

For the classical PSO with uniformly distributed accelerations the mean trajectory is stable for |omega| < 1 and 0 < phi < 2 (omega + 1), and the variance for phi < 12 (1 - omega^2) / (7 - 5 omega) (Poli, 2009; García-Gonzalo & Fernández-Martínez, 2014, for general distributions). phi is the mean total acceleration (a_g + a_l) / 2.

Source code in pinnrl/uq/swarm.py 
66
67
68
69
70
71
72
73
74
75
def second_order_stable(omega: float, phi: float) -> bool:
    """Whether a PSO parameter pair is inside the second-order (variance) stability region.

    For the classical PSO with uniformly distributed accelerations the mean
    trajectory is stable for ``|omega| < 1`` and ``0 < phi < 2 (omega + 1)``, and
    the variance for ``phi < 12 (1 - omega^2) / (7 - 5 omega)`` (Poli, 2009;
    García-Gonzalo & Fernández-Martínez, 2014, for general distributions).
    ``phi`` is the mean total acceleration ``(a_g + a_l) / 2``.
    """
    return abs(omega) < 1 and 0 < phi < 12 * (1 - omega**2) / (7 - 5 * omega)

ensemble_summary 

ensemble_summary(estimates)

Summarize the estimates of an ensemble of independently trained members.

Source code in pinnrl/uq/__init__.py 
71
72
73
74
75
76
77
78
79
80
81
def ensemble_summary(estimates: Sequence[Sequence[float]]) -> EnsembleSummary:
    """Summarize the estimates of an ensemble of independently trained members."""
    e = np.atleast_2d(np.asarray(estimates, dtype=float))
    med = np.median(e, axis=0)
    scale = float(np.mean(np.abs(med))) or 1.0
    return EnsembleSummary(
        estimates=e,
        mean=e.mean(axis=0),
        std=e.std(axis=0, ddof=1) if e.shape[0] > 1 else np.zeros(e.shape[1]),
        dispersion=float(np.median(np.linalg.norm(e - med, axis=1))) / scale,
    )

EnsembleSummary dataclass 

EnsembleSummary(estimates, mean, std, dispersion)

Point estimates of an ensemble and what they imply.

Attributes:

Name Type Description
estimates ndarray

(n_members, p).
mean ndarray

Per parameter.
std ndarray

Per parameter (the ensemble spread).
dispersion float

Median distance of the members to their median, relative to the mean absolute value of the parameters (a collapse diagnostic).

interval 

interval(level=0.9, gaussian=False)

Central interval per parameter, (2, p): empirical, or Gaussian from mean and std.

Source code in pinnrl/uq/__init__.py 
61
62
63
64
65
66
67
68
def interval(self, level: float = 0.9, gaussian: bool = False) -> np.ndarray:
    """Central interval per parameter, ``(2, p)``: empirical, or Gaussian from mean and std."""
    if gaussian:
        from pinnrl.uq.metrics import _normal_quantile

        z = _normal_quantile(0.5 + level / 2)
        return np.stack([self.mean - z * self.std, self.mean + z * self.std])
    return np.percentile(self.estimates, [50 * (1 - level), 50 * (1 + level)], axis=0)

coverage 

coverage(lower, upper, truth)

Fraction of intervals [lower, upper] that contain the truth.

Source code in pinnrl/uq/metrics.py 
11
12
13
14
def coverage(lower: Sequence[float], upper: Sequence[float], truth: Sequence[float]) -> float:
    """Fraction of intervals ``[lower, upper]`` that contain the truth."""
    lo, hi, t = (np.asarray(v, dtype=float) for v in (lower, upper, truth))
    return float(np.mean((lo <= t) & (t <= hi)))

mean_width 

mean_width(lower, upper)

Mean width of the intervals.

Source code in pinnrl/uq/metrics.py 
17
18
19
def mean_width(lower: Sequence[float], upper: Sequence[float]) -> float:
    """Mean width of the intervals."""
    return float(np.mean(np.asarray(upper, dtype=float) - np.asarray(lower, dtype=float)))

gaussian_interval 

gaussian_interval(mean, std, level=0.9)

Central interval of a Gaussian at the given level.

Source code in pinnrl/uq/metrics.py 
22
23
24
25
def gaussian_interval(mean: float, std: float, level: float = 0.9) -> tuple[float, float]:
    """Central interval of a Gaussian at the given level."""
    z = _normal_quantile(0.5 + level / 2)
    return mean - z * std, mean + z * std

sample_interval 

sample_interval(samples, level=0.9)

Central empirical interval of a set of samples at the given level.

Source code in pinnrl/uq/metrics.py 
28
29
30
31
32
def sample_interval(samples: Sequence[float], level: float = 0.9) -> tuple[float, float]:
    """Central empirical interval of a set of samples at the given level."""
    s = np.asarray(samples, dtype=float)
    lo, hi = np.percentile(s, [50 * (1 - level), 50 * (1 + level)])
    return float(lo), float(hi)

crps_gaussian 

crps_gaussian(mean, std, truth)

Continuous ranked probability score of a Gaussian forecast (Gneiting & Raftery, 2007).

CRPS = sigma [ z (2 Phi(z) - 1) + 2 phi(z) - 1/sqrt(pi) ] with z = (y - mu) / sigma. Lower is better; it rewards both calibration and sharpness.

Source code in pinnrl/uq/metrics.py 
35
36
37
38
39
40
41
42
43
44
45
46
def crps_gaussian(mean: float, std: float, truth: float) -> float:
    """Continuous ranked probability score of a Gaussian forecast (Gneiting & Raftery, 2007).

    ``CRPS = sigma [ z (2 Phi(z) - 1) + 2 phi(z) - 1/sqrt(pi) ]`` with ``z = (y - mu) / sigma``.
    Lower is better; it rewards both calibration and sharpness.
    """
    if std <= 0:
        return abs(truth - mean)
    z = (truth - mean) / std
    phi = math.exp(-0.5 * z * z) / math.sqrt(2 * math.pi)
    Phi = 0.5 * (1 + math.erf(z / math.sqrt(2)))
    return std * (z * (2 * Phi - 1) + 2 * phi - 1 / math.sqrt(math.pi))

crps_samples 

crps_samples(samples, truth)

CRPS of an empirical forecast: E|X - y| - E|X - X'| / 2.

Source code in pinnrl/uq/metrics.py 
49
50
51
52
53
54
def crps_samples(samples: Sequence[float], truth: float) -> float:
    """CRPS of an empirical forecast: ``E|X - y| - E|X - X'| / 2``."""
    s = np.asarray(samples, dtype=float)
    term1 = float(np.mean(np.abs(s - truth)))
    term2 = float(np.mean(np.abs(s[:, None] - s[None, :])))
    return term1 - 0.5 * term2