u_analytics/capability/
percentile.rs

1//! Percentile-based process capability (ISO 22514-2).
2//!
3//! When process data is non-normal, capability indices based on mean and
4//! standard deviation can be misleading. The percentile method replaces the
5//! ±3σ spread with the distance between the 0.135 % and 99.865 % quantiles of
6//! the process distribution, and the mean with its median.
7//!
8//! Those quantiles are taken from a **distribution fitted to the data**, not
9//! from the data's own order statistics. With fewer than about 740 values the
10//! sample's 0.135 % point is its minimum (and the 99.865 % point its maximum),
11//! so an empirical spread is the sample range -- which is shorter than the
12//! process spread, and the shorter the smaller the sample: on normal data with
13//! a true Cp of 1.0 it reported about 1.37 at n = 50. The fitted model is the
14//! one [`fit_best`](crate::distribution::fit_best) ranks first by AIC, or the
15//! family the caller names.
16//!
17//! # Indices
18//!
19//! - **Cp\*** — (USL − LSL) / (X₉₉.₈₆₅ − X₀.₁₃₅)
20//! - **Cpk\*** — min(Cpu\*, Cpl\*)
21//! - **Cpu\*** — (USL − X₅₀) / (X₉₉.₈₆₅ − X₅₀)
22//! - **Cpl\*** — (X₅₀ − LSL) / (X₅₀ − X₀.₁₃₅)
23//!
24//! # References
25//!
26//! - ISO 22514-2, "Statistical methods in process management — Capability and
27//!   performance — Part 2: Process capability and performance of time-dependent
28//!   process models" (quantile method on an identified distribution).
29//! - Clements, J.A. (1989). "Process capability calculations for non-normal
30//!   distributions." *Quality Progress*, 22(9), 95–100 (quantiles of a fitted
31//!   distribution, there Pearson curves).
32
33use std::fmt;
34
35use crate::distribution::{
36    fit_best, fit_beta, fit_exponential, fit_gamma, fit_lognormal, fit_normal, fit_weibull,
37    fitted_quantile, FitResult,
38};
39
40/// Lower tail probability: the normal 3σ point, Φ(−3).
41const P_LOWER: f64 = 0.001_349_898_031_630_094_6;
42/// Upper tail probability: Φ(3).
43const P_UPPER: f64 = 1.0 - P_LOWER;
44
45/// Fewest values the method accepts.
46pub const MIN_PERCENTILE_DATA: usize = 20;
47
48/// The distribution family a percentile capability is computed from.
49#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
50#[non_exhaustive]
51pub enum CapabilityModel {
52    /// The family [`fit_best`] ranks first by AIC.
53    #[default]
54    Auto,
55    /// Normal.
56    Normal,
57    /// Log-normal (positive data).
58    LogNormal,
59    /// Two-parameter Weibull (positive data).
60    Weibull,
61    /// Gamma (positive data).
62    Gamma,
63    /// Exponential (positive data).
64    Exponential,
65    /// Beta (data within (0, 1)).
66    Beta,
67}
68
69impl CapabilityModel {
70    /// The fitted families, in the names [`FitResult::distribution`] uses.
71    fn of_fit(name: &str) -> Option<Self> {
72        Some(match name {
73            "Normal" => Self::Normal,
74            "LogNormal" => Self::LogNormal,
75            "Weibull" => Self::Weibull,
76            "Gamma" => Self::Gamma,
77            "Exponential" => Self::Exponential,
78            "Beta" => Self::Beta,
79            _ => return None,
80        })
81    }
82
83    pub(crate) fn fit(self, data: &[f64]) -> Option<FitResult> {
84        match self {
85            Self::Auto => fit_best(data)
86                .into_iter()
87                .find(|f| Self::of_fit(&f.distribution).is_some()),
88            Self::Normal => fit_normal(data),
89            Self::LogNormal => fit_lognormal(data),
90            Self::Weibull => fit_weibull(data),
91            Self::Gamma => fit_gamma(data),
92            Self::Exponential => fit_exponential(data),
93            Self::Beta => fit_beta(data),
94        }
95    }
96}
97
98/// Why a percentile capability was refused.
99#[derive(Debug, Clone, PartialEq)]
100#[non_exhaustive]
101pub enum PercentileCapabilityError {
102    /// Neither USL nor LSL was given.
103    NoSpecification,
104    /// USL is not above LSL.
105    InvalidSpecification,
106    /// Fewer than [`MIN_PERCENTILE_DATA`] values.
107    InsufficientData {
108        /// How many were given.
109        got: usize,
110    },
111    /// A value is NaN or infinite.
112    NonFiniteData {
113        /// Position of the first one.
114        index: usize,
115    },
116    /// Every value is the same.
117    NoVariation,
118    /// The named family cannot describe the data (log-normal, Weibull, gamma
119    /// and exponential need positive values, beta values within (0, 1)).
120    ModelNotFitted(CapabilityModel),
121}
122
123impl fmt::Display for PercentileCapabilityError {
124    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
125        match self {
126            Self::NoSpecification => {
127                write!(f, "at least one specification limit (USL or LSL) is required")
128            }
129            Self::InvalidSpecification => write!(f, "USL must be greater than LSL"),
130            Self::InsufficientData { got } => write!(
131                f,
132                "at least {MIN_PERCENTILE_DATA} data points are required for percentile capability, got {got}"
133            ),
134            Self::NonFiniteData { index } => {
135                write!(f, "data[{index}] is not a finite number")
136            }
137            Self::NoVariation => write!(f, "data has no variation"),
138            Self::ModelNotFitted(model) => {
139                write!(f, "the {model:?} distribution cannot be fitted to this data")
140            }
141        }
142    }
143}
144
145impl std::error::Error for PercentileCapabilityError {}
146
147/// Result of percentile-based process capability analysis.
148#[derive(Debug, Clone)]
149pub struct PercentileCapabilityResult {
150    /// Cp* = (USL − LSL) / (X₉₉.₈₆₅ − X₀.₁₃₅). Requires both limits.
151    pub cp_star: Option<f64>,
152    /// Cpk* = min(Cpu*, Cpl*). Requires at least one limit.
153    pub cpk_star: Option<f64>,
154    /// Cpu* = (USL − X₅₀) / (X₉₉.₈₆₅ − X₅₀). Requires USL.
155    pub cpu_star: Option<f64>,
156    /// Cpl* = (X₅₀ − LSL) / (X₅₀ − X₀.₁₃₅). Requires LSL.
157    pub cpl_star: Option<f64>,
158    /// Median of the fitted distribution, X₅₀.
159    pub median: f64,
160    /// 0.135 % quantile of the fitted distribution (lower natural process limit).
161    pub percentile_lower: f64,
162    /// 99.865 % quantile of the fitted distribution (upper natural process limit).
163    pub percentile_upper: f64,
164    /// The family the quantiles come from — never [`CapabilityModel::Auto`].
165    pub model: CapabilityModel,
166    /// The fit the quantiles come from: its parameters (as
167    /// [`FitResult::parameters`] names them) and likelihood. Pass it to
168    /// [`goodness_of_fit`](crate::distribution::goodness_of_fit) to judge
169    /// whether the family describes the data — the indices mean nothing when
170    /// it does not.
171    pub fit: FitResult,
172}
173
174/// Computes percentile-based process capability indices.
175///
176/// Fits `model` to `data` (or, for [`CapabilityModel::Auto`], the family
177/// [`fit_best`] ranks first) and takes the 0.135 %, 50 % and 99.865 %
178/// quantiles from it.
179///
180/// # Errors
181///
182/// [`PercentileCapabilityError`] — no limit, USL ≤ LSL, fewer than
183/// [`MIN_PERCENTILE_DATA`] values, a value that is not finite, no variation,
184/// or a named family the data does not admit.
185///
186/// # Examples
187///
188/// ```
189/// use u_analytics::capability::{percentile_capability, CapabilityModel};
190///
191/// let data: Vec<f64> = (0..100).map(|i| 10.0 + (i as f64) * 0.1).collect();
192/// let result = percentile_capability(&data, Some(5.0), Some(15.0), CapabilityModel::Auto).unwrap();
193/// assert!(result.cp_star.is_some());
194/// assert_ne!(result.model, CapabilityModel::Auto);
195/// ```
196pub fn percentile_capability(
197    data: &[f64],
198    lsl: Option<f64>,
199    usl: Option<f64>,
200    model: CapabilityModel,
201) -> Result<PercentileCapabilityResult, PercentileCapabilityError> {
202    if usl.is_none() && lsl.is_none() {
203        return Err(PercentileCapabilityError::NoSpecification);
204    }
205    if let (Some(u), Some(l)) = (usl, lsl) {
206        if u <= l {
207            return Err(PercentileCapabilityError::InvalidSpecification);
208        }
209    }
210    if data.len() < MIN_PERCENTILE_DATA {
211        return Err(PercentileCapabilityError::InsufficientData { got: data.len() });
212    }
213    if let Some(index) = data.iter().position(|v| !v.is_finite()) {
214        return Err(PercentileCapabilityError::NonFiniteData { index });
215    }
216    if data.iter().all(|&v| v == data[0]) {
217        return Err(PercentileCapabilityError::NoVariation);
218    }
219
220    let fit = model
221        .fit(data)
222        .ok_or(PercentileCapabilityError::ModelNotFitted(model))?;
223    let fitted = CapabilityModel::of_fit(&fit.distribution)
224        .ok_or(PercentileCapabilityError::ModelNotFitted(model))?;
225    let (p_lower, median, p_upper) =
226        quantiles(&fit).ok_or(PercentileCapabilityError::ModelNotFitted(model))?;
227
228    let spread = p_upper - p_lower;
229    if spread <= 0.0 {
230        return Err(PercentileCapabilityError::NoVariation);
231    }
232    let upper_spread = p_upper - median;
233    let lower_spread = median - p_lower;
234
235    let cp_star = match (usl, lsl) {
236        (Some(u), Some(l)) => Some((u - l) / spread),
237        _ => None,
238    };
239    let cpu_star = usl.and_then(|u| (upper_spread > 0.0).then(|| (u - median) / upper_spread));
240    let cpl_star = lsl.and_then(|l| (lower_spread > 0.0).then(|| (median - l) / lower_spread));
241    let cpk_star = match (cpu_star, cpl_star) {
242        (Some(u), Some(l)) => Some(u.min(l)),
243        (one, other) => one.or(other),
244    };
245
246    Ok(PercentileCapabilityResult {
247        cp_star,
248        cpk_star,
249        cpu_star,
250        cpl_star,
251        median,
252        percentile_lower: p_lower,
253        percentile_upper: p_upper,
254        model: fitted,
255        fit,
256    })
257}
258
259/// The 0.135 %, 50 % and 99.865 % quantiles of a fitted family.
260fn quantiles(fit: &FitResult) -> Option<(f64, f64, f64)> {
261    let q = |p| fitted_quantile(fit, p);
262    let (lo, mid, hi) = (q(P_LOWER)?, q(0.5)?, q(P_UPPER)?);
263    (lo.is_finite() && mid.is_finite() && hi.is_finite()).then_some((lo, mid, hi))
264}
265
266// ---------------------------------------------------------------------------
267// Tests
268// ---------------------------------------------------------------------------
269
270#[cfg(test)]
271mod tests {
272    use super::*;
273    use u_numflow::distributions::Sample;
274    use u_numflow::distributions::{LogNormal, Normal};
275
276    /// Generate uniform-like data for testing.
277    fn uniform_data(n: usize) -> Vec<f64> {
278        (0..n)
279            .map(|i| 10.0 + (i as f64) / (n as f64 - 1.0) * 10.0)
280            .collect()
281    }
282
283    fn pc(
284        data: &[f64],
285        lsl: Option<f64>,
286        usl: Option<f64>,
287    ) -> Result<PercentileCapabilityResult, PercentileCapabilityError> {
288        percentile_capability(data, lsl, usl, CapabilityModel::Auto)
289    }
290
291    #[test]
292    fn basic_two_sided() {
293        let data = uniform_data(100);
294        let result = pc(&data, Some(5.0), Some(25.0)).unwrap();
295
296        assert!(result.cp_star.is_some());
297        assert!(result.cpk_star.is_some());
298        assert!(result.cpu_star.is_some());
299        assert!(result.cpl_star.is_some());
300        assert!(result.cp_star.unwrap() > 0.0);
301        assert!(
302            (result.median - 15.0).abs() < 0.5,
303            "median should be ~15, got {}",
304            result.median
305        );
306    }
307
308    #[test]
309    fn usl_only() {
310        let data = uniform_data(50);
311        let result = pc(&data, None, Some(25.0)).unwrap();
312
313        assert!(result.cp_star.is_none(), "Cp* requires both limits");
314        assert!(result.cpu_star.is_some());
315        assert!(result.cpl_star.is_none());
316        assert_eq!(result.cpk_star, result.cpu_star);
317    }
318
319    #[test]
320    fn lsl_only() {
321        let data = uniform_data(50);
322        let result = pc(&data, Some(5.0), None).unwrap();
323
324        assert!(result.cp_star.is_none());
325        assert!(result.cpu_star.is_none());
326        assert!(result.cpl_star.is_some());
327        assert_eq!(result.cpk_star, result.cpl_star);
328    }
329
330    #[test]
331    fn refusals_name_their_reason() {
332        let few: Vec<f64> = (0..19).map(|i| i as f64).collect();
333        assert_eq!(
334            pc(&few, Some(0.0), Some(20.0)).unwrap_err(),
335            PercentileCapabilityError::InsufficientData { got: 19 }
336        );
337        let data = uniform_data(30);
338        assert_eq!(
339            pc(&data, None, None).unwrap_err(),
340            PercentileCapabilityError::NoSpecification
341        );
342        assert_eq!(
343            pc(&data, Some(10.0), Some(5.0)).unwrap_err(),
344            PercentileCapabilityError::InvalidSpecification
345        );
346        assert_eq!(
347            pc(&data, Some(10.0), Some(10.0)).unwrap_err(),
348            PercentileCapabilityError::InvalidSpecification
349        );
350        let mut nan = uniform_data(30);
351        nan[10] = f64::NAN;
352        assert_eq!(
353            pc(&nan, Some(0.0), Some(30.0)).unwrap_err(),
354            PercentileCapabilityError::NonFiniteData { index: 10 }
355        );
356        assert_eq!(
357            pc(&[3.0; 25], Some(0.0), Some(30.0)).unwrap_err(),
358            PercentileCapabilityError::NoVariation
359        );
360        let mut negative = uniform_data(30);
361        negative[0] = -1.0;
362        assert_eq!(
363            percentile_capability(
364                &negative,
365                Some(-5.0),
366                Some(30.0),
367                CapabilityModel::LogNormal
368            )
369            .unwrap_err(),
370            PercentileCapabilityError::ModelNotFitted(CapabilityModel::LogNormal)
371        );
372    }
373
374    #[test]
375    fn cpk_star_equals_min_of_cpu_cpl() {
376        let data = uniform_data(100);
377        let result = pc(&data, Some(5.0), Some(25.0)).unwrap();
378        let (cpu, cpl, cpk) = (
379            result.cpu_star.unwrap(),
380            result.cpl_star.unwrap(),
381            result.cpk_star.unwrap(),
382        );
383        assert!((cpk - cpu.min(cpl)).abs() < 1e-10);
384    }
385
386    #[test]
387    fn symmetric_data_gives_balanced_cpu_cpl() {
388        let data: Vec<f64> = (0..200).map(|i| 40.0 + (i as f64) * 0.1).collect();
389        let result = pc(&data, Some(30.0), Some(70.0)).unwrap();
390        let (cpu, cpl) = (result.cpu_star.unwrap(), result.cpl_star.unwrap());
391        assert!(
392            (cpu - cpl).abs() < cpu * 0.2,
393            "symmetric data should give similar Cpu* ({cpu}) and Cpl* ({cpl})"
394        );
395    }
396
397    /// The defect this module was rewritten for: on normal data with a true Cp
398    /// of 1.0 (limits at ±3σ), the empirical 0.135 % / 99.865 % points are the
399    /// sample minimum and maximum, so Cp* averaged 1.37 at n = 50 and 1.68 at
400    /// n = 20. Quantiles of the fitted normal stay near 1.0 at every size.
401    #[test]
402    fn small_samples_do_not_inflate_cp_star() {
403        let process = Normal::new(10.0, 1.0).unwrap();
404        for (n, empirical_at_least) in [(20_usize, 1.5), (50, 1.25), (100, 1.15)] {
405            let reps = 400;
406            let (mut sum, mut empirical) = (0.0, 0.0);
407            for rep in 0..reps {
408                let data = process.sample_n(n, 1_000 + rep as u64 * 7 + n as u64);
409                let r =
410                    percentile_capability(&data, Some(7.0), Some(13.0), CapabilityModel::Normal)
411                        .unwrap();
412                sum += r.cp_star.unwrap();
413                // What the empirical points gave on the same samples.
414                let lo = u_numflow::stats::quantile(&data, P_LOWER).unwrap();
415                let hi = u_numflow::stats::quantile(&data, P_UPPER).unwrap();
416                empirical += 6.0 / (hi - lo);
417            }
418            let mean = sum / reps as f64;
419            let empirical = empirical / reps as f64;
420            assert!(
421                empirical > empirical_at_least,
422                "n = {n}: the empirical Cp* {empirical} no longer shows the bias this guards"
423            );
424            // The MLE σ is low by about (n−1)/n, and 1/σ̂ is biased up; both
425            // are a few percent at these sizes. The empirical method was 37 %
426            // high at n = 50.
427            assert!(
428                (mean - 1.0).abs() < 0.08,
429                "n = {n}: mean Cp* {mean} (true 1.0)"
430            );
431        }
432    }
433
434    #[test]
435    fn quantiles_are_the_fitted_model_s() {
436        let process = Normal::new(50.0, 2.0).unwrap();
437        let data = process.sample_n(200, 42);
438        let r =
439            percentile_capability(&data, Some(40.0), Some(60.0), CapabilityModel::Normal).unwrap();
440        let mu = r.fit.parameters.iter().find(|(n, _)| n == "mu").unwrap().1;
441        let sigma = r
442            .fit
443            .parameters
444            .iter()
445            .find(|(n, _)| n == "sigma")
446            .unwrap()
447            .1;
448        assert_eq!(r.model, CapabilityModel::Normal);
449        assert!((r.median - mu).abs() < 1e-9);
450        assert!((r.percentile_lower - (mu - 3.0 * sigma)).abs() < 1e-6);
451        assert!((r.percentile_upper - (mu + 3.0 * sigma)).abs() < 1e-6);
452    }
453
454    #[test]
455    fn auto_reports_the_family_it_used() {
456        let process = LogNormal::new(1.0, 0.6).unwrap();
457        let data = process.sample_n(300, 7);
458        let r = pc(&data, None, Some(30.0)).unwrap();
459        assert_ne!(r.model, CapabilityModel::Auto);
460        assert_eq!(
461            CapabilityModel::of_fit(&fit_best(&data)[0].distribution),
462            Some(r.model)
463        );
464    }
465}
466
467/// Against SciPy 1.17.1: the family fitted by maximum likelihood
468/// (`lognorm.fit(x, floc=0)`; the Weibull's from its likelihood equation), its
469/// `ppf` at `norm.cdf(-3)`, 0.5 and `norm.cdf(3)` -- the equivalents of
470/// ±3 sigma that ISO 22514-2 rounds to 0.135 % and 99.865 %, which this
471/// function takes exactly -- and the ISO 22514-2 ratios taken from them.
472#[cfg(test)]
473mod scipy_reference {
474    use super::*;
475
476    const DATA: [f64; 40] = [
477        3.195, 5.365, 1.217, 5.806, 2.661, 1.896, 1.894, 1.67, 2.458, 3.956, 3.535, 3.623, 1.268,
478        1.341, 5.47, 4.204, 7.26, 4.685, 1.714, 4.847, 3.607, 2.994, 1.88, 2.722, 2.544, 4.062,
479        2.875, 5.546, 1.918, 3.36, 3.402, 1.842, 1.597, 4.289, 3.431, 2.527, 3.246, 2.006, 3.026,
480        2.07,
481    ];
482
483    fn close(got: f64, want: f64, rel: f64, what: &str) {
484        assert!(
485            (got - want).abs() <= rel * want.abs(),
486            "{what}: got {got}, SciPy {want}"
487        );
488    }
489
490    fn check(model: CapabilityModel, rel: f64, want: [f64; 7]) {
491        let r = percentile_capability(&DATA, Some(0.9), Some(9.0), model).expect("valid");
492        assert_eq!(r.model, model);
493        let [lo, med, hi, cp, cpk, cpu, cpl] = want;
494        close(r.percentile_lower, lo, rel, "X0.135");
495        close(r.median, med, rel, "X50");
496        close(r.percentile_upper, hi, rel, "X99.865");
497        close(r.cp_star.expect("both limits"), cp, rel, "Cp*");
498        close(r.cpk_star.expect("limits"), cpk, rel, "Cpk*");
499        close(r.cpu_star.expect("USL"), cpu, rel, "Cpu*");
500        close(r.cpl_star.expect("LSL"), cpl, rel, "Cpl*");
501    }
502
503    #[test]
504    fn log_normal_percentile_capability_reproduces_scipy() {
505        // The log-normal MLE is closed form: tight.
506        check(
507            CapabilityModel::LogNormal,
508            1e-9,
509            [
510                0.7542968892328896,
511                2.8798090991738787,
512                10.994743006455773,
513                0.7909811650077455,
514                0.754188631817962,
515                0.754188631817962,
516                0.9314503534321473,
517            ],
518        );
519    }
520
521    #[test]
522    fn weibull_percentile_capability_reproduces_scipy() {
523        // The exact MLE: the shape's likelihood equation solved to 1e-15
524        // (`scipy.optimize.brentq`), the scale from it. `weibull_min.fit` stops
525        // short of it -- its shape is 2.413418 against 2.413394 -- so its
526        // quantiles are not the reference.
527        check(
528            CapabilityModel::Weibull,
529            1e-9,
530            [
531                0.23258289429884546,
532                3.0872615285886305,
533                7.858145308508201,
534                1.0622167336676158,
535                0.7662023676906101,
536                1.2393381906090049,
537                0.7662023676906101,
538            ],
539        );
540    }
541}
u_analytics/capability/percentile.rs

u_analytics/capability/
percentile.rs
