"""Class for ControlChart: robust control charts with a balance between CUSUM and Shewhart properties."""
import logging
from typing import ClassVar
import numpy as np
import pandas as pd
from ..regression.methods import repeated_median_slope
from ..univariate.metrics import median_absolute_deviation
logger = logging.getLogger(__name__)
#: Consistency constant for the bounded biweight rho with cutoff k = 2.52,
#: c_k = 1 / E[rho_norm(Z)] for Z ~ N(0, 1), where rho_norm is the biweight
#: rho normalised to a maximum of 1. Computed via numerical integration
#: (scipy.integrate.quad of rho_norm(z) * phi(z); E = 0.3061160, so
#: c_k = 3.266736). This makes E[rho(Z)] = 1, the condition for the scale
#: estimates built from rho to be consistent for sigma on Gaussian data.
BIWEIGHT_RHO_CONSISTENCY = 3.266736
[docs]
def rho(x: float, k: float = 2.52) -> float:
"""
Bi-weight rho function.
Fixed cutoff of k=2.52 is from p 289 of the paper
https://onlinelibrary.wiley.com/doi/abs/10.1002/for.1125
The multiplier is the consistency constant c_k (chosen so that
``E[rho(Z)] = 1`` for standard-normal Z), NOT the cutoff k. The paper
treats the two as separate constants; an earlier version of this code
conflated them and used k = 2.52 as the multiplier, which made every
scale estimate derived from rho a factor ``sqrt(2.52 * 0.30612) = 0.878``
too small, i.e. +/-3S control limits that were really +/-2.63 sigma
(a ~3x inflation of the false-alarm rate).
"""
c_k = BIWEIGHT_RHO_CONSISTENCY
return c_k if np.abs(x) > k else c_k * (1 - np.power(1 - np.power(x / k, 2), 3))
def _finite(values: pd.Series | np.ndarray) -> bool:
"""Return ``True`` if at least one entry of ``values`` is finite."""
return bool(np.isfinite(np.asarray(values, dtype=float)).any())
def _training_error_radicand(future_errors: pd.Series, rho_func: np.vectorize) -> float:
"""
Return tau squared for the training-sample errors, or NaN when it is undefined.
Equation 16 of the Holt-Winters paper. This is the shared computation behind the
lambda grid search, which treats an undefined cell as a non-contender, and
:func:`_tau_from_training_errors`, which rejects it. Returning NaN rather than
routing an all-NaN slice through ``np.nanmedian`` / ``np.nanmean`` keeps a
``RuntimeWarning`` from escaping ``calculate_limits``. (#557)
Parameters
----------
future_errors : pd.Series
One-step-ahead errors over the training samples (those after the warm-up).
rho_func : np.vectorize
Element-wise bounded loss applied to the standardised errors.
Returns
-------
float
The radicand, or ``np.nan`` when the errors hold nothing finite or carry no
usable spread.
"""
if not _finite(future_errors):
return float("nan")
s_t_median_error = 1.48 * future_errors.abs().median()
if not np.isfinite(s_t_median_error) or s_t_median_error <= 0:
return float("nan")
return float(np.power(s_t_median_error, 2) * np.nanmean(rho_func(future_errors / s_t_median_error)))
def _tau_from_training_errors(future_errors: pd.Series, rho_func: np.vectorize) -> float:
"""
Return the robust scale estimate (tau) for the training-sample errors.
Equation 16 of the Holt-Winters paper, with each way the estimate can come out
undefined reported rather than silently reduced to zero. The previous
``np.sqrt(max(0.0, resids))`` was written to stop a negative radicand reaching
``sqrt``, but ``max(0.0, nan)`` returns ``0.0``: ``nan > 0.0`` is ``False``, so
``max`` keeps its first argument. A NaN therefore became a zero scale, which the
caller carried outward as control limits of zero width. (#557)
Parameters
----------
future_errors : pd.Series
One-step-ahead errors over the training samples (those after the warm-up).
rho_func : np.vectorize
Element-wise bounded loss applied to the standardised errors.
Returns
-------
float
The scale estimate, strictly positive.
Raises
------
ValueError
If no training error is finite, or if the resulting radicand is non-finite
or not positive (the errors carry no usable spread).
"""
if not _finite(future_errors):
# Missing values at or near the start of the series propagate through the
# Holt-Winters recursion and leave every training error NaN.
raise ValueError(
f"The control-chart scale estimate is undefined: none of the {future_errors.size} "
"training-sample errors is finite. Missing values at or near the start of the "
"series propagate through the Holt-Winters recursion and leave every error NaN. "
"Supply a series without leading gaps, or pass an explicit positive 's'."
)
resids = _training_error_radicand(future_errors, rho_func)
if not np.isfinite(resids) or resids <= 0:
raise ValueError(
f"The control-chart scale estimate is undefined (tau^2 = {resids}), so the control "
"limits would have zero width. Supply more representative data, or pass an "
"explicit positive 's'."
)
return float(np.sqrt(resids))
[docs]
def psi(x: float, k: float = 2.0) -> float:
"""
Pre-clean based on the Huber psi function.
Can be interpreted as replacing unexpected high or low values by a more likely value.
From p 288 of the paper
https://onlinelibrary.wiley.com/doi/abs/10.1002/for.1125
"""
return x if abs(x) < k else k * np.sign(x)
[docs]
class ControlChart:
"""Create control chart instance objects."""
[docs]
def __init__(self, style: str = "robust", variant: str = "HW") -> None:
"""
Create/initialize a control chart.
Args: style (str, optional): Which style control chart to calculate. Defaults to "robust".
Other choice is 'regular' (i.e. not-robust) calculations. User should then ensure that
no outliers are present in the data.
variant (str, optional): Only two variants are currently accepted:
``'hw'`` (the default) and ``'xbar.no.subgroup'``. Any other value,
including ``'cusum'``, raises ``ValueError`` at construction time.
The variant string is compared case-insensitively (it is normalised
via ``.strip().lower()`` on assignment), so ``'HW'``, ``'hw'``, and
``'Hw'`` are all equivalent.
The default is a Holt-Winters (`'hw'`) chart, with automatic determination of
control chart parameters. This chart is a blend of infinite history (CUSUM)
charts, and an instantaneous (no history taken into account) Shewhart chart. The
exact blend is specified by parameters `ld_1` (lambda 1) and `ld_2` (lambda 2).
The other accepted variant is:
'xbar.no.subgroup' [Shewhart chart, with no subgroups]. In other words, each
observation is independently plotted on the control chart.
A pure 'cusum' (CUmulative SUM) chart is a planned future variant but
is not currently implemented; passing ``variant='cusum'`` raises
``ValueError``. The Holt-Winters ('hw') default already blends
CUSUM-style infinite history with Shewhart-style instantaneous
behaviour via its lambda parameters.
"""
self.style = style.strip()
self.variant = variant.strip().lower()
# An unknown variant previously slipped through every fit branch and
# surfaced much later as a misleading "input is likely constant or too
# short" error from calculate_limits. Reject it up front instead.
_supported = {"hw", "xbar.no.subgroup"}
if self.variant not in _supported:
raise ValueError(
f"Control chart variant {variant!r} is not implemented; supported variants are {sorted(_supported)}. "
"(A standalone CUSUM chart is a possible future variant; the 'hw' chart already blends "
"CUSUM-style history with Shewhart behaviour.)"
)
# Will be calculated by the self.calculate_limits() function
self.target: float | None = None
self._given_target: float | None = None
self._given_s: float | None = None
self.s: float | None = None
# index of elements which are found to be outside +/- 3S
self.idx_outside_3S: list[int] = []
self.warm_up: dict[str, float | np.ndarray | pd.Series] = {}
self.warm_up_M: int = 0
columns = [
"y",
"psi_input",
"rho_input",
"y_star",
"alpha_hat",
"beta_hat",
"sigma_hat",
"error",
]
self.df = pd.DataFrame(columns=columns, dtype=np.float64)
[docs]
def calculate_limits( # noqa: C901 - branch count is mostly simple input-validation guard clauses
self,
y: np.ndarray | pd.Series,
target: float | None = None,
s: float | None = None,
**kwargs,
) -> None:
"""
Find for a given vector `y`, the control chart target and limits.
Works for both the Holt-Winters ('hw') and 'xbar.no.subgroup' variants.
For the Holt-Winters variant, when there are fewer than
min(20, max(10, np.ceil(0.10 * N)))
measurements (where N is the length of the input vector), the target and
standard deviation are estimated directly from the data and any provided
`target` / `s` are ignored for that small-sample case.
Otherwise, if `target` and `s` are numeric, those values are used;
if not, they are estimated.
"""
self._given_target = target
self._given_s = s
logger.debug(
"ControlChart.calculate_limits: variant=%s, style=%s, given target=%s, s=%s",
self.variant,
self.style,
target,
s,
)
if s is not None:
self.s = float(s)
if not 0.0 < s < 1e300:
raise ValueError(
f"The given standard deviation must be positive and not excessively large (0 < s < 1e300); got {s}."
)
if target is not None:
self.target = float(target)
self.df["y"] = y.ravel() if isinstance(y, np.ndarray) else pd.Series(y).values.ravel()
self.N = self.df.shape[0]
# Between M = 10 and 20 samples required to warm-up (calculate summary statistics)
self.warm_up["M"] = self.warm_up_M = int(min(20, max(10, np.ceil(0.10 * self.N))))
if (self.warm_up_M > self.N) and self.variant.strip().lower() == "hw":
# TO CHECK: Completely handle the case with very few samples. Is everything filled in?
# Also check case when some of these samples are NAN, you might have even fewer still.
self.target = self._target_calculated_best = self.df["y"].median()
self.s = self._tau = self.df["y"].std()
self.df["y_star"] = self.df["y"].values
self.df["alpha_hat"] = self.target
self.df["beta_hat"] = 0
self.df["sigma_hat"] = np.nan
self.df["error"] = np.nan
return
# Check if there are enough training samples:
if 2 * self.warm_up_M > self.N:
self.train_samples: list[int] = [int(i) for i in np.arange(0, self.N)]
else:
self.train_samples = [int(i) for i in np.arange(self.warm_up_M, self.N)]
self._apply_tuning_kwargs(kwargs)
if self.variant.strip().lower() == "hw":
if not hasattr(self, "ld_1"):
self.ld_1 = None
if not hasattr(self, "ld_2"):
self.ld_2 = None
self._holt_winters_parameter_fit()
if self.variant.strip().lower() == "xbar.no.subgroup":
self._xbar_no_subgroup_fit()
# After whichever fit is completed, check which are outside +/- 3S.
# Explicit validation (not assert) so the guard survives `python -O`
# and surfaces as a documented ValueError at the tool boundary (SEC-17).
if self.target is None or self.s is None:
raise ValueError(
"Control chart limits could not be estimated; the input is likely "
"constant or too short to fit the chosen variant."
)
idx_bool = (self.df["y"] - self.target).abs() > 3.0 * self.s
self.idx_outside_3S = np.nonzero(idx_bool.to_numpy())[0].tolist()
#: Instance attributes a caller may pin via ``calculate_limits(**kwargs)``:
#: the Holt-Winters smoothing lambdas. Everything else is internal state.
_TUNING_KWARGS: ClassVar[frozenset[str]] = frozenset({"ld_1", "ld_2"})
def _apply_tuning_kwargs(self, kwargs: dict[str, object]) -> None:
"""Set caller-pinned tuning parameters, rejecting anything off the allowlist.
A blanket ``setattr(self, key, val)`` over ``**kwargs`` would let a caller
silently overwrite internal state (``self.s``, ``self.target``,
``self.train_samples``, even a bound method) and would swallow typos. We
therefore accept only the documented Holt-Winters smoothing lambdas and
raise a clear ``ValueError`` otherwise.
"""
unknown = set(kwargs) - self._TUNING_KWARGS
if unknown:
raise ValueError(
f"calculate_limits() got unexpected keyword argument(s) {sorted(unknown)}; "
f"only {sorted(self._TUNING_KWARGS)} (Holt-Winters smoothing lambdas) are accepted."
)
for key, val in kwargs.items():
setattr(self, key, val)
def _xbar_no_subgroup_fit(self) -> None:
"""
Fit the control chart from the data samples, assuming each sample is its own subgroup.
The `style` attribute ('regular' | 'robust') switches how the average and standard
deviation are calculated.
Control chart limits assume the data are normally distributed and independent. In
particular, this last assumption can have consequences if not actually met. Limits may be
too wide, or too narrow.
"""
if self.style == "regular":
self.target = self.df["y"].mean()
self.s = self.df["y"].std()
elif self.style == "robust":
self.target = self.df["y"].median()
self.s = (self.df["y"] - self.target).abs().median() * 1.4826
def _holt_winters_parameter_fit(self) -> None:
"""
Recommended in the paper: not to fit the lambda_s value, but to use a grid search for the
lambda_1 and lambda_2 values. This is done in a 5x5 grid in the code below.
"""
self.ld_s = ld_s = 0.2
rho_func = np.vectorize(rho)
# ``is not None``: an explicit ld_1=0.0 (or ld_2=0.0) is a legitimate
# user choice, but 0.0 is falsy and a plain truthiness test silently
# discarded it and ran the grid search instead.
if self.ld_1 is not None and self.ld_2 is not None:
# User has provided their own lambda_1 and lambda_2 values.
for _name, _val in (("Lambda_1", self.ld_1), ("Lambda_2", self.ld_2), ("Lambda_s", self.ld_s)):
if _val < 0.0:
raise ValueError(f"{_name} must be greater than or equal to zero.")
if _val > 1.0:
raise ValueError(f"{_name} must be less than or equal to 1.0.")
self._holt_winters_warmup_fit(ld_1=self.ld_1, ld_2=self.ld_2, ld_s=self.ld_s)
else:
# User wants to find an value for ld_1 and ld_2 that best fits the data
ld_1_index = np.linspace(0.1, 0.9, num=5, endpoint=True)
ld_2_index = np.linspace(0.1, 0.9, num=5, endpoint=True)
residuals, _ = np.meshgrid(ld_1_index, ld_2_index)
for i, ld_1 in enumerate(ld_1_index):
for j, ld_2 in enumerate(ld_2_index):
self._holt_winters_warmup_fit(ld_1=ld_1, ld_2=ld_2, ld_s=ld_s)
# Apply equation 16 from the paper to the residuals in the 'training' period,
# that is the samples after the warm-up period. NaN-aware
# statistics are required: row 0 never receives an "error"
# value, and for small samples (2 * warm_up_M > N) the
# training window includes row 0. With plain
# np.median/np.average every grid cell became NaN and the
# search silently "chose" (0.1, 0.1) via argmin-of-NaN.
future_errors = self.df["error"].iloc[np.asarray(self.train_samples, dtype=int)]
# An unusable cell records NaN and cannot win the search; the
# `np.all(np.isnan(residuals))` check below still catches the case
# where every cell is unusable. (#557)
residuals[i, j] = _training_error_radicand(future_errors, rho_func)
if np.all(np.isnan(residuals)):
raise ValueError(
"The Holt-Winters lambda grid search produced no usable residuals; "
"the input is likely constant, too short, or entirely missing."
)
min_idx = np.nanargmin(residuals)
best_ld_1 = ld_1_index[np.unravel_index(min_idx, residuals.shape)[0]]
best_ld_2 = ld_2_index[np.unravel_index(min_idx, residuals.shape)[1]]
# Store the parameters that were calculated, even if a `target` or `s` were provided.
self.ld_1 = best_ld_1
self.ld_2 = best_ld_2
self._residuals_HW = residuals
self._holt_winters_warmup_fit(ld_1=best_ld_1, ld_2=best_ld_2, ld_s=ld_s)
# Common code for both branches of if-else above
future_errors = self.df["error"].iloc[np.asarray(self.train_samples, dtype=int)]
self._tau = _tau_from_training_errors(future_errors, rho_func)
if self.target is None:
# Estimate the target as the median of the y-star (cleaned) y-values
self.target = self.df["y_star"].median()
else:
self._target_calculated_best = self.df["y_star"].median()
if self.s is None:
self.s = self._tau # or an alternative: self.df["sigma_hat"] is approximately OK
# The "delta" emphasizes that it is the deviation from the target.
self._delta_UCL_3sigma = +3.0 * self._tau
self._delta_LCL_3sigma = -3.0 * self._tau
def _holt_winters_warmup_fit(self, ld_1: float = 0.5, ld_2: float = 0.8, ld_s: float = 0.2) -> None:
"""
See paper: https://onlinelibrary.wiley.com/doi/abs/10.1002/for.1125.
Calculates the Holt-Winters fitting and control chart parameters, for given values of the
smoothing parameters lambda_1 (how must local history for the level is used, with values
approaching 1.0 implying that less history is used), and lambda_2 (history for the trend
that is used, with lambda_2 approaching 1.0 implying that historical data is less
interesting), and lambda_s, a similar parameter for the moving variance of the sequence.
lambda_1 = ld_1 = 0.5 (default): value must be between 0 <= ld_1 <= 1.0
lambda_2 = ld_2 = 0.8 (default): value must be between 0 <= ld_2 <= 1.0
lambda_s = ld_s = 0.2 (default): value must be between 0 <= ld_2 <= 1.0, based on values
used in the paper, recommended on page 291.
The ideal lambda values (ld_1, ld_2, ld_s) can be found from a grid search.
"""
df = self.df
y_warm_up = df["y"].iloc[0 : self.warm_up_M]
self.warm_up["y_zero_robust"] = y_warm_up.median()
if isinstance(self.target, float):
self.warm_up["alpha_0"] = self.target
self.warm_up["beta_0"] = 0.0
else:
# p 290 of the paper, https://onlinelibrary.wiley.com/doi/abs/10.1002/for.1125
self.warm_up["beta_0"] = repeated_median_slope(np.arange(self.warm_up_M), y_warm_up.to_numpy())
self.warm_up["alpha_0"] = np.nanmedian(y_warm_up - self.warm_up["beta_0"] * np.arange(self.warm_up_M))
if isinstance(self.s, float):
self.warm_up["sigma_0"] = self.s
else:
# p 290 of the paper, https://onlinelibrary.wiley.com/doi/abs/10.1002/for.1125
# The residual is y_t - alpha_0 - beta_0 * t (alpha_0 above is the
# median of exactly that de-trended series). Subtracting beta_0 as
# a constant, as an earlier version did, leaves the whole warm-up
# trend inside the residuals and inflates sigma_0 whenever the
# window drifts - precisely the situation this chart is for.
warm_up_residuals = y_warm_up - self.warm_up["alpha_0"] - self.warm_up["beta_0"] * np.arange(self.warm_up_M)
# Some other method that does not rely on SciPy for 1 function.
self.warm_up["sigma_0"] = median_absolute_deviation(np.asarray(warm_up_residuals), nan_policy="omit")
self.warm_up["residuals"] = warm_up_residuals
# A constant (zero-variance) warm-up window gives sigma_0 = MAD = 0, which
# would make rho/psi infinite and silently poison every downstream control
# limit with 0/NaN. Fail loudly instead of returning meaningless limits.
_residuals_array = warm_up_residuals.dropna().to_numpy()
_is_constant = _residuals_array.shape[0] == 0 or (_residuals_array[0] == _residuals_array).all()
if not _is_constant and self.warm_up["sigma_0"] == 0:
# Corner case: if there are multiple unique values in the warm-up residuals, but the MAD is zero,
# then sigma_0 is set to regular standard deviation instead, which is non-zero.
# This can happen when the warm-up residuals are symmetrically distributed around the median,
# leading to a MAD of zero, but still have variability that can be captured by the standard deviation.
self.warm_up["sigma_0"] = warm_up_residuals.std()
# A constant (zero-variance) warm-up window gives sigma_0 = MAD = 0, which
# would make rho/psi infinite and silently poison every downstream control
# limit with 0/NaN. Fail loudly instead of returning meaningless limits.
if not np.isfinite(self.warm_up["sigma_0"]) or self.warm_up["sigma_0"] <= 0:
raise ValueError(
"The Holt-Winters warm-up window has zero (or non-finite) variance "
"(sigma_0 = 0), so the control-chart limits would be undefined. "
"Supply more representative warm-up data, or pass a positive 's'."
)
df.loc[0, "rho_input"] = (
self.warm_up["y_zero_robust"] - self.warm_up["alpha_0"] - self.warm_up["beta_0"]
) / self.warm_up["sigma_0"]
df.loc[0, "psi_input"] = df["rho_input"][0]
df.loc[0, "y_star"] = self.warm_up["y_zero_robust"]
df.loc[0, "alpha_hat"] = self.warm_up["alpha_0"]
df.loc[0, "beta_hat"] = self.warm_up["beta_0"]
df.loc[0, "sigma_hat"] = self.warm_up["sigma_0"]
for i in range(1, self.N):
# Cover the warm-up period, and the rest of the data set. We need that for the residual
# calculation later anyway.
# Error = observed - predicted. Predicted = one-step-ahead prediction
error_i = df["y"][i] - (df["alpha_hat"][i - 1] + df["beta_hat"][i - 1])
if np.isnan(error_i):
# If there is an error, replace it with the median of the last 10 error estimates
# or as many points as available. When the gap is at the very start of the series
# there is no finite history to take that median over, so `error_i` stays NaN and
# propagates through alpha_hat / beta_hat / sigma_hat for the rest of the
# recursion; `calculate_limits` raises on that downstream. Filter to the finite
# history explicitly rather than letting an all-NaN slice reach pandas' median,
# which returns NaN but emits "Mean of empty slice". (#557)
recent_errors = df["error"].iloc[max(i - 10, 0) : i].abs().to_numpy(dtype=float)
recent_errors = recent_errors[np.isfinite(recent_errors)]
error_i = float(np.median(recent_errors)) if recent_errors.size else np.nan
rho_i = error_i / df["sigma_hat"][i - 1]
prior_variance = np.power(df["sigma_hat"][i - 1], 2)
sigma_i = np.sqrt(rho(rho_i) * ld_s * prior_variance + (1.0 - ld_s) * prior_variance)
psi_i = error_i / sigma_i
y_star_i = psi(psi_i) * sigma_i + df["alpha_hat"][i - 1] + df["beta_hat"][i - 1]
alpha_i = ld_1 * y_star_i + (1 - ld_1) * (df["alpha_hat"][i - 1] + df["beta_hat"][i - 1])
beta_i = ld_2 * (alpha_i - df["alpha_hat"][i - 1]) + (1 - ld_2) * df["beta_hat"][i - 1]
df.loc[i] = [
df["y"][i],
psi_i,
rho_i,
y_star_i,
alpha_i,
beta_i,
sigma_i,
error_i,
]
# Checks: for algorithm debugging:
# future_errors = df["error"]
# S_T_median_error = 1.48 * future_errors.abs().median() # must handle NaNs!
# resids = np.power(S_T_median_error, 2) * \
# np.nanmean((future_errors / S_T_median_error).apply(rho))
# print(np.sqrt(max(0.0, resids)))