Skip to main content

fdars_core/spm/
control.rs

1//! Control limits for SPM monitoring statistics.
2//!
3//! Provides upper control limits (UCL) for T-squared and SPE statistics:
4//! - T-squared: chi-squared distribution quantile
5//! - SPE: moment-matched chi-squared approximation (Box, 1954)
6
7use super::chi_squared::chi2_quantile;
8use crate::error::FdarError;
9
10/// A control limit for a monitoring statistic.
11#[derive(Debug, Clone, PartialEq)]
12#[non_exhaustive]
13pub struct ControlLimit {
14    /// Upper control limit value.
15    pub ucl: f64,
16    /// Significance level used.
17    pub alpha: f64,
18    /// Human-readable description of how the limit was computed.
19    pub description: String,
20}
21
22/// Compute the T-squared control limit based on the chi-squared distribution.
23///
24/// UCL = chi2_quantile(1 - alpha, ncomp)
25///
26/// # Arguments
27/// * `ncomp` - Number of principal components (degrees of freedom)
28/// * `alpha` - Significance level (e.g. 0.05)
29///
30/// # Errors
31///
32/// Returns [`FdarError::InvalidParameter`] if `ncomp` is 0 or `alpha` is not in (0, 1).
33pub fn t2_control_limit(ncomp: usize, alpha: f64) -> Result<ControlLimit, FdarError> {
34    if ncomp == 0 {
35        return Err(FdarError::InvalidParameter {
36            parameter: "ncomp",
37            message: "ncomp must be at least 1".to_string(),
38        });
39    }
40    if alpha <= 0.0 || alpha >= 1.0 {
41        return Err(FdarError::InvalidParameter {
42            parameter: "alpha",
43            message: format!("alpha must be in (0, 1), got {alpha}"),
44        });
45    }
46
47    let ucl = chi2_quantile(1.0 - alpha, ncomp);
48
49    Ok(ControlLimit {
50        ucl,
51        alpha,
52        description: format!("T2 ~ chi2({ncomp}), alpha={alpha}"),
53    })
54}
55
56/// Compute the SPE control limit using moment-matched chi-squared approximation.
57///
58/// Estimates parameters a and b such that SPE ~ a * chi2(b):
59///   a = var(spe) / (2 * mean(spe))
60///   b = 2 * mean(spe)^2 / var(spe)
61///   UCL = a * chi2_quantile(1 - alpha, round(b))
62///
63/// # Arguments
64/// * `spe_values` - In-control SPE values from calibration data
65/// * `alpha` - Significance level (e.g. 0.05)
66///
67/// # Errors
68///
69/// Returns [`FdarError::InvalidDimension`] if `spe_values` has fewer than 2 elements.
70/// Returns [`FdarError::InvalidParameter`] if `alpha` is not in (0, 1).
71/// Returns [`FdarError::ComputationFailed`] if the estimated variance is zero or negative.
72pub fn spe_control_limit(spe_values: &[f64], alpha: f64) -> Result<ControlLimit, FdarError> {
73    let n = spe_values.len();
74    if n < 2 {
75        return Err(FdarError::InvalidDimension {
76            parameter: "spe_values",
77            expected: "at least 2 values".to_string(),
78            actual: format!("{n} values"),
79        });
80    }
81    if alpha <= 0.0 || alpha >= 1.0 {
82        return Err(FdarError::InvalidParameter {
83            parameter: "alpha",
84            message: format!("alpha must be in (0, 1), got {alpha}"),
85        });
86    }
87
88    let mean: f64 = spe_values.iter().sum::<f64>() / n as f64;
89    let var: f64 = spe_values.iter().map(|&s| (s - mean).powi(2)).sum::<f64>() / (n as f64 - 1.0);
90
91    if mean <= 0.0 || var <= 0.0 {
92        return Err(FdarError::ComputationFailed {
93            operation: "spe_control_limit",
94            detail: format!(
95                "SPE mean={mean:.6}, var={var:.6}; cannot estimate chi-squared parameters from non-positive moments"
96            ),
97        });
98    }
99
100    let a = var / (2.0 * mean);
101    let b = 2.0 * mean * mean / var;
102    let b_int = b.round().max(1.0) as usize;
103
104    let ucl = a * chi2_quantile(1.0 - alpha, b_int);
105
106    Ok(ControlLimit {
107        ucl,
108        alpha,
109        description: format!("SPE ~ {a:.4} * chi2({b_int}), alpha={alpha}"),
110    })
111}