Statistical Process Control — Control Charts

Statistical process control (SPC) uses control charts to monitor whether a process remains in a state of statistical control — stable over time with only common-cause variation present.

The Control Chart Framework

X-bar Chart (Process Mean)

R Chart (Process Variability)

p Chart (Proportion Defective)

c Chart (Count of Defects)

Western Electric Rules

Process Capability Indices

means the process variation fits within the specification width. means 6σ exactly equals the specification range.

Average Run Length (ARL)

CUSUM Chart

EWMA Chart

Python Implementation

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(42)

# --- Simulate in-control and out-of-control process ---
n_subgroups = 50
subgroup_size = 5
mu0 = 50.0
sigma = 2.0

# In-control: first 30 subgroups
X_in = np.random.randn(30, subgroup_size) * sigma + mu0
# Out-of-control: shift of 1.5σ after subgroup 30
X_out = np.random.randn(20, subgroup_size) * sigma + (mu0 + 1.5 * sigma)
X = np.vstack([X_in, X_out])

xbar = X.mean(axis=1)
R = X.max(axis=1) - X.min(axis=1)
R_bar = R[:30].mean()  # Use in-control phase for limits
xbar_bar = xbar[:30].mean()

# --- X-bar chart ---
A2 = 0.577  # For n=5
UCL_xbar = xbar_bar + A2 * R_bar
LCL_xbar = xbar_bar - A2 * R_bar

# --- R chart ---
D4 = 2.114  # For n=5
D3 = 0      # For n=5
UCL_R = D4 * R_bar
LCL_R = D3 * R_bar

# --- Western Electric rules detection ---
def detect_western_electric(xbar, center, ucl, lcl):
    sigma_est = (ucl - center) / 3
    signals = np.zeros(len(xbar), dtype=bool)
    for i in range(len(xbar)):
        # Rule 1: Beyond 3σ
        if xbar[i] > ucl or xbar[i] < lcl:
            signals[i] = True
        # Rule 2: 2 of 3 beyond 2σ
        if i >= 2:
            above_2s = np.sum(xbar[i-2:i+1] > center + 2*sigma_est) >= 2
            below_2s = np.sum(xbar[i-2:i+1] < center - 2*sigma_est) >= 2
            if above_2s or below_2s:
                signals[i] = True
        # Rule 3: 4 of 5 beyond 1σ
        if i >= 4:
            above_1s = np.sum(xbar[i-4:i+1] > center + sigma_est) >= 4
            below_1s = np.sum(xbar[i-4:i+1] < center - sigma_est) >= 4
            if above_1s or below_1s:
                signals[i] = True
        # Rule 4: 8 consecutive same side
        if i >= 7:
            if np.all(xbar[i-7:i+1] > center) or np.all(xbar[i-7:i+1] < center):
                signals[i] = True
    return signals

signals = detect_western_electric(xbar, xbar_bar, UCL_xbar, LCL_xbar)

# --- CUSUM chart ---
def cusum_chart(x, mu0, sigma, K_factor=0.5, H=5):
    K = K_factor * sigma
    S_pos = np.zeros(len(x) + 1)
    S_neg = np.zeros(len(x) + 1)
    for t in range(len(x)):
        S_pos[t+1] = max(0, S_pos[t] + (x[t] - mu0) - K)
        S_neg[t+1] = max(0, S_neg[t] + (mu0 - x[t]) - K)
    return S_pos[1:], S_neg[1:]

S_pos, S_neg = cusum_chart(xbar, mu0, sigma, K_factor=0.5, H=5*sigma)

# --- EWMA chart ---
def ewma_chart(x, mu0, sigma, lam=0.2, L=3):
    z = np.zeros(len(x) + 1)
    z[0] = mu0
    for t in range(len(x)):
        z[t+1] = lam * x[t] + (1 - lam) * z[t]
    z = z[1:]
    ucl = mu0 + L * sigma * np.sqrt(lam / (2 - lam))
    lcl = mu0 - L * sigma * np.sqrt(lam / (2 - lam))
    return z, ucl, lcl

ewma_z, ewma_ucl, ewma_lcl = ewma_chart(xbar, mu0, sigma, lam=0.2)

# --- Process capability ---
USL = mu0 + 3 * sigma
LSL = mu0 - 3 * sigma
Cp = (USL - LSL) / (6 * sigma)
Cpk = min((USL - mu0) / (3 * sigma), (mu0 - LSL) / (3 * sigma))
print(f"Cp = {Cp:.3f}, Cpk = {Cpk:.3f}")

# --- Plotting ---
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# X-bar chart
axes[0, 0].plot(range(1, n_subgroups+1), xbar, 'b-o', markersize=4)
axes[0, 0].axhline(UCL_xbar, color='red', linestyle='--', label='UCL')
axes[0, 0].axhline(xbar_bar, color='green', linestyle='-', label='CL')
axes[0, 0].axhline(LCL_xbar, color='red', linestyle='--', label='LCL')
axes[0, 0].axvline(30.5, color='orange', linestyle=':', label='Shift')
if signals.any():
    axes[0, 0].scatter(np.where(signals)[0]+1, xbar[signals], color='red', s=80, zorder=5, label='Signal')
axes[0, 0].set_title('X-bar Chart')
axes[0, 0].set_xlabel('Subgroup')
axes[0, 0].set_ylabel('X-bar')
axes[0, 0].legend()

# R chart
axes[0, 1].plot(range(1, n_subgroups+1), R, 'b-o', markersize=4)
axes[0, 1].axhline(UCL_R, color='red', linestyle='--', label='UCL')
axes[0, 1].axhline(R_bar, color='green', linestyle='-', label='CL')
axes[0, 1].axhline(LCL_R, color='red', linestyle='--', label='LCL')
axes[0, 1].set_title('R Chart')
axes[0, 1].set_xlabel('Subgroup')
axes[0, 1].set_ylabel('Range')
axes[0, 1].legend()

# CUSUM
axes[1, 0].plot(range(1, n_subgroups+1), S_pos, 'b-', label='CUSUM⁺')
axes[1, 0].plot(range(1, n_subgroups+1), S_neg, 'r-', label='CUSUM⁻')
axes[1, 0].axhline(5*sigma, color='gray', linestyle='--', label='Decision interval H')
axes[1, 0].axvline(30.5, color='orange', linestyle=':', label='Shift')
axes[1, 0].set_title('CUSUM Chart')
axes[1, 0].set_xlabel('Subgroup')
axes[1, 0].set_ylabel('CUSUM statistic')
axes[1, 0].legend()

# EWMA
axes[1, 1].plot(range(1, n_subgroups+1), ewma_z, 'b-', linewidth=2, label='EWMA')
axes[1, 1].axhline(ewma_ucl, color='red', linestyle='--', label='UCL')
axes[1, 1].axhline(mu0, color='green', linestyle='-', label='CL')
axes[1, 1].axhline(ewma_lcl, color='red', linestyle='--', label='LCL')
axes[1, 1].axvline(30.5, color='orange', linestyle=':', label='Shift')
axes[1, 1].set_title('EWMA Chart (λ=0.2)')
axes[1, 1].set_xlabel('Subgroup')
axes[1, 1].set_ylabel('EWMA statistic')
axes[1, 1].legend()

plt.tight_layout()
plt.savefig('spc_control_charts.png', dpi=150)
plt.show()

# --- ARL calculation via simulation ---
def simulate_arl(mu_shift, sigma, mu0, n=5, n_sims=10000, k=3):
    arl_count = 0
    for _ in range(n_sims):
        t = 0
        while True:
            t += 1
            x = np.random.randn(n) * sigma + mu0 + mu_shift
            xbar_t = x.mean()
            R_t = x.max() - x.min()
            ucl = mu0 + k * sigma / np.sqrt(n)
            lcl = mu0 - k * sigma / np.sqrt(n)
            if xbar_t > ucl or xbar_t < lcl:
                arl_count += t
                break
            if t > 10000:
                arl_count += t
                break
    return arl_count / n_sims

print("\n=== Average Run Length (X-bar chart, n=5) ===")
for shift in [0, 0.5, 1.0, 1.5, 2.0, 3.0]:
    arl = simulate_arl(shift * sigma, sigma, mu0)
    print(f"Shift = {shift}σ: ARL ≈ {arl:.1f}")

Related Topics

• See Six Sigma for DMAIC methodology
• See Design of Experiments for process optimization
• See Hypothesis Testing for the inference framework underlying control charts

Key Takeaways