Stochasticity and Uncertainty › Worked Example: Uncertainty in a Noisy Nonlinear System
World Models HandbookAll chapters Stochasticity and Uncertainty › Worked Example: Uncertainty in a Noisy Nonlinear SystemOne section · Sections in this chapter8 / 11 ↓
Reading tools Definitions & inline links
Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.
Article links
Make inline references clickable
Let us build these ideas in one small, executable study. The mean transition is with additive noise of standard deviation . Fitting a high-degree polynomial only inside a limited interval makes out-of-interval behavior underdetermined and intentionally exposes ensemble disagreement. The construction examines data-dependent disagreement, recursive propagation, a bimodal transition, and empirical interval coverage. Its numerical conclusions apply to this fixed seed, model class, and sampling design.
First, the system and the data.
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm
def true_transition(s, a, rng=None, noise=0.05):
"""Stochastic one-step dynamics: nonlinear mean, additive noise."""
mean = np.sin(s) + 0.4 * a
if rng is None:
return mean
return mean + noise * rng.standard_normal(np.shape(mean))
rng = np.random.default_rng(0)
def collect_data(n):
s = rng.uniform(-1.5, 1.5, size=n)
a = rng.uniform(-1.0, 1.0, size=n)
return s, a, true_transition(s, a, rng)
s_small, a_small, y_small = collect_data(25)
s_large, a_large, y_large = collect_data(500)Both datasets are drawn from the same state interval, from to ; the intended design difference is sample count, although the independently drawn finite datasets also differ in their realized states, actions, and transition noise. In this polynomial ensemble, additional interior data should reduce bootstrap disagreement most clearly inside the interval. That is an experiment-specific expectation, not a general locality law: data can also reduce exterior uncertainty when shared parameters or structure justify the transfer. The deliberately high-degree polynomial basis is transparent, fast, and prone to unstable extrapolation in this setup.
def poly_features(s, a, degree=7):
"""Polynomial basis in state, plus a linear action term."""
s = np.asarray(s, dtype=float)
a = np.asarray(a, dtype=float)
cols = [np.ones_like(s)] + [s**d for d in range(1, degree + 1)] + [a]
return np.column_stack(cols)
def fit_ensemble(s, a, y, n_members=25):
"""Bootstrap an ensemble of ridge-regularized polynomial models."""
X = poly_features(s, a)
models = []
for _ in range(n_members):
idx = rng.integers(0, len(y), size=len(y)) # bootstrap resample
lam = 10 ** rng.uniform(-4, -2) # jitter regularization
Xk, yk = X[idx], y[idx]
w = np.linalg.solve(Xk.T @ Xk + lam * np.eye(X.shape[1]), Xk.T @ yk)
models.append(w)
return models
def predict(models, s, a):
"""Per-member predictions, shape (n_members, n_points)."""
X = poly_features(s, a)
return np.stack([X @ w for w in models])
ens_small = fit_ensemble(s_small, a_small, y_small)
ens_large = fit_ensemble(s_large, a_large, y_large)Each member sees a different bootstrap resample and regularization strength. The regularization jitter adds diversity but is not a calibrated posterior over hyperparameters; it is an explicit experimental perturbation.
Reading ensemble disagreement as epistemic uncertainty
We evaluate both ensembles on a grid of states that extends well beyond the training region, with the action held at zero, and summarize each ensemble by its mean and standard deviation across members.
grid = np.linspace(-3.0, 3.0, 200)
a_zero = np.zeros_like(grid)
preds_small = predict(ens_small, grid, a_zero)
preds_large = predict(ens_large, grid, a_zero)
mean_small, std_small = preds_small.mean(0), preds_small.std(0)
mean_large, std_large = preds_large.mean(0), preds_large.std(0)
true_curve = np.sin(grid) # true mean dynamics with a = 0
inside = np.abs(grid) < 1.5
print(
f"25 transitions: mean spread inside data region: {std_small[inside].mean():.4f}"
)
print(
f"25 transitions: mean spread outside data region: {std_small[~inside].mean():.4f}"
)
print(
f"500 transitions: mean spread inside data region: {std_large[inside].mean():.4f}"
)
print(
f"500 transitions: mean spread outside data region: {std_large[~inside].mean():.4f}"
)25 transitions: mean spread inside data region: 0.1136 25 transitions: mean spread outside data region: 80.0959 500 transitions: mean spread inside data region: 0.0063 500 transitions: mean spread outside data region: 4.8433
For this seed and polynomial ensemble, the printed values quantify how member disagreement changes with dataset size inside and outside the sampled interval. The comparison is evidence about this construction, not a universal map from input support to epistemic uncertainty. In models with shared global parameters or structural priors, interior data can justifiably reduce uncertainty elsewhere.
The teal traces are the fitted member functions. Shared data-dependent symmetric-log limits include every shown trace and both ±2-standard-deviation bands with five percent padding; the bands summarize member disagreement rather than calibrated predictive coverage. Outside the sampled state interval, several high-degree polynomial fits diverge while remaining unconstrained by the training loss there. A planner can use that pattern as a warning signal after validating that ensemble spread tracks error for the deployment regime.
Watching uncertainty compound over rollouts
Next we iterate each fitted member deterministically from the same initial state under one open-loop action sequence. The resulting spread measures between-member divergence only; it omits the system's within-model noise and therefore is not the total predictive distribution.
horizon = 12
actions = 0.8 * np.sin(0.7 * np.arange(horizon))
s0 = 0.3
def rollout(w, s0, actions):
s = s0
traj = [s0]
for a in actions:
s = (poly_features(np.array([s]), np.array([a])) @ w).item()
traj.append(s)
return np.array(traj)
def true_rollout(s0, actions):
s = s0
traj = [s0]
for a in actions:
s = float(np.sin(s) + 0.4 * a)
traj.append(s)
return np.array(traj)
trajs = np.stack([rollout(w, s0, actions) for w in ens_small])
truth_roll = true_rollout(s0, actions)
spread = trajs.std(axis=0)
steps = np.arange(horizon + 1)
spread_ratio = spread[-1] / spread[1]
print(f"Ensemble spread at step 1: {spread[1]:.4f}")
print(f"Ensemble spread at step {horizon}: {spread[-1]:.4f}")
print(f"Step-{horizon}/step-1 spread ratio: {spread_ratio:.3f}")Ensemble spread at step 1: 0.0193 Ensemble spread at step 12: 0.2334 Step-12/step-1 spread ratio: 12.093
For this fixed ensemble and action sequence, the printed ratio is about . The recursive rollout feeds each member's prediction back as its next input, so local differences can be amplified. The magnitude is specific to the fitted models and should not be generalized to every system.
Short horizons, feedback, and frequent replanning can reduce exposure to accumulated model error when state estimates are adequate, the closed loop is stabilizing, model and optimization errors remain controlled, and the horizon suits the dynamics. Receding-horizon control replans from the latest state estimate, but it does not reset estimation error, model bias, or optimization failure. Recursive prediction permits error growth in locally expanding regimes; contractive dynamics can damp it.
A multimodal transition that a Gaussian cannot represent
Now we switch from epistemic to aleatoric structure. In this synthetic transition, an object at state receives an impulse of with equal probability plus Gaussian noise. The two branches put very little density in the middle, although Gaussian noise gives nonzero density at every finite location. We sample this fixed construction and fit two candidate distributions.
n = 4000
s_fixed, a_fixed = 0.1, 0.0
impulse = np.where(rng.random(n) < 0.5, 0.6, -0.6)
y_multi = s_fixed + impulse + 0.03 * rng.standard_normal(n)
## Candidate 1: a single Gaussian
gauss_mean, gauss_std = y_multi.mean(), y_multi.std()
## Candidate 2: a 2-component Gaussian mixture fit by EM
def fit_gmm_1d(y, n_iter=300):
mu = np.array([-0.2, 0.6])
var = np.array([0.05, 0.05])
pi = np.array([0.5, 0.5])
for _ in range(n_iter):
logp = (
np.log(pi)[:, None]
- 0.5 * np.log(2 * np.pi * var)[:, None]
- (y[None] - mu[:, None]) ** 2 / (2 * var[:, None])
)
logp -= logp.max(axis=0, keepdims=True)
r = np.exp(logp)
r /= r.sum(axis=0, keepdims=True)
nk = r.sum(axis=1)
pi = nk / len(y)
mu = (r @ y) / nk
var = np.maximum((r @ y**2) / nk - mu**2, 1e-6)
return pi, mu, var
pi_fit, mu_fit, var_fit = fit_gmm_1d(y_multi)
print(f"Single Gaussian: mean = {gauss_mean:.3f}, std = {gauss_std:.3f}")
print(f"Mixture weights: {pi_fit.round(3)}")
print(f"Mixture means: {mu_fit.round(3)}")
print(f"Mixture std devs: {np.sqrt(var_fit).round(3)}")Single Gaussian: mean = 0.096, std = 0.601 Mixture weights: [0.503 0.496] Mixture means: [-0.5 0.7] Mixture std devs: [0.03 0.029]
For this sample, the mixture fit places components near with near-equal weights. The single Gaussian has approximately the correct mean and variance but allocates its highest density near the low-density middle. Its error is distributional shape: it assigns much more mass between the branches than the generating mixture does.
This is the unimodal-distribution trap; its mean also illustrates point averaging. In this construction, the mixture places density near both sampled branches, while the Gaussian mode lies near the middle. Before assuming optimization or data are the only problems, check whether the output family can represent the observed conditional shape.
Measuring calibration on held-out transitions
Finally, we measure empirical coverage on one fixed 300-point test sample spanning the training interval and extrapolation. We compare intervals using the known noise standard deviation with intervals that also add member-mean variance under a Gaussian approximation. Both use the same ensemble mean. The comparison isolates the numerical effect of that added variance term in this construction; it does not establish universal calibration.
n_test = 300
s_test = rng.uniform(-2.5, 2.5, n_test)
a_test = rng.uniform(-1.0, 1.0, n_test)
y_test = true_transition(s_test, a_test, rng)
preds_test = predict(ens_small, s_test, a_test)
mean_ep = preds_test.mean(0)
alea_std = 0.05
std_single = np.full(n_test, alea_std)
std_ensemble = np.sqrt(alea_std**2 + preds_test.var(0))
levels = np.linspace(0.1, 0.9, 9)
def coverage(y, mean, std, level):
z = norm.ppf((1 + level) / 2)
return float(np.mean(np.abs(y - mean) <= z * std))
cov_single = np.array(
[coverage(y_test, mean_ep, std_single, p) for p in levels]
)
cov_ensemble = np.array(
[coverage(y_test, mean_ep, std_ensemble, p) for p in levels]
)
print(
f"Nominal 50% coverage: aleatoric-only = {cov_single[4]:.2f}, with epistemic = {cov_ensemble[4]:.2f}"
)
print(
f"Nominal 90% coverage: aleatoric-only = {cov_single[-1]:.2f}, with epistemic = {cov_ensemble[-1]:.2f}"
)Nominal 50% coverage: aleatoric-only = 0.28, with epistemic = 0.81 Nominal 90% coverage: aleatoric-only = 0.51, with epistemic = 0.98
The variance sum follows the law of total variance for an equal-weight mixture of members, while converting that variance to Gaussian intervals is an additional approximation. The identity decomposes predictive variance, not every difference in distributional shape. In this seed, the added member variance improves mean absolute coverage error but overcovers at several nominal levels. Bayesian models, quantile regression, and conformal prediction offer different predictive distributions or calibrated sets; they do not all identify the same aleatoric/epistemic decomposition.