Process Monitoring#

Class for ControlChart: robust control charts with a balance between CUSUM and Shewhart properties.

process_improve.monitoring.control_charts.BIWEIGHT_RHO_CONSISTENCY = 3.266736#

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.

Departure from the paper, kept deliberately: equation (13) of Gelper, Fried and Croux (2010, p 289) pairs the cutoff k = 2 with c_k = 2.52 (“for the common choice of k = 2 we have c_k = 2.52”; numerically 1 / E[rho_norm(Z)] = 2.5153 for k = 2). This module uses 2.52 as the cutoff instead, with the matching c_k above. Both pairs give consistent scale estimates on Gaussian data, but this biweight saturates later (at |x| > 2.52 rather than |x| > 2), so it down-weights large forecast errors less than the paper’s. The constants are left unchanged because existing results and callers depend on them.

process_improve.monitoring.control_charts.rho(x, k=2.52)[source]#

Bi-weight rho function.

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).

Parameters:

x (float) – Value at which to evaluate the bi-weight rho function.
k (float, optional) – Bi-weight cutoff. The default of 2.52 is kept for compatibility. In the paper (equation 13, p 289) 2.52 is instead the consistency constant c_k for the cutoff k = 2; see the note on BIWEIGHT_RHO_CONSISTENCY.

Returns:

c_k * (1 - (1 - (x / k)**2)**3) for |x| <= k; for |x| > k the function saturates at the consistency constant c_k.

Return type:

float

References

https://onlinelibrary.wiley.com/doi/abs/10.1002/for.1125

process_improve.monitoring.control_charts.psi(x, k=2.0)[source]#

Pre-clean based on the Huber psi function.

Can be interpreted as replacing unexpected high or low values by a more likely value.

Parameters:

x (float) – Value to clean.
k (float, optional) – Huber tuning constant: the threshold beyond which x is clipped to k * sign(x). The default of 2.0 follows p 288 of the referenced paper.

Returns:

x itself when |x| < k; otherwise k * sign(x).

Return type:

float

References

https://onlinelibrary.wiley.com/doi/abs/10.1002/for.1125

process_improve.monitoring.control_charts.cusum_about_target(values, target)[source]#

Return the cumulative sum of the deviations from a target, as Box, Hunter and Hunter plot it.

Unlike the tabular CUSUM of ControlChart(variant="cusum"), it has no reference value and no decision interval: plotted against the sample number, its slope is the mean’s distance from the target, and a change of slope marks a change of mean (Box, Hunter and Hunter 1978, section 17.1).

Parameters:

values (Sequence[float] or pd.Series) – The observations, in the order they were made. A missing one (NaN) holds the sum.
target (float) – The target the deviations are taken from.

Returns:

One row per observation, with the index of values when it is a Series: y, the deviation from the target, and the cusum of the deviations so far.

Return type:

pd.DataFrame

Examples

>>> cusum_about_target([61.0, 59.5, 62.0], target=60.0)["cusum"].tolist()
[1.0, 0.5, 2.5]

class process_improve.monitoring.control_charts.ControlChart(style='robust', variant='HW')[source]#

Bases: object

Create control chart instance objects.

Parameters:

style (str)
variant (str)

ld_1: float | None#

each exists on an instance only once a caller pins it in calculate_limits or a fit sets it, which hasattr relies on. A Holt-Winters weight left as None is searched for.

Type:: Tuning parameters, declared but not assigned here

ld_2: float | None#

k: float#

h: float#

__init__(style='robust', variant='HW')[source]#

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): Four variants are accepted: 'hw' (the default),

'xbar.no.subgroup', 'ewma' and 'cusum'. Any other value 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 variants are:

‘xbar.no.subgroup’ [Shewhart chart, with no subgroups]. In other words, each observation is independently plotted on the control chart.

‘ewma’ [exponentially weighted moving average chart]. The charted statistic is z_t = ld_1 * y_t + (1 - ld_1) * z_{t-1}, started at the target, and it is compared with its exact time-varying limits, target +/- 3 s sqrt(ld_1 / (2 - ld_1) * (1 - (1 - ld_1) ** (2 t))). Pass the weight ld_1 (lambda, default 0.2) to calculate_limits. Averaging over recent samples makes it detect a small sustained shift that a Shewhart chart misses; ld_1 = 1 gives the Shewhart individuals chart.

‘cusum’ [tabular cumulative sum chart]. Two one-sided sums, C+_t = max(0, C+_(t-1) + (y_t - target) - k s) and C-_t = max(0, C-_(t-1) - (y_t - target) - k s), start at zero, and an alarm is raised when either exceeds the decision interval h s; that sum then restarts from zero. Pass k (the reference value: half the shift to detect, in units of s, default 0.5) and h (default 5) to calculate_limits. Each alarm also estimates when the shift started and the new mean (cusum_alarms).

Parameters:

style (str)
variant (str)

Return type:

None

calculate_limits(y, target=None, sigma=None, **kwargs)[source]#

Find for a given vector y, the control chart target and limits.

Works for the Holt-Winters (‘hw’), ‘xbar.no.subgroup’, ‘ewma’ and ‘cusum’ variants. For ‘ewma’, idx_outside_3S lists the samples where the EWMA statistic is outside its limits, and df gains the columns ewma, ewma_ucl and ewma_lcl.

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), there is no Holt-Winters fit: a target and standard deviation that are not given are the median and the standard deviation of the data. In every case, if target and sigma (a known standard deviation) are numeric, those values are used, and the fitted standard deviation s is sigma; if not, they are estimated. When the data carry no usable spread (a constant series, for example), a given sigma is also the scale of the Holt-Winters limits, and the smoothing weights that were not pinned are 0.5 and 0.8; without sigma such a series raises ValueError.

The keyword s is a deprecated alias for sigma.

Missing values (NaN) in y are bridged by imputing a one-step-ahead error equal to the median absolute error of the previous 10 rows. That error is never negative, so a gap pulls the Holt-Winters level upward, increasingly so for long gaps; fill or drop long runs of missing values before fitting. A gap that starts at index 1, the first row with a one-step-ahead error, has no earlier error to impute from: those rows carry the forecast forward instead (the level follows the trend, the trend and scale are held), and their error and cleaned value stay NaN.

Every call starts afresh: nothing fitted by an earlier call on the same chart (target, s, lambdas, warm-up statistics) is carried into this one.

Parameters:

y (ndarray | Series)
target (float | None)
sigma (float | None)

Return type:

None

EWMA_DEFAULT_WEIGHT: ClassVar[float] = 0.2#: Default EWMA weight, inside the 0.05 to 0.25 range usually recommended for small shifts.

CUSUM_DEFAULT_K: ClassVar[float] = 0.5#

k = 0.5 tunes the chart to a one-sigma shift, and h = 5 gives an in-control average run length of 465.

Type:: Default CUSUM reference value and decision interval, in units of s

CUSUM_DEFAULT_H: ClassVar[float] = 5.0#

process_improve.monitoring.metrics.calculate_cpk(df, which_column, specifications=(nan, nan), trim_percentile=2.5)[source]#

Calculate the process capability, Cpk, near either the lower or the upper limit [will be automatically determined which].

Process capability, nearer the lower limit = (avg - lower_spec)/(3 x std deviation) Process capability, nearer the upper limit = (upper_spec - avg)/(3 x std deviation)

Parameters:

df (pd.DataFrame) – Raw data, at least one column is numeric.
which_column (str) – Indicates which is the column of data that should be used for the Cpk calculation.
specifications (tuple of (lower, upper), optional) –
A 2-tuple (lower_spec, upper_spec) of the lower and upper specification limits. Each element may be:
- a numeric value, when the specification is constant over time;
- a string, interpreted as a column name in df whose values give the per-row specification (use this when the specification changes over time);
- None, in which case the corresponding spec is estimated from the data using trim_percentile (a percentile-based robust limit).
Default is (np.nan, np.nan), which treats both specs as numeric NaN and yields NaN for the corresponding side of the Cpk calculation.
trim_percentile (float, optional) – Controls two things. (1) When a specification limit is missing, trim_percentile is used as a percentile on the data (in percent) to estimate that limit: the lower spec is set to np.nanpercentile(data, trim_percentile) and the upper spec to np.nanpercentile(data, 100 - trim_percentile). Default 2.5 therefore yields the 2.5th and 97.5th percentiles. (2) When trim_percentile > 0 the centre/spread used in the Cpk formula switch from mean/std to robust alternatives (median and Sn); when 0 the classical mean/std are used.

Returns:

A bunch with the following fields:

cpk: the Cpk value (the limiting, i.e. smaller, of the two sides).
center: the center (mean or median) of the limiting side, i.e. of the distance-to-specification metric; the lower side when the two sides tie.
spread: the spread (standard deviation or Sn) of the limiting side.
rsd: the relative standard deviation OF THE DATA COLUMN, as a percentage: spread / center-of-the-data * 100. (Earlier versions divided by the distance-to-spec centre, which is not an RSD of anything and changed value when the spec limit moved.)
defects_per_million: the parts per million outside the specifications of a normal process with this capability, (Phi(-3 Cpk_lower) + Phi(-3 Cpk_upper)) x 10^6, 2 Phi(-3 Cp) x 10^6 when the process is centred; a side with no capability (NaN) adds nothing.
defects_per_million_nearer: those beyond the nearer specification alone, Phi(-3 Cpk) x 10^6.

Return type:

sklearn.utils.Bunch

Notes

The spread is estimated from ALL the values in the column (the overall long-term variation), not from a within-subgroup range or moving range. Strictly this makes the statistic a performance index (Ppk-style) rather than a short-term capability index (Cpk with sigma-hat = Rbar/d2); on a stable process the two coincide, and on an unstable one this value is the more honest of the two.

Evolutionary operation#

Evolutionary operation (EVOP): the calculation sheet of a \(2^2\) design with a centre point.

In evolutionary operation a full-scale process runs, cycle after cycle, at five conditions close to its usual one: the centre (point 1) and the corners of a small square in two factors A and B (point 2 with A and B both low, 3 both high, 4 A high and B low, 5 A low and B high). The changes are small enough to keep the product saleable, and repeating the cycle lets their effects emerge from the noise. After each cycle the calculation sheet updates:

each point’s sum and average over the cycles so far;
the effects of A and B, their interaction, and the change in mean (the corners against the centre);
an estimate of the noise, from the range of the differences between each point’s previous average and its new reading, times \(f_{k,n} = \sqrt{(n-1)/n} / d_2(k)\) (see evop_f());
the two-standard-deviation error limits of the averages, the effects and the change in mean.

The method and its sheet are Box and Draper’s (Evolutionary Operation, 1969); Montgomery works an example in section 11.7.

Examples

Montgomery’s Example 11.7, three cycles of yields at points 1 to 5, the noise factor rounded to two decimals as his Table 11.27 prints it:

>>> from process_improve.monitoring import evop_worksheet
>>> cycles = [[84.5, 84.2, 84.9, 84.5, 84.3], [84.9, 84.6, 85.9, 83.5, 84.0], [85.0, 84.0, 86.6, 84.9, 85.2]]
>>> board = evop_worksheet(cycles, factors=("T", "P"), range_factor_decimals=2)[-1]
>>> {name: round(effect, 2) for name, effect in board["effects"].items()}
{'T': 0.67, 'P': 0.87, 'T:P': 0.63, 'change_in_mean': -0.07}
>>> round(board["average_s"], 2), {name: round(limit, 2) for name, limit in board["limits"].items()}
(0.58, {'averages': 0.67, 'effects': 0.67, 'change_in_mean': 0.6})

process_improve.monitoring.evop.POINTS = ((0, 0), (-1, -1), (1, 1), (1, -1), (-1, 1))#

the centre, then the corners as (A, B) in coded units.

Type:: The points of the design in the sheet’s order

process_improve.monitoring.evop.CHANGE_IN_MEAN_LIMIT = 1.7888543819998317#: Two standard deviations of the change in mean, (y2 + y3 + y4 + y5 - 4 y1) / 5, in units of sigma / sqrt(n): its variance is (4 + 16) sigma^2 / (25 n). Box and Draper’s sheets print 1.78.

process_improve.monitoring.evop.range_d2(k)[source]#

Return \(d_2(k)\), the mean range of k independent standard normal observations.

It is \(\int_{-\infty}^{\infty} [1 - \Phi(x)^k - (1 - \Phi(x))^k]\,dx\): 1.128 for two observations, 2.326 for five. A range of k observations divided by it estimates their standard deviation.

Parameters:: k (int) – The number of observations in the range, 2 or more.
Returns:: The mean range, in units of the observations’ standard deviation.
Return type:: float
Raises:: ValueError – If k is not a whole number of at least 2.

process_improve.monitoring.evop.evop_f(k, n, decimals=None)[source]#

Return \(f_{k,n} = \sqrt{(n-1)/n} / d_2(k)\), which turns the range of a cycle’s differences into a sigma.

In cycle \(n\), each of the k differences between a point’s previous average and its new reading has variance \(\sigma^2 (1 + 1/(n-1)) = \sigma^2 n/(n-1)\). Their range, divided by \(d_2(k)\) and multiplied by \(\sqrt{(n-1)/n}\), estimates \(\sigma\).

Montgomery’s Table 11.27 prints it to two decimals, and matches that rounding in all but two entries: \(f_{5,9}\) and \(f_{10,3}\), 0.4054 and 0.2653, print as 0.40 and 0.26. All of its entries are reproduced with \(d_2\) rounded to two decimals and \(\sqrt{(n-1)/n}\) to three.

Parameters:

k (int) – The number of points in the design, whose differences give the range: 5 for a \(2^2\) design with a centre point.
n (int) – The cycle, 2 or more: the first cycle has no previous average to differ from.
decimals (int or None) – Round \(f_{k,n}\) to this many decimals, as a printed table does; None leaves it unrounded.

Returns:

\(f_{k,n}\).

Return type:

float

Raises:

ValueError – If n is not a whole number of at least 2, or k is not one of at least 2.

process_improve.monitoring.evop.evop_worksheet(cycles, factors=('A', 'B'), range_factor_decimals=None)[source]#

Work the EVOP calculation sheet of a \(2^2\) design with a centre point, one cycle after another.

Parameters:

cycles (Sequence[Sequence[float]]) – The readings, a row per cycle in the order they were run, each row the five points in the order of POINTS: the centre, then A and B both low, both high, A high and B low, A low and B high.
factors (tuple[str, str]) – The names of the two factors, A and B, which name the effects.
range_factor_decimals (int or None) – Round the factor that turns a range into a standard deviation, Box and Draper’s \(f_{k,n}\), to this many decimals, as a printed table does; None uses it unrounded.

Returns:

A sheet per cycle, the last of them the information board:

cycle, \(n\); readings, the cycle’s five readings;
sums and averages of each point over the cycles so far;
effects: the main effects of A and B, (y3 + y4 - y2 - y5) / 2 and (y3 + y5 - y2 - y4) / 2, their interaction "A:B", (y2 + y3 - y4 - y5) / 2, and change_in_mean, (y2 + y3 + y4 + y5 - 4 y1) / 5, all from the averages;
from the second cycle: differences, each point’s previous average less its new reading, their range, s (the range times \(f_{5,n}\)), the running sum_s of these and its average average_s over the cycles that gave one, and limits, two standard deviations of the averages and the effects, \(2 S/\sqrt{n}\), and of the change_in_mean, \(1.79\,S/\sqrt{n}\), with \(S\) the average. In the first cycle these are None.

Return type:

list[dict[str, Any]]

Raises:

ValueError – If there is no cycle, a cycle does not have five readings, a reading is not a finite number, or the factors are not two different names (neither of them "change_in_mean").