fdars_core/spm/
frcc.rs

1//! Functional Regression Control Chart (FRCC).
2//!
3//! Monitors a functional response after adjusting for known scalar covariates
4//! using function-on-scalar regression (FOSR). The residuals are then monitored
5//! via FPCA-based T-squared and SPE statistics.
6//!
7//! This is useful when the process output (functional) depends on known inputs
8//! (scalar predictors), and we want to detect deviations beyond what the inputs
9//! explain.
10
11use crate::error::FdarError;
12use crate::function_on_scalar::{fosr, predict_fosr, FosrResult};
13use crate::matrix::FdMatrix;
14use crate::regression::{fdata_to_pc_1d, FpcaResult};
15
16use super::control::{spe_control_limit, t2_control_limit, ControlLimit};
17use super::stats::{hotelling_t2, spe_univariate};
18
19/// Configuration for FRCC chart construction.
20#[derive(Debug, Clone, PartialEq)]
21pub struct FrccConfig {
22    /// Number of principal components for residual FPCA (default 5).
23    pub ncomp: usize,
24    /// FOSR smoothing parameter; use small positive value (default 1e-4).
25    pub fosr_lambda: f64,
26    /// Significance level (default 0.05).
27    pub alpha: f64,
28    /// Fraction of data for tuning (default 0.5).
29    pub tuning_fraction: f64,
30    /// Random seed (default 42).
31    pub seed: u64,
32}
33
34impl Default for FrccConfig {
35    fn default() -> Self {
36        Self {
37            ncomp: 5,
38            fosr_lambda: 1e-4,
39            alpha: 0.05,
40            tuning_fraction: 0.5,
41            seed: 42,
42        }
43    }
44}
45
46/// Phase I FRCC chart.
47#[derive(Debug, Clone, PartialEq)]
48#[non_exhaustive]
49pub struct FrccChart {
50    /// Fitted FOSR model.
51    pub fosr: FosrResult,
52    /// FPCA on calibration residuals.
53    pub residual_fpca: FpcaResult,
54    /// Eigenvalues from residual FPCA.
55    pub eigenvalues: Vec<f64>,
56    /// T-squared control limit.
57    pub t2_limit: ControlLimit,
58    /// SPE control limit.
59    pub spe_limit: ControlLimit,
60    /// Configuration used.
61    pub config: FrccConfig,
62}
63
64/// Result of FRCC monitoring.
65#[derive(Debug, Clone, PartialEq)]
66#[non_exhaustive]
67pub struct FrccMonitorResult {
68    /// T-squared values.
69    pub t2: Vec<f64>,
70    /// SPE values.
71    pub spe: Vec<f64>,
72    /// T-squared alarm flags.
73    pub t2_alarm: Vec<bool>,
74    /// SPE alarm flags.
75    pub spe_alarm: Vec<bool>,
76    /// Residual scores.
77    pub residual_scores: FdMatrix,
78}
79
80/// Split indices into tuning and calibration sets (same logic as phase.rs).
81fn split_indices(n: usize, tuning_fraction: f64, seed: u64) -> (Vec<usize>, Vec<usize>) {
82    let n_tune = ((n as f64 * tuning_fraction).round() as usize)
83        .max(2)
84        .min(n - 1);
85
86    let mut indices: Vec<usize> = (0..n).collect();
87    let mut rng_state: u64 = seed;
88    for i in (1..n).rev() {
89        rng_state = rng_state
90            .wrapping_mul(6_364_136_223_846_793_005)
91            .wrapping_add(1);
92        let j = (rng_state >> 33) as usize % (i + 1);
93        indices.swap(i, j);
94    }
95
96    let tune_indices: Vec<usize> = indices[..n_tune].to_vec();
97    let cal_indices: Vec<usize> = indices[n_tune..].to_vec();
98    (tune_indices, cal_indices)
99}
100
101/// Compute residuals: observed - predicted.
102fn compute_residuals(observed: &FdMatrix, predicted: &FdMatrix) -> FdMatrix {
103    let (n, m) = observed.shape();
104    let mut residuals = FdMatrix::zeros(n, m);
105    for i in 0..n {
106        for j in 0..m {
107            residuals[(i, j)] = observed[(i, j)] - predicted[(i, j)];
108        }
109    }
110    residuals
111}
112
113/// Compute centered reconstruction (without adding back mean) for SPE.
114fn centered_reconstruct(fpca: &FpcaResult, scores: &FdMatrix, ncomp: usize) -> FdMatrix {
115    let n = scores.nrows();
116    let m = fpca.mean.len();
117    let ncomp = ncomp.min(fpca.rotation.ncols()).min(scores.ncols());
118
119    let mut recon = FdMatrix::zeros(n, m);
120    for i in 0..n {
121        for j in 0..m {
122            let mut val = 0.0;
123            for k in 0..ncomp {
124                val += scores[(i, k)] * fpca.rotation[(j, k)];
125            }
126            recon[(i, j)] = val;
127        }
128    }
129    recon
130}
131
132/// Center data by subtracting mean.
133fn center_data(data: &FdMatrix, mean: &[f64]) -> FdMatrix {
134    let (n, m) = data.shape();
135    let mut centered = FdMatrix::zeros(n, m);
136    for i in 0..n {
137        for j in 0..m {
138            centered[(i, j)] = data[(i, j)] - mean[j];
139        }
140    }
141    centered
142}
143
144/// Build a Functional Regression Control Chart from Phase I data.
145///
146/// 1. Splits data into tuning and calibration sets
147/// 2. Fits FOSR on tuning set
148/// 3. Computes residuals on calibration set
149/// 4. Runs FPCA on calibration residuals
150/// 5. Computes T-squared and SPE control limits
151///
152/// # Arguments
153/// * `y_curves` - Functional response (n x m)
154/// * `predictors` - Scalar predictors (n x p)
155/// * `argvals` - Grid points (length m)
156/// * `config` - FRCC configuration
157///
158/// # Errors
159///
160/// Returns errors from FOSR fitting, FPCA, or control limit estimation.
161#[must_use = "expensive computation whose result should not be discarded"]
162pub fn frcc_phase1(
163    y_curves: &FdMatrix,
164    predictors: &FdMatrix,
165    argvals: &[f64],
166    config: &FrccConfig,
167) -> Result<FrccChart, FdarError> {
168    let (n, m) = y_curves.shape();
169    let p = predictors.ncols();
170
171    if n < 6 {
172        return Err(FdarError::InvalidDimension {
173            parameter: "y_curves",
174            expected: "at least 6 observations".to_string(),
175            actual: format!("{n} observations"),
176        });
177    }
178    if predictors.nrows() != n {
179        return Err(FdarError::InvalidDimension {
180            parameter: "predictors",
181            expected: format!("{n} rows"),
182            actual: format!("{} rows", predictors.nrows()),
183        });
184    }
185    if argvals.len() != m {
186        return Err(FdarError::InvalidDimension {
187            parameter: "argvals",
188            expected: format!("{m}"),
189            actual: format!("{}", argvals.len()),
190        });
191    }
192
193    // Split
194    let (tune_idx, cal_idx) = split_indices(n, config.tuning_fraction, config.seed);
195
196    let tune_y = crate::cv::subset_rows(y_curves, &tune_idx);
197    let tune_x = crate::cv::subset_rows(predictors, &tune_idx);
198    let cal_y = crate::cv::subset_rows(y_curves, &cal_idx);
199    let cal_x = crate::cv::subset_rows(predictors, &cal_idx);
200
201    let n_tune = tune_y.nrows();
202    let n_cal = cal_y.nrows();
203
204    // Ensure enough observations for FOSR (needs n >= p + 2)
205    if n_tune < p + 2 {
206        return Err(FdarError::InvalidDimension {
207            parameter: "y_curves",
208            expected: format!("tuning set with at least {} observations", p + 2),
209            actual: format!("{n_tune} observations in tuning set"),
210        });
211    }
212
213    // FOSR on tuning set
214    let fosr_lambda = if config.fosr_lambda < 0.0 {
215        1e-4
216    } else {
217        config.fosr_lambda
218    };
219    let fosr_result = fosr(&tune_y, &tune_x, fosr_lambda)?;
220
221    // Predict on calibration set and compute residuals
222    let cal_predicted = predict_fosr(&fosr_result, &cal_x);
223    let cal_residuals = compute_residuals(&cal_y, &cal_predicted);
224
225    // FPCA on calibration residuals
226    let ncomp = config.ncomp.min(n_cal - 1).min(m);
227    let residual_fpca = fdata_to_pc_1d(&cal_residuals, ncomp)?;
228    let actual_ncomp = residual_fpca.scores.ncols();
229
230    // Eigenvalues
231    let eigenvalues: Vec<f64> = residual_fpca
232        .singular_values
233        .iter()
234        .take(actual_ncomp)
235        .map(|&sv| sv * sv / (n_cal as f64 - 1.0))
236        .collect();
237
238    // T-squared on calibration residual scores
239    let _t2_cal = hotelling_t2(&residual_fpca.scores, &eigenvalues)?;
240
241    // SPE on calibration residuals
242    let cal_resid_centered = center_data(&cal_residuals, &residual_fpca.mean);
243    let cal_resid_recon = centered_reconstruct(&residual_fpca, &residual_fpca.scores, actual_ncomp);
244    let spe_cal = spe_univariate(&cal_resid_centered, &cal_resid_recon, argvals)?;
245
246    // Control limits
247    let t2_limit = t2_control_limit(actual_ncomp, config.alpha)?;
248    let spe_limit = spe_control_limit(&spe_cal, config.alpha)?;
249
250    Ok(FrccChart {
251        fosr: fosr_result,
252        residual_fpca,
253        eigenvalues,
254        t2_limit,
255        spe_limit,
256        config: config.clone(),
257    })
258}
259
260/// Monitor new data against a Functional Regression Control Chart.
261///
262/// 1. Predicts functional response from FOSR model
263/// 2. Computes residuals
264/// 3. Projects residuals through FPCA
265/// 4. Computes T-squared and SPE
266///
267/// # Arguments
268/// * `chart` - Phase I FRCC chart
269/// * `new_y` - New functional response (n_new x m)
270/// * `new_predictors` - New scalar predictors (n_new x p)
271/// * `argvals` - Grid points (length m)
272///
273/// # Errors
274///
275/// Returns errors from FOSR prediction, FPCA projection, or statistic computation.
276#[must_use = "monitoring result should not be discarded"]
277pub fn frcc_monitor(
278    chart: &FrccChart,
279    new_y: &FdMatrix,
280    new_predictors: &FdMatrix,
281    argvals: &[f64],
282) -> Result<FrccMonitorResult, FdarError> {
283    let m = chart.residual_fpca.mean.len();
284    if new_y.ncols() != m {
285        return Err(FdarError::InvalidDimension {
286            parameter: "new_y",
287            expected: format!("{m} columns"),
288            actual: format!("{} columns", new_y.ncols()),
289        });
290    }
291    if new_y.nrows() != new_predictors.nrows() {
292        return Err(FdarError::InvalidDimension {
293            parameter: "new_predictors",
294            expected: format!("{} rows", new_y.nrows()),
295            actual: format!("{} rows", new_predictors.nrows()),
296        });
297    }
298
299    let ncomp = chart.eigenvalues.len();
300
301    // Predict and compute residuals
302    let predicted = predict_fosr(&chart.fosr, new_predictors);
303    let residuals = compute_residuals(new_y, &predicted);
304
305    // Project residuals through FPCA
306    let residual_scores = chart.residual_fpca.project(&residuals)?;
307
308    // T-squared
309    let t2 = hotelling_t2(&residual_scores, &chart.eigenvalues)?;
310
311    // SPE
312    let resid_centered = center_data(&residuals, &chart.residual_fpca.mean);
313    let resid_recon = centered_reconstruct(&chart.residual_fpca, &residual_scores, ncomp);
314    let spe = spe_univariate(&resid_centered, &resid_recon, argvals)?;
315
316    // Alarms
317    let t2_alarm: Vec<bool> = t2.iter().map(|&v| v > chart.t2_limit.ucl).collect();
318    let spe_alarm: Vec<bool> = spe.iter().map(|&v| v > chart.spe_limit.ucl).collect();
319
320    Ok(FrccMonitorResult {
321        t2,
322        spe,
323        t2_alarm,
324        spe_alarm,
325        residual_scores,
326    })
327}
fdars_core/spm/frcc.rs

fdars_core/spm/
frcc.rs
