Propagating Measurement Uncertainty Through Aggregation #

An hourly mean of sixty readings is not eight times more accurate than one reading, and treating it as though it were is how a network’s published figures come to claim a precision its hardware never had. The correction is not complicated — it is the standard separation of random from systematic error — but it has to be applied at every aggregation step, including the spatial ones. This guide does that for the three aggregations an environmental pipeline actually performs: over time, over sensors, and over space. It continues uncertainty quantification and data quality scoring.

Random Averages Down, Systematic Does Not #

Split every uncertainty budget into two parts before doing anything else.

Random (independent between readings): electronic noise, quantisation, short-term turbulence. Averaging n readings reduces this by a factor of sqrt(n), because the errors partially cancel.

Systematic (shared across readings): calibration offset, drift at the current age, an unmodelled humidity response. Every reading in the hour carries the same error, so the mean carries it too, unchanged.

The consequence is a floor. For a low-cost PM2.5 sensor with 1.8 µg/m³ of noise and 5.7 µg/m³ of systematic error, a single reading has a combined uncertainty near 6.0 µg/m³ and an hourly mean of sixty samples has 5.7 µg/m³ — a 4% improvement, not the 87% that sqrt(60) would suggest. Recognising that floor changes what a network is used for: more sampling does not fix it, but a co-location campaign that shrinks the calibration residual does.

Production-Ready Implementation #

# python 3.11 · numpy==1.26.4 · pandas==2.2.2
from __future__ import annotations

from dataclasses import dataclass

import numpy as np
import pandas as pd


@dataclass(frozen=True)
class Uncertainty:
    """A measurement uncertainty split into the parts that behave differently.

    random:     independent between readings; averages down as 1/sqrt(n)
    systematic: shared by every reading from this sensor; survives averaging
    """

    random: float
    systematic: float

    @property
    def combined(self) -> float:
        return (self.random ** 2 + self.systematic ** 2) ** 0.5

    def after_averaging(self, n: int) -> "Uncertainty":
        if n <= 0:
            raise ValueError("n must be positive")
        return Uncertainty(random=self.random / n ** 0.5, systematic=self.systematic)


def temporal_mean(values: np.ndarray, u: Uncertainty) -> tuple[float, Uncertainty]:
    """Mean of n readings from ONE sensor, with correctly propagated uncertainty."""
    n = int(np.count_nonzero(~np.isnan(values)))
    return float(np.nanmean(values)), u.after_averaging(n)

Combining different sensors is an inverse-variance weighted mean, which is where the quality score earns its keep — a reference instrument and a low-cost node should not contribute equally:

def inverse_variance_mean(
    values: np.ndarray, uncertainties: np.ndarray
) -> tuple[float, float]:
    """Minimum-variance combination of readings with differing uncertainties.

    Each sensor's weight is 1/u^2, so a reference instrument with a third the
    uncertainty of a low-cost node carries nine times its weight.
    Returns (weighted mean, uncertainty of that mean).
    """
    mask = np.isfinite(values) & np.isfinite(uncertainties) & (uncertainties > 0)
    v, u = values[mask], uncertainties[mask]
    if v.size == 0:
        return float("nan"), float("nan")
    w = 1.0 / u ** 2
    mean = float(np.sum(w * v) / np.sum(w))
    return mean, float(1.0 / np.sqrt(np.sum(w)))

And through IDW interpolation, where the weights are known so the propagation is exact:

def idw_with_uncertainty(
    target_xy: tuple[float, float],
    sample_xy: np.ndarray,
    values: np.ndarray,
    uncertainties: np.ndarray,
    power: float = 2.0,
) -> tuple[float, float]:
    """IDW estimate at one target point, with propagated uncertainty.

    The interpolated value is sum(w*v)/sum(w); because the weights are constants,
    its uncertainty is sqrt(sum((w*u)^2))/sum(w) — assuming sensor errors are
    independent, which is true for noise and NOT true for a shared calibration
    bias across identical hardware.
    """
    d = np.hypot(sample_xy[:, 0] - target_xy[0], sample_xy[:, 1] - target_xy[1])
    d = np.maximum(d, 1e-9)
    w = 1.0 / d ** power
    total = float(np.sum(w))
    value = float(np.sum(w * values) / total)
    u = float(np.sqrt(np.sum((w * uncertainties) ** 2)) / total)
    return value, u

The caveat in that docstring is the one that matters in practice: a fleet of identical low-cost sensors calibrated against the same reference shares a calibration bias, so their errors are not independent and the propagated figure is optimistic. Where that applies, add the shared bias back as a systematic term after the interpolation rather than pretending it averaged away.

Parameter Tuning Guide #

Read the last column as a research priority. When averaging buys 4%, the only route to a materially better figure is reducing the systematic term — better calibration, humidity correction, more frequent service — and no amount of sampling will substitute.

Verification and Testing #

Simulation is the right test here, because the true value is known only in a simulation.

# python 3.11 · numpy==1.26.4 · pytest==8.2.0
import numpy as np


def test_random_error_averages_down_and_systematic_does_not():
    rng = np.random.default_rng(42)
    truth, n_trials, n_samples = 20.0, 4000, 60
    u = Uncertainty(random=1.8, systematic=5.7)

    means = []
    for _ in range(n_trials):
        bias = rng.normal(0, u.systematic)            # one bias per sensor-hour
        readings = truth + bias + rng.normal(0, u.random, n_samples)
        means.append(readings.mean())

    observed = float(np.std(means))
    predicted = u.after_averaging(n_samples).combined
    assert abs(observed - predicted) / predicted < 0.06     # within 6%


def test_inverse_variance_beats_the_plain_mean():
    values = np.array([21.0, 19.0, 26.0])
    uncertainties = np.array([1.0, 1.0, 8.0])              # the third sensor is poor
    weighted, u = inverse_variance_mean(values, uncertainties)

    assert abs(weighted - 20.0) < abs(values.mean() - 20.0)  # closer to the reference pair
    assert u < uncertainties.min()                            # combining beats any single input


def test_uncertainty_of_a_mean_never_falls_below_the_systematic_floor():
    u = Uncertainty(random=1.8, systematic=5.7)
    assert u.after_averaging(100_000).combined >= 5.7

The first test is worth keeping in CI permanently: it is the empirical demonstration that the model is right, and it fails loudly if someone “simplifies” the propagation into a single sqrt(n) division.

Gotchas #

Dividing the whole uncertainty by sqrt(n). The single most common error in this area. It produces an hourly mean claiming 0.8 µg/m³ uncertainty from sensors that are 6 µg/m³ off, and because the number is small it is rarely questioned.

Assuming independence across identical hardware. A hundred nodes of the same model calibrated in the same campaign share a systematic error. Averaging them does not remove it, and a spatial mean across the network inherits it in full.

Propagating uncertainty through a non-linear transform linearly. If the reported value is a function of the raw measurement — a humidity correction, a log transform — the uncertainty must be scaled by the derivative at that point, not carried across unchanged.

Reporting a coverage factor without saying so. A “±12 µg/m³” figure means nothing unless the reader knows whether it is k=1 (68%) or k=2 (95%). State the coverage factor in the column description and keep the stored value at k=1.

FAQ #

Related #

• Uncertainty Quantification and Data Quality Scoring — the uncertainty budget this propagation starts from
• Computing a Composite Data Quality Score per Reading — the dimensionless companion to these physical error bars
• Ordinary Kriging vs IDW for Sparse Sensor Networks — the interpolation choice that decides whether uncertainty survives to the map