Control charts (Shewhart, CUSUM, EWMA)

Deep dive

How it works

Parts:

• Phase I sample — m observations from a period without known disturbance; gives centre mu and spread sigma (individuals: mean moving range / 1.128; subgroups: R-bar/d2 or s-bar/c4).
• Shewhart rule — signal when one standardized value has |z| > 3.
• CUSUM — two running sums S+, S- with allowance k and decision limit h.
• EWMA — smoothed statistic w with weight lambda and limit factor L.
• Count charts — c, u (Poisson counts), p (binomial fractions).
• Run rules (optional) — extra patterns such as 8 points on one side of the centre line.

Operations:

1. Phase I: chart a stable period, drop points with a known cause, repeat until none signals. mu = mean; sigma = mean(|x_i - x_(i-1)|) / 1.128, where 1.128 is d2 for n = 2 (NIST 6.3.2.2). Subgroups: plot the subgroup mean against mu +- 3 sigma/sqrt(n).
2. Phase II, per observation: z = (x - mu) / sigma.
3. Shewhart: signal if |z| > 3.
4. CUSUM: S+ = max(0, S+ + z - k), S- = max(0, S- - z - k); signal if either exceeds h.
5. EWMA: w = lambda z + (1 - lambda) w from w = 0; signal if |w| > L sqrt(lambda / (2 - lambda)).
6. On a signal: fix the cause, reset sums and w; refit Phase I only after a deliberate process change.
7. Design: fix ARL0 from the false-alarm budget, then check ARL1 at the smallest shift worth catching.

 Phase I (once)                           Phase II (each observation)

 x_1..x_m --> drop known causes           x_t --> z = (x_t - mu) / sigma
          --> mu, sigma = mean(MR)/1.128 ---+--> Shewhart  |z| > 3 ----------+
                                            +--> CUSUM     S+ or S- > h -----+--> signal
                                            +--> EWMA      |w| > limit ------+      |
     ^                                                                              v
     +---- refit only after a deliberate change <------------------------------ find cause, fix

ARL. Each point is a hypothesis test; with signal probability p per independent point the run length is geometric, P(RL = t) = (1-p)^(t-1) p, mean 1/p. Normal data, 3 sigma: p = 2(1 - Phi(3)) = 0.0027, ARL0 = 370.4. After a +1 sigma shift the upper limit is 2 sigma away: p = 1 - Phi(2) = 0.02275, ARL1 = 43.96 (NIST 6.3.2.3.1, ignoring the 0.000032 lower tail). The geometric median ln 0.5 / ln(1-p) is 256: half of all false alarms come before observation 256.

Counts (NIST 6.3.3). c chart: Poisson defects per unit, mean and variance c, limits c +- 3 sqrt(c). u chart, unit size n varying: plot count/n, variance u-bar/n (derived), limits u-bar +- 3 sqrt(u-bar/n). p chart: p +- 3 sqrt(p(1-p)/n). Skew costs false alarms: at c = 0.8 the upper limit is 3.48 and P(count > 3) = 0.009, not 0.00135 (NIST 6.3.1).

Run rules. Western Electric rules add 2 of 3 beyond 2 sigma, 4 of 5 beyond 1 sigma, 8 in a row on one side, trends: more sensitive to drifts, but false alarms come every 91.75 points instead of 371 (Champ and Woodall 1987, via NIST 6.3.2).

CUSUM design (Page 1954). The log-likelihood-ratio increment of N(delta, 1) against N(0, 1) is delta (z - delta/2): CUSUM is the SPRT recursion restarted at zero, k = delta/2. Rule of thumb k = 0.5, h = 4 to 5 (Montgomery, via NIST 6.3.2.3); beats Shewhart for shifts of 2 sigma or less (NIST).

EWMA design (Roberts 1959). The variance of w tends to lambda/(2-lambda) times that of z; lambda = 1 is Shewhart; usual lambda 0.2 to 0.3 (Hunter 1986, via NIST 6.3.2.4); Lucas and Saccucci (1990) tabulate L. It is exponential smoothing of z.

Autocorrelation. All limits assume independent points (NIST 6.3.2.4). For AR(1) with coefficient phi, Var(x_t - x_(t-1)) = 2 sigma^2 (1 - phi), so the moving range estimates sigma sqrt(1 - phi): 0.707 sigma at phi = 0.5. Limits shrink to 2.12 sigma, p = 2(1 - Phi(2.12)) = 0.0339, ARL0 near 29.5 for independent points; measured below: 33.8 (alarms cluster). Remedy: chart residuals of a time-series model (ARIMA).

Capability (stable process only; NIST 6.1.6): Cp = (USL - LSL) / (6 sigma), Cpk = min(USL - mu, mu - LSL) / (3 sigma); needs normal data, about 50 independent values.

Code

Std-only, seeded. Moving-range Phase I, then Shewhart, CUSUM and EWMA; ARL0 and ARL1 (1 sigma shift) from 20 000 runs each (standard error about 0.7 %). Asserts: Shewhart within 5 % of 370.4 and 43.96, CUSUM and EWMA at least 3x faster at the shift, autocorrelated ARL0 below 50. h = 4.8 and L = 2.7 came from a grid search with this simulation until ARL0 was near 370.

struct Rng(u64);
impl Rng {
    fn next(&mut self) -> f64 { // xorshift64*, uniform in (0,1)
        self.0 ^= self.0 >> 12; self.0 ^= self.0 << 25; self.0 ^= self.0 >> 27;
        ((self.0.wrapping_mul(0x2545F4914F6CDD1D) >> 11) as f64 + 0.5) / (1u64 << 53) as f64
    }
    fn normal(&mut self) -> f64 { // Box-Muller
        (-2.0 * self.next().ln()).sqrt() * (std::f64::consts::TAU * self.next()).cos()
    }
}

/// Phase I: centre = mean, sigma = mean moving range / d2 (d2 = 1.128 for n = 2).
fn phase1(x: &[f64]) -> (f64, f64) {
    let mu = x.iter().sum::<f64>() / x.len() as f64;
    let mr = x.windows(2).map(|w| (w[1] - w[0]).abs()).sum::<f64>() / (x.len() - 1) as f64;
    (mu, mr / 1.128)
}

#[derive(Clone, Copy)]
enum Chart { Shewhart, Cusum { k: f64, h: f64 }, Ewma { lam: f64, l: f64 } }

/// Phase II stream: mean `shift`, unit sd, AR(1) coefficient `phi`; limits use (mu, sigma).
#[derive(Clone, Copy)]
struct Stream { shift: f64, phi: f64, mu: f64, sigma: f64 }

fn run_length(c: Chart, s: Stream, rng: &mut Rng) -> u64 {
    let (mut e, mut hi, mut lo, mut w) = (rng.normal(), 0.0f64, 0.0f64, 0.0f64);
    for t in 1..10_000_000u64 {
        e = s.phi * e + (1.0 - s.phi * s.phi).sqrt() * rng.normal(); // unit marginal variance
        let z = (s.shift + e - s.mu) / s.sigma;
        let alarm = match c {
            Chart::Shewhart => z.abs() > 3.0,
            Chart::Cusum { k, h } => {
                hi = (hi + z - k).max(0.0); // one-sided sums, reset at 0
                lo = (lo - z - k).max(0.0);
                hi > h || lo > h
            }
            Chart::Ewma { lam, l } => {
                w = lam * z + (1.0 - lam) * w;
                w.abs() > l * (lam / (2.0 - lam)).sqrt() // asymptotic limit
            }
        };
        if alarm { return t; }
    }
    u64::MAX
}

fn arl(c: Chart, s: Stream, n: u64, seed: u64) -> f64 {
    let mut rng = Rng(seed);
    (0..n).map(|_| run_length(c, s, &mut rng) as f64).sum::<f64>() / n as f64
}

fn sample(n: usize, phi: f64, rng: &mut Rng) -> Vec<f64> {
    let mut e = rng.normal();
    (0..n).map(|_| { e = phi * e + (1.0 - phi * phi).sqrt() * rng.normal(); e }).collect()
}

fn main() {
    let (n, mut rng) = (20_000u64, Rng(7)); // runs per ARL estimate
    let ideal = Stream { shift: 0.0, phi: 0.0, mu: 0.0, sigma: 1.0 };
    let one = Stream { shift: 1.0, ..ideal };
    let charts = [("Shewhart", Chart::Shewhart), ("CUSUM", Chart::Cusum { k: 0.5, h: 4.8 }),
                  ("EWMA", Chart::Ewma { lam: 0.1, l: 2.7 })];
    let mut a1 = vec![];
    for (name, c) in charts {
        let t = std::time::Instant::now();
        let (a0, a) = (arl(c, ideal, n, 1), arl(c, one, n, 2));
        let ns = t.elapsed().as_nanos() as f64 / ((a0 + a) * n as f64);
        println!("{name:9} ARL0 {a0:6.1}  ARL1 {a:5.2}  {ns:.1} ns/obs");
        a1.push((a0, a));
    }
    assert!((a1[0].0 / 370.4 - 1.0).abs() < 0.05, "Shewhart ARL0 vs 1/0.0027");
    assert!((a1[0].1 / 43.96 - 1.0).abs() < 0.05, "Shewhart ARL1 vs 1/0.02275");
    for i in 1..3 { assert!(a1[i].1 < a1[0].1 / 3.0 && a1[i].0 > 300.0); }

    // Phase I: sigma = mean moving range / 1.128 on an AR(1) process (phi = 0.5).
    let (mu, s_ar) = phase1(&sample(100_000, 0.5, &mut rng));
    let a_ar = arl(Chart::Shewhart, Stream { phi: 0.5, mu, sigma: s_ar, ..ideal }, n, 3);
    println!("AR(1) phi=.5: sigma-hat {s_ar:.3}, Shewhart ARL0 {a_ar:.1}");
    assert!((s_ar - 0.5f64.sqrt()).abs() < 0.02 && a_ar < 50.0);
}

Output (fixed seeds; ns/obs varies):

Shewhart  ARL0  368.9  ARL1 43.74  15.3 ns/obs
CUSUM     ARL0  382.1  ARL1 10.02  15.5 ns/obs
EWMA      ARL0  365.1  ARL1  9.77  15.6 ns/obs
AR(1) phi=.5: sigma-hat 0.705, Shewhart ARL0 33.8

Complexity and constants

• O(1) time per observation; state 0 floats (Shewhart), 2 (CUSUM), 1 (EWMA). Phase I: one O(m) pass.
• Measured on this Mac (Apple M4 Pro), release build, one thread: 15 to 16 ns per observation including the random generator (17 to 46 ns while the machine was busy).

Efficiency

• Cost: about 15 ns per observation, measured on this Mac (Apple M4 Pro), random generator included; the chart update itself is a few flops. 1e6 series updated once a minute: 1e6 x 15 ns = 15 ms per minute (arithmetic). Bound: input rate.
• Levers:
  • CUSUM or EWMA beside Shewhart — ARL1 10 vs 44 at equal ARL0 (measured here) — shifts of 2 sigma or less.
  • Subgroups of n — a 1 sigma shift becomes sqrt(n) sigma of the mean; n = 4: Shewhart ARL 6.30 subgroups = 25 observations (NIST 6.3.2.3.1) — sampling cheap, delay costly.
  • Count per interval, not per event (c, u) — fewer points to chart — event streams.
• Naive → standard → state of the art:
  • Standard deviation of all history, 3 sigma per point · limits widen with past disturbances
  • Moving-range sigma, Shewhart plus tuned CUSUM or EWMA · about 15 ns, 1 to 3 floats
  • Limits designed from exact ARL (R package spc) · same runtime, design offline

Implementation notes

• Version mu, sigma and parameters per series; reset CUSUM and EWMA after a handled alarm.
• Subgroups: within-subgroup spread, not the spread of all points. A robust sigma (robust statistics) needs its own constant, not 1.128.
• Low counts: Poisson or binomial quantile limits (qcc, NEWS 2.7).

Pitfalls

• Phase I data that already contains the shift: the limits absorb it.
• Reading ARL0 = 370 as "one false alarm per 370 points": half come before point 256.
• Stacking run rules without recomputing ARL0 (371 to 91.75, NIST).

Variants and alternatives

In production

Factories chart dimensions and defect counts (NIST examples: wafers, flow rates); hospitals plot infection and complication rates (Mohammed et al. 2008); clinical chemistry laboratories use a multi-rule Shewhart chart (Westgard et al. 1981).

Implementations

Sources

• NIST/SEMATECH e-Handbook of Statistical Methods, 6.3 Univariate and Multivariate Control Charts and 6.1.6 Process Capability. 6.3 · 6.1.6
• E. S. Page, Continuous inspection schemes, Biometrika 41(1-2), 1954. paper
• S. W. Roberts, Control chart tests based on geometric moving averages, Technometrics 1(3), 1959. paper
• James M. Lucas, Michael S. Saccucci, Exponentially weighted moving average control schemes: properties and enhancements, Technometrics 32(1), 1990. paper
• Charles W. Champ, William H. Woodall, Exact results for Shewhart control charts with supplementary runs rules, Technometrics 29(4), 1987. paper
• J. O. Westgard, P. L. Barry, M. R. Hunt, T. Groth, A multi-rule Shewhart chart for quality control in clinical chemistry, Clinical Chemistry 27(3), 1981. paper
• M. A. Mohammed, P. Worthington, W. H. Woodall, Plotting basic control charts: tutorial notes for healthcare practitioners, Quality and Safety in Health Care 17(2), 2008. paper