Streaming Limit Updates with Incremental Statistics

This is the recipe for keeping Individual and Moving Range control limits current on a live measurement stream by updating running statistics in constant time per point, instead of re-reading the whole history on every new observation. It is the streaming counterpart to the batch approach in vectorized control limit calculation engines: where the batch engine reduces a bounded table in one pass, a stream never ends, so recomputing from scratch grows without bound and eventually cannot keep pace with arrivals. Welford's algorithm gives a numerically stable running mean and variance, and an incremental moving range gives the dispersion estimate an I-MR chart needs — both updated with a fixed amount of work per point.

Prerequisites

  • Python 3.9+ with numpy installed (pip install numpy); the core algorithm needs only the standard library, and NumPy is used for the batch cross-check.
  • Familiarity with the I-MR chart and the span-2 moving range constant $d_2 = 1.128$ — see individual moving range (I-MR) charts.
  • Data arriving one observation at a time (or in small batches) from a live source, already cleaned on ingest so gaps are explicit rather than silent.
  • A decision on the Phase I baseline window: streaming updates are for Phase II monitoring against frozen limits, or for a governed baseline refresh, not for silently re-anchoring limits every point.
  • Understanding that I-MR is the n = 1 case; if you can form rational subgroups, prefer the subgrouped engine instead.

Step 1 — Track a running mean and variance with Welford's algorithm

The naive route accumulates a running sum and sum-of-squares and derives variance as $E[x^2] - E[x]^2$. On a long stream those two large numbers grow until their difference loses most of its significant digits — catastrophic cancellation — and the variance can even come out negative. Welford's method avoids this by updating the mean and the sum of squared deviations $M_2$ directly, so every quantity stays on the scale of the data.

from __future__ import annotations

from dataclasses import dataclass, field


@dataclass
class RunningStats:
    """Welford running mean and variance, updated one observation at a time."""

    count: int = 0
    mean: float = 0.0
    m2: float = 0.0  # running sum of squared deviations from the mean

    def update(self, x: float) -> None:
        """Fold a single observation into the running statistics."""
        if x != x:  # NaN guard: never let a gap poison the state
            raise ValueError("Received NaN; resolve missing values before updating.")
        self.count += 1
        delta = x - self.mean
        self.mean += delta / self.count
        delta2 = x - self.mean
        self.m2 += delta * delta2

    @property
    def variance(self) -> float:
        """Sample variance (n - 1 denominator); NaN until two points seen."""
        if self.count < 2:
            return float("nan")
        return self.m2 / (self.count - 1)

    @property
    def std(self) -> float:
        return self.variance ** 0.5

Each update is a handful of floating-point operations regardless of how long the stream has run, and because updates fold in one at a time the mean tracks the true running average exactly.

Step 2 — Maintain an incremental moving range

The I-MR chart estimates variation from the span-2 moving range, $MR_i = |x_i - x_{i-1}|$. Averaging those does not need history either — keep the previous value and a running mean of the moving ranges, again with Welford so the moving-range mean is stable.

@dataclass
class IncrementalMovingRange:
    """Running mean of the span-2 moving range |x_i - x_(i-1)|."""

    _prev: float | None = None
    _mr_stats: RunningStats = field(default_factory=RunningStats)

    def update(self, x: float) -> float | None:
        """Fold in one observation; return the new moving range, or None first."""
        if self._prev is None:
            self._prev = x
            return None
        mr = abs(x - self._prev)
        self._prev = x
        self._mr_stats.update(mr)
        return mr

    @property
    def mr_bar(self) -> float:
        """Mean moving range; NaN until at least one MR has been seen."""
        if self._mr_stats.count < 1:
            return float("nan")
        return self._mr_stats.mean

The first observation seeds _prev and produces no moving range; from the second point on, each arrival yields exactly one moving range and one constant-time update.

Step 3 — Assemble live I-MR limits

I-MR limits use the mean moving range and the span-2 constant $d_2 = 1.128$. The estimated process sigma is $\hat{\sigma} = \overline{MR}/d_2$, giving the Individuals limits $\bar{X} \pm 3\hat{\sigma}$, which reduce to $\bar{X} \pm 2.66\,\overline{MR}$. The moving-range chart uses $\text{UCL} = D_4\overline{MR} = 3.267\,\overline{MR}$ with a lower limit of zero at span 2.

_D2_SPAN2 = 1.128
_D4_SPAN2 = 3.267  # D4 for the span-2 moving range (n = 2)


@dataclass
class StreamingIMR:
    """Live I-MR control limits maintained incrementally on a stream."""

    _x: RunningStats = field(default_factory=RunningStats)
    _mr: IncrementalMovingRange = field(default_factory=IncrementalMovingRange)

    def update(self, x: float) -> None:
        """Fold one observation into both the Individuals and MR statistics."""
        self._mr.update(x)   # uses previous value, so update MR before overwriting
        self._x.update(x)

    def limits(self) -> dict[str, float]:
        """Return current I-MR limits from the running statistics.

        Raises
        ------
        ValueError
            If fewer than two observations have been seen, so no moving
            range and hence no dispersion estimate exists yet.
        """
        if self._x.count < 2:
            raise ValueError("Need at least two observations for I-MR limits.")
        x_bar = self._x.mean
        mr_bar = self._mr.mr_bar
        sigma_hat = mr_bar / _D2_SPAN2
        return {
            "count": self._x.count,
            "x_bar": x_bar,
            "mr_bar": mr_bar,
            "i_ucl": x_bar + 3.0 * sigma_hat,   # X-bar + 2.66 * MR-bar
            "i_lcl": x_bar - 3.0 * sigma_hat,   # X-bar - 2.66 * MR-bar
            "mr_ucl": _D4_SPAN2 * mr_bar,
            "mr_lcl": 0.0,                      # span-2 MR lower limit is zero
        }

Update the moving range before the Individuals statistic on each point, because the moving range needs the previous value; overwriting _prev first would compute $|x_i - x_i| = 0$ and quietly deflate every limit.

Step 4 — Drive the stream and freeze the baseline

In Phase II the limits come from a frozen Phase I baseline; the incremental state is what lets you monitor cheaply against them and, when a governed change calls for it, propose a refreshed baseline without a full recompute.

def run_stream(baseline: list[float], live: list[float]) -> dict[str, float]:
    """Freeze limits from a baseline window, then monitor a live stream.

    The baseline window establishes the frozen Phase I limits; live points are
    checked against them. Promoting a new baseline is a governed action, not an
    automatic side effect of monitoring.
    """
    engine = StreamingIMR()
    for value in baseline:
        engine.update(value)
    frozen = engine.limits()  # Phase I: freeze and stop moving the limits

    for value in live:
        # Phase II: evaluate against frozen limits without recomputing them.
        if value > frozen["i_ucl"] or value < frozen["i_lcl"]:
            # Signal handling / alerting hooks in here; limits stay frozen.
            pass
    return frozen

Freezing means capturing engine.limits() once and comparing subsequent points against that snapshot, exactly as the batch engine emits an immutable frame. Refreshing the baseline under tool wear is the job of rolling window limit recalibration, triggered by an engineering change, not by the monitoring loop deciding on its own.

Verification

Prove the incremental limits equal a from-scratch batch recompute over the same points, and prove Welford beats the naive sum-of-squares on a hard case.

import numpy as np


def batch_imr_limits(x: list[float]) -> dict[str, float]:
    """Reference I-MR limits computed the slow, from-scratch way."""
    arr = np.asarray(x, dtype="float64")
    mr = np.abs(np.diff(arr))
    x_bar, mr_bar = arr.mean(), mr.mean()
    sigma = mr_bar / 1.128
    return {"i_ucl": x_bar + 3 * sigma, "i_lcl": x_bar - 3 * sigma,
            "mr_ucl": 3.267 * mr_bar}


rng = np.random.default_rng(11)
sample = (50_000 + rng.normal(0, 2, size=5_000)).tolist()  # large offset stresses variance

engine = StreamingIMR()
for v in sample:
    engine.update(v)
stream_limits = engine.limits()
ref = batch_imr_limits(sample)

for key in ("i_ucl", "i_lcl", "mr_ucl"):
    assert np.isclose(stream_limits[key], ref[key], atol=1e-6), key

# Welford variance stays correct where naive E[x^2] - E[x]^2 loses precision:
naive_var = (np.mean(np.square(sample)) - np.mean(sample) ** 2)
assert np.isclose(engine._x.variance, np.var(sample, ddof=1), atol=1e-6)
print("streaming I-MR matches batch recompute:", stream_limits)

The incremental and batch limits agree to floating-point tolerance because they compute the same statistics in a different order. The large 50,000 offset is deliberate: it is exactly where the naive sum-of-squares formula sheds significant digits while Welford's running variance stays accurate. Once the limits are stable, the same live points flow into out-of-control rule detection for signal evaluation.

Root-Cause Table

Symptom Cause Fix
Streaming variance drifts from the batch value Naive running sum-of-squares suffered catastrophic cancellation Use Welford's mean and M2 update; never compute variance as E[x^2] minus E[x] squared on a long stream
Every limit is too tight and MR-bar is near zero Individuals value overwrote the previous value before the moving range read it Update the moving range first, then the Individuals statistic, each point
Limits keep moving on every new point The monitoring loop recomputes the baseline instead of freezing it Freeze Phase I limits once; refresh only through governed recalibration
First call to limits() raises Fewer than two observations, so no moving range exists Wait for at least two points; the span-2 moving range is undefined for a single value
Variance briefly reads NaN Only one observation folded in so far Expected: sample variance needs two points; guard callers accordingly
A gap silently deflates the moving range A missing reading was folded in as zero or skipped inconsistently Resolve missing values on ingest; the NaN guard rejects them rather than poisoning the state

This how-to belongs to vectorized control limit calculation engines, within automated control chart generation and calculation.