Advanced Statistical Process Monitoring: CUSUM, EWMA & Multivariate SPC

Standard Shewhart control charts are "memoryless" filters: the decision regarding whether a process is in control at inspection time $t$ depends exclusively on the single sample statistic $W_t$ observed at that instantaneous moment. All prior historical data points $W_{t-1}, W_{t-2}, \dots, W_1$ are discarded.

While memorylessness makes Shewhart charts extraordinarily robust against large, abrupt catastrophic shocks (shifts of magnitude $\delta \ge 2.5\sigma$, where detection occurs on the very first sample with near $100\%$ probability), it severely degrades their statistical power against small, persistent process shifts ($0.5\sigma \le \delta \le 1.5\sigma$).

As derived in Unit 1, when a process experiences a $0.5\sigma$ mean shift with subgroup size $n = 4$:

• The probability of detection on a single sample is only $p = 1 - \beta = 0.0228$.
• The Average Run Length is $ARL_1 = 1 / 0.0228 \approx 44$ subgroups.
• If subgroups are drawn hourly, forty-four hours of production—nearly two full days—will elapse before the shift is detected, allowing tens of thousands of off-target parts to enter the supply chain.

To eliminate this operational vulnerability, statistical quality engineering incorporates cumulative memory into the charting statistic. The two preeminent memory-based control charts are the Cumulative Sum (CUSUM) chart and the Exponentially Weighted Moving Average (EWMA) chart.

Originally introduced by E. S. Page in 1954, the CUSUM chart directly incorporates all past process history by plotting cumulative sums of deviations of sample values from a target value $\mu_0$.

Let $\bar{X}_i$ denote the sample mean of the $i$-th rational subgroup of size $n$, with in-control mean $\mu_0$ and known standard deviation $\sigma$. The standardized subgroup mean is:

The Two-Sided Tabular (Algorithmic) CUSUM

While early practitioners used graphical V-masks, modern industrial practice universally employs the tabular CUSUM. The tabular CUSUM maintains two separate one-sided accumulators:

• $C_i^+$: The Upper CUSUM, which accumulates positive deviations above target (detecting upward mean shifts).
• $C_i^-$: The Lower CUSUM, which accumulates negative deviations below target (detecting downward mean shifts).

Starting with initial conditions $C_0^+ = 0$ and $C_0^- = 0$, the recursive update equations are:

where:

• $K$: The reference value (slack value), conventionally chosen halfway between the in-control mean $\mu_0$ and the out-of-control target $\mu_1 = \mu_0 + \delta \sigma / \sqrt{n}$:

In standardized units, $k = K / (\sigma / \sqrt{n}) = \delta / 2$.

• $H$: The decision interval (threshold), above which an out-of-control alarm is triggered:

In standardized units, $h = H / (\sigma / \sqrt{n})$. Typically, quality engineers set $k = 0.5$ and $h = 4.0$ or $5.0$.

Whenever the accumulation terms become negative, the $\max(0, \cdot)$ operator resets the accumulator to zero. This resetting mechanism prevents the accumulation of "favorable" common-cause noise, ensuring that the CUSUM reacts rapidly as soon as an adverse shift commences.

Counter Statistics for Shift Diagnosis

The tabular CUSUM also maintains run counters $N_i^+$ and $N_i^-$, recording the number of consecutive periods since the respective accumulator was last zero:

When an alarm occurs at sample $t$, the estimated shift onset time is $t - N_t^+ + 1$, and the estimated magnitude of the shifted mean is:

Historically, CUSUM charts were plotted as raw cumulative sums $S_m = \sum_{i=1}^m (\bar{X}_i - \mu_0)$ against a V-shaped cardboard mask placed over the most recent data point. The V-mask is defined by two geometric parameters:

• Lead distance $d$: The distance between the current vertex and the plotted point.
• Angle $\theta$: The half-angle of the V-mask opening.

George A. Barnard (1959) and E. S. Page proved that the V-mask CUSUM is mathematically identical to Wald’s Sequential Probability Ratio Test (SPRT) for testing the composite hypothesis:

The decision boundaries of the SPRT translate directly into the linear slopes of the V-mask arms:

where $A$ is the horizontal-to-vertical plotting scale factor. If any historical point falls behind the arms of the V-mask, the null hypothesis $H_0$ is rejected, triggering an out-of-control alarm. Because the tabular CUSUM executes identical mathematics algorithmically without drafting errors, it has superseded physical V-masks.

Evaluating the performance of a CUSUM chart requires computing its Average Run Length ($ARL$). Because the CUSUM accumulator is continuous and bounded below by zero, run lengths do not follow a simple geometric distribution.

Siegmund's Continuous Brownian Approximation

David Siegmund (1985) derived an accurate analytical continuous approximation for the ARL of a one-sided CUSUM using Brownian motion boundary-crossing theory:

where:

The constant $1.166$ represents an overshoot correction factor accounting for the discrete sampling steps across the continuous decision threshold $h$.

When $\Delta = 0$ (the shift magnitude exactly matches the reference value $k$), applying L'Hôpital's rule yields:

For a two-sided CUSUM, the overall average run length $ARL_{\text{two-sided}}$ is obtained from the two individual one-sided run lengths via the parallel hazard formula:

Brook and Evans Markov Chain Method

D. Brook and D. A. Evans (1972) formulated the exact numerical evaluation of CUSUM ARL by discretizing the continuous interval $[0, h]$ into $t$ discrete transient states of width $w = h / (t - 0.5)$, with state $t + 1$ representing the absorbing out-of-control state.

The transition probability matrix $\mathbf{P}$ of the Markov chain has the canonical form:

where $\mathbf{R}$ is the $t \times t$ submatrix of transition probabilities between transient states. The fundamental matrix of the Markov chain is:

The vector of expected run lengths starting from each initial state is:

The standard zero-start $ARL$ corresponds to the first element $\mu_{\text{ARL}, 1}$, representing the expected steps to absorption starting from $C_0 = 0$.

S. W. Roberts (1959) introduced the EWMA control chart, which monitors an exponentially smoothed moving average that assigns geometrically decaying weights to historical observations.

Let $X_t$ denote the $t$-th observation (or subgroup mean $\bar{X}_t$) at time $t = 1, 2, 3, \dots$. The EWMA statistic $Z_t$ is defined recursively:

where:

• $\lambda \in (0, 1]$: The smoothing constant (weight assigned to the newest observation).
• $Z_0$: The starting value, initialized to the in-control process target: $Z_0 = \mu_0$.

Expanding the recursion reveals the explicit geometric weighting:

Notice two profound mathematical properties:

1. If $\lambda = 1$, $Z_t = X_t$, and the EWMA collapses into the memoryless Shewhart control chart.
2. As $\lambda \to 0$, the EWMA assigns virtually equal weight to all past history, behaving like a cumulative average.

Variance and Control Limits

Assuming independent observations $X_t$ with common variance $\sigma_X^2 = \sigma^2 / n$:

The exact, time-varying $3\sigma$ control limits for the EWMA chart are:

where $L$ is the control limit width multiplier (typically $L \approx 2.7$ to $3.0$).

As time progresses ($t \to \infty$), the term $(1 - \lambda)^{2t} \to 0$, and the limits converge to steady-state horizontal asymptotic limits:

Industrial practice typically uses exact time-varying limits for the first $t \le 10$ samples to prevent false alarms during startup, switching to steady-state limits thereafter.

The design parameters of an EWMA chart are the smoothing constant $\lambda$ and limit multiplier $L$. J. M. Lucas and M. S. Saccucci (1990) published extensive design tables:

• To optimize detection of small shifts ($\delta \approx 0.5\sigma$ to $1.0\sigma$), select small smoothing values: $\lambda \in [0.05, 0.20]$. A universal industry standard is $\lambda = 0.10$ and $L = 2.70$, which delivers $ARL_0 \approx 370$ and $ARL_1(\delta = 1.0) \approx 3.8$ (compared to Shewhart’s $ARL_1 \approx 44$).
• To optimize detection of larger shifts ($\delta \ge 2.0\sigma$), select $\lambda \in [0.20, 0.40]$.

Fast Initial Response (FIR)

If a process is already out of control when monitoring commences (for example, after maintenance shutdown or raw material reload), standard charts starting at $C_0 = 0$ or $Z_0 = \mu_0$ require multiple periods to accumulate enough evidence to signal.

To resolve this, Lucas proposed the Fast Initial Response (FIR) mechanism:

• In CUSUM: Initialize accumulators at half the decision interval: $C_0^+ = H / 2$ and $C_0^- = H / 2$ (a "head-start"). If the process is in control, the accumulator rapidly decays to zero; if the process is out of control, it triggers an alarm within 1 to 2 samples.
• In EWMA: Construct narrow, expanding parabolic limits during initial startup.

Modern industrial processes monitor dozens or hundreds of correlated process variables simultaneously (e.g., temperatures, pressures, flow rates, and vibration spectra in a chemical cracking reactor). Monitoring each variable with an individual Shewhart chart commits two severe statistical errors:

1. Distortion of False Alarm Probability: If $p$ independent quality characteristics are monitored simultaneously at $\alpha = 0.0027$, the overall joint false alarm probability across the system inflates to:

For $p = 10$ variables, $\alpha_{\text{joint}} \approx 2.67\%$, dropping the in-control ARL from $370$ down to $37$!

2. Failure to Detect Covariance Shifts: Individual univariate charts monitor rectangular control regions:

When variables are correlated, the true in-control distribution occupies an inclined hyper-ellipsoid. A point may fall comfortably within the univariate limits of every individual variable, yet represent a severe multivariate outlier that violates the physical correlation structure.

Hotelling's $T^2$ Statistic

Harold Hotelling (1947) formulated the multivariate generalization of Student’s $t$-statistic. Let $\mathbf{X}_i = (X_{i1}, \dots, X_{ip})^T$ denote a $p \times 1$ vector of observations from a multivariate normal distribution $\mathcal{N}_p(\boldsymbol{\mu}, \boldsymbol{\Sigma})$.

When baseline parameters are estimated from $m$ historical subgroups of size $n$, with sample mean vector $\bar{\bar{\mathbf{X}}}$ and pooled sample covariance matrix $\mathbf{S}$:

For a new subgroup mean vector $\bar{\mathbf{X}}$, the Hotelling's $T^2$ statistic is:

Control Limits for Hotelling's $T^2$

The $T^2$ statistic is a non-negative scalar representing the squared Mahalanobis distance from the process center. The lower control limit is naturally zero ($LCL = 0$).

The Upper Control Limit depends on the operational phase:

• Phase I (Retrospective Analysis):

• Phase II (Prospective Monitoring):

For individual observations ($n = 1$), the covariance matrix is estimated from successive differences, and Phase II limits follow a Beta or scaled $F$-distribution.

Out-of-Control Signal Decomposition (MYT Method)

When Hotelling’s $T^2$ breaches $UCL$, the scalar statistic indicates that an assignable cause has occurred, but provides no insight into which physical variable or covariance relationship caused the alarm.

Mason, Young, and Tracy (MYT, 1995) formulated an exact orthogonal decomposition of $T^2$ into unconditional and conditional components:

where $T_j^2$ is the univariate squared $t$-statistic for variable $j$, and $T_{-j \mid j}^2$ is the conditional Hotelling statistic of all remaining variables conditioned upon the observed value of variable $j$. By systematically computing these orthogonal sub-components, engineers isolate the exact sensor responsible for the disturbance.

Just as univariate Shewhart charts miss small shifts, Hotelling’s $T^2$ is memoryless and insensitive to small multivariate shifts. C. A. Lowry et al. (1992) developed the Multivariate EWMA (MEWMA).

The MEWMA accumulates a $p \times 1$ vector $\mathbf{Z}_t$:

where $\boldsymbol{\Lambda} = \text{diag}(\lambda_1, \dots, \lambda_p)$, usually with a common smoothing parameter $\lambda_1 = \dots = \lambda_p = \lambda$. The MEWMA plotting statistic is:

High-Dimensional Sensor Monitoring via PCA

When processes feature hundreds of collinear sensors ($p \gg n$), the covariance matrix $\mathbf{S}$ becomes ill-conditioned or singular. Modern industrial informatics applies Principal Component Analysis (PCA):

1. Decompose the standardized data matrix into $k$ dominant principal components ($k \ll p$) capturing $85-95\%$ of total variance:

1. Monitor systematic process shifts in the model subspace using Hotelling's $T^2$ on scores:

1. Monitor abnormal breakdowns in sensor correlations using the Squared Prediction Error ($SPE$ / $Q$-statistic) in the residual subspace:

This dual-statistic framework ($T^2$ and $Q$) forms the core of modern Industry 4.0 predictive maintenance and process anomaly detection.

A foundational assumption underlying both Shewhart and CUSUM/EWMA control charts is that successive observations $X_1, X_2, \dots, X_t$ are statistically independent. In modern automated manufacturing—such as continuous chemical distillation, petroleum refining, glass sheet floating, and paper mill web rollers—measurements are logged by high-frequency digital sensors every few seconds. Because the physical volume of fluid or thermal mass in the vessel cannot change instantaneously, consecutive measurements exhibit strong positive autocorrelation ($\rho_1 = \text{Corr}(X_t, X_{t-1}) > 0$).

Impact of Autocorrelation on Classical SPC

When positive autocorrelation is present, applying standard Shewhart limits computed from sample variance generates catastrophic false alarm rates:

1. The sample variance $S^2$ severely underestimates the true long-term process variance $\sigma_X^2$:

1. The true variance of the sample mean is inflated:

1. As a result, the control limits are set artificially too narrow, and the process wanders in smooth serial waves. For an autoregressive process $AR(1)$ with $\rho_1 = 0.70$, the actual in-control ARL collapses from $370$ down to below 15! Plant personnel are overwhelmed with continuous false alarms and eventually disable the monitoring software.

The Residual Control Chart (ARIMA Modeling)

To restore statistical integrity to autocorrelated process monitoring, George Box and T. Luceño (1997) introduced the Residual Control Chart (also called the Special Cause Chart, SCC):

1. Fit a time series model (typically an autoregressive integrated moving average model, $\text{ARIMA}(p, d, q)$) to the in-control baseline data:

1. Compute the one-step-ahead forecasts $\widehat{X}_{t \mid t-1}$ in real time.
2. Compute the one-step-ahead forecast errors (residuals):

1. Because the ARIMA filter captures all deterministic serial correlation, the residuals $e_t$ are independent and identically distributed Gaussian white noise: $e_t \stackrel{\text{i.i.d.}}{\sim} \mathcal{N}(0, \sigma_e^2)$.
2. Plot the residuals $e_t$ on a standard Shewhart, CUSUM, or EWMA control chart with limits:

An assignable cause in the original physical process manifests as a sudden bias or structural breakdown in the residuals, triggering a statistically valid, uncorrupted alarm.

When industrial process data violate normality assumptions and cannot be normalized via mathematical transformations (for instance, multimodal distributions resulting from mixed supply batches, or Cauchy-like heavy-tailed distributions in electronic noise), standard parametric control limits fail. In such environments, quality statisticians deploy Nonparametric (Distribution-Free) Control Charts.

The Shewhart Sign Control Chart

To monitor the median $\widetilde{\mu}$ of an arbitrary continuous distribution without distributional assumptions:

1. In subgroup $i$ of size $n$, count the number of observations strictly exceeding the target median $\widetilde{\mu}_0$:

1. Under the in-control null hypothesis, each observation has probability $P(X_{ij} > \widetilde{\mu}_0) = 0.5$, so $S_i \sim \text{Binomial}(n, 0.5)$.
2. Because the distribution of $S_i$ is completely known and independent of the parent continuous distribution $F(x)$, exact symmetric binomial control limits are established.

The Wilcoxon Signed-Rank and Rank-Sum Charts

To increase statistical power, nonparametric charts incorporate rank magnitude:

• For paired or symmetric continuous data, the Wilcoxon Signed-Rank Statistic $W_i$ sums the ranks of absolute deviations $|X_{ij} - \widetilde{\mu}_0|$ with sign indicators.
• For two-sample Phase II monitoring against a Phase I reference sample, the Mann-Whitney-Wilcoxon Rank-Sum Statistic combines reference and prospective samples, counting rank inversions.
• Nonparametric CUSUM and EWMA charts replace raw values with standardized rank scores, achieving stable, invariant in-control $ARL_0$ values across any continuous distribution family.

In modern automated manufacturing and smart machining systems, the quality of a process is increasingly characterized not by a scalar measurement or a fixed vector, but by a functional relationship—a profile—relating a dependent response variable $Y$ to one or more independent explanatory variables $X$:

Examples of industrial profiles include:

• The force-displacement signature curve of an automated metal stamping press.
• The torque-angle rotation profile of a robotic nut-runner tightening automotive cylinder head bolts.
• The spectrophotometric absorbance calibration curve of an in-line optical sensor.

Linear Profile Monitoring

In the baseline Phase I formulation, suppose the profile can be modeled as a simple linear function:

where the explanatory design points $X_1, X_2, \dots, X_n$ remain fixed across all inspection cycles. For each sample $i$, ordinary least squares (OLS) regression yields sample estimates of the intercept $\widehat{a}_{0i}$, slope $\widehat{a}_{1i}$, and residual mean square error $MSE_i$.

Let $\mathbf{B}_i = (\widehat{a}_{0i}, \widehat{a}_{1i})^T$. Because the intercept and slope estimators are correlated ($\text{Cov}(\widehat{a}_0, \widehat{a}_1) = -\bar{X} \sigma^2 / S_{xx}$), univariate monitoring of slope or intercept independently causes severe blind spots.

J. D. Kang and M. B. Albin (2000) formulated the bivariate $T^2$ profile control chart:

where $\boldsymbol{\beta}_0 = (A_0, A_1)^T$ is the in-control parameter vector, and:

Simultaneously, the residual variance $\sigma^2$ is monitored using a sample variance chart based on $MSE_i / \sigma^2 \sim \chi^2_{n-2} / (n - 2)$. If the stamping die wears out or cracking occurs, the torque-angle slope shifts or residual curvature develops, triggering an immediate profile anomaly alarm.

Classical Shewhart, CUSUM, and EWMA control charts employ rigid, fixed sampling intervals $h$ (e.g., take a sample of size $n = 5$ every $h = 1.0\text{ hour}$). This static cadence is statistically inefficient: when the monitored sample mean falls very close to the center line (indicating strong evidence of stability), the process is inspected just as frequently as when the sample mean plots dangerously close to the warning limit.

To maximize statistical responsiveness while preserving long-term inspection budgets, Marion R. Reynolds et al. (1988) formulated Adaptive Control Charts:

Variable Sampling Interval (VSI) Charts

In a VSI control chart, the inter-sample inspection interval $h$ switches dynamically between two values:

• $h_1$ (Short sampling interval, e.g., $15\text{ minutes}$): Applied whenever the current sample statistic plots in a warning region ($w \le |Z_t| < 3$), indicating heightened suspicion of an incipient assignable cause.
• $h_2$ (Long sampling interval, e.g., $2.0\text{ hours}$): Applied whenever the current sample statistic plots in the central region ($|Z_t| < w$), confirming robust in-control operation.

Let $p_0 = P(|Z_t| < w \mid \text{in-control})$ be the probability that a sample falls in the central region under in-control conditions. The average in-control sampling interval is:

By setting $\bar{h}_0$ equal to the standard fixed sampling interval $h_{\text{fixed}}$, the VSI chart consumes the exact same long-term inspection budget as a classical chart, but reduces the Average Time to Signal ($ATS_1$) for small and moderate shifts by over $50\%$!

Variable Sample Size (VSS) Charts

Similarly, Variable Sample Size (VSS) charts switch the subgroup size between a small sample $n_1$ (when $|Z_t| < w$) and a large sample $n_2$ (when $w \le |Z_t| < 3$). Combining both approaches yields the Variable Sample Size and Sampling Interval (VSSI) chart, achieving near-optimal sequential detection without continuous-time sensor hardware.