fdars_core/alignment/
shape_ci.rs

1//! Bootstrap confidence intervals for curve shapes in the elastic metric.
2
3use rand::Rng;
4
5use super::karcher::karcher_mean;
6use super::pairwise::elastic_align_pair;
7use crate::error::FdarError;
8use crate::iter_maybe_parallel;
9use crate::matrix::FdMatrix;
10#[cfg(feature = "parallel")]
11use rayon::iter::ParallelIterator;
12
13// ─── Types ──────────────────────────────────────────────────────────────────
14
15/// Configuration for shape bootstrap confidence intervals.
16///
17/// Construct via `ShapeCiConfig::default()`, then assign the fields you need (e.g. `let mut c = ShapeCiConfig::default(); c.field = …;`). This struct is `#[non_exhaustive]`, so external crates cannot build it with a struct literal — not even functional-update `..Default::default()` form.
18#[non_exhaustive]
19#[derive(Debug, Clone, PartialEq)]
20pub struct ShapeCiConfig {
21    /// Number of bootstrap resamples.
22    pub n_bootstrap: usize,
23    /// Confidence level (e.g., 0.95 for 95% CI).
24    pub confidence_level: f64,
25    /// Roughness penalty for elastic alignment.
26    pub lambda: f64,
27    /// Maximum Karcher mean iterations.
28    pub max_iter: usize,
29    /// Convergence tolerance for the Karcher mean.
30    pub tol: f64,
31    /// Random seed for reproducibility.
32    pub seed: u64,
33}
34
35impl Default for ShapeCiConfig {
36    fn default() -> Self {
37        Self {
38            n_bootstrap: 200,
39            confidence_level: 0.95,
40            lambda: 0.0,
41            max_iter: 15,
42            tol: 1e-3,
43            seed: 42,
44        }
45    }
46}
47
48/// Result of shape bootstrap confidence interval computation.
49#[derive(Debug, Clone, PartialEq)]
50#[non_exhaustive]
51pub struct ShapeCiResult {
52    /// Karcher mean of the full sample.
53    pub mean: Vec<f64>,
54    /// Pointwise lower confidence band (length m).
55    pub lower_band: Vec<f64>,
56    /// Pointwise upper confidence band (length m).
57    pub upper_band: Vec<f64>,
58    /// Bootstrap Karcher means (n_bootstrap x m).
59    pub bootstrap_means: FdMatrix,
60}
61
62// ─── Public API ─────────────────────────────────────────────────────────────
63
64/// Compute bootstrap confidence intervals for the elastic Karcher mean.
65///
66/// Resamples the input curves with replacement, computes the Karcher mean
67/// of each bootstrap sample, aligns each bootstrap mean to the full-sample
68/// mean, and derives pointwise confidence bands from the empirical quantiles.
69///
70/// # Arguments
71/// * `data`    - Functional data matrix (n x m).
72/// * `argvals` - Evaluation points (length m).
73/// * `config`  - Bootstrap configuration.
74///
75/// # Errors
76/// Returns [`FdarError::InvalidDimension`] if `n < 3` or `argvals` length
77/// does not match `m`.
78/// Returns [`FdarError::InvalidParameter`] if `confidence_level` is not in
79/// `(0, 1)` or `n_bootstrap < 1`.
80#[must_use = "expensive computation whose result should not be discarded"]
81pub fn shape_confidence_interval(
82    data: &FdMatrix,
83    argvals: &[f64],
84    config: &ShapeCiConfig,
85) -> Result<ShapeCiResult, FdarError> {
86    let (n, m) = data.shape();
87
88    // ── Validation ──
89    if argvals.len() != m {
90        return Err(FdarError::InvalidDimension {
91            parameter: "argvals",
92            expected: format!("{m}"),
93            actual: format!("{}", argvals.len()),
94        });
95    }
96    if n < 3 {
97        return Err(FdarError::InvalidDimension {
98            parameter: "data",
99            expected: "at least 3 rows".to_string(),
100            actual: format!("{n} rows"),
101        });
102    }
103    if config.confidence_level <= 0.0 || config.confidence_level >= 1.0 {
104        return Err(FdarError::InvalidParameter {
105            parameter: "confidence_level",
106            message: format!("must be in (0, 1), got {}", config.confidence_level),
107        });
108    }
109    if config.n_bootstrap < 1 {
110        return Err(FdarError::InvalidParameter {
111            parameter: "n_bootstrap",
112            message: format!("must be >= 1, got {}", config.n_bootstrap),
113        });
114    }
115
116    // ── Full-sample Karcher mean ──
117    let full_karcher = karcher_mean(data, argvals, config.max_iter, config.tol, config.lambda);
118
119    // ── Bootstrap loop ──
120    let boot_means: Vec<Vec<f64>> = iter_maybe_parallel!(0..config.n_bootstrap)
121        .map(|b| {
122            let mut rng = crate::helpers::seed_for_thread(config.seed, b);
123
124            // Resample n indices with replacement
125            let indices: Vec<usize> = (0..n).map(|_| rng.gen_range(0..n)).collect();
126
127            // Build bootstrap matrix
128            let mut boot_data = FdMatrix::zeros(n, m);
129            for (row, &idx) in indices.iter().enumerate() {
130                for j in 0..m {
131                    boot_data[(row, j)] = data[(idx, j)];
132                }
133            }
134
135            // Compute bootstrap Karcher mean
136            let boot_karcher = karcher_mean(
137                &boot_data,
138                argvals,
139                config.max_iter,
140                config.tol,
141                config.lambda,
142            );
143
144            // Align bootstrap mean to full-sample mean
145            let aligned = elastic_align_pair(
146                &full_karcher.mean,
147                &boot_karcher.mean,
148                argvals,
149                config.lambda,
150            );
151
152            aligned.f_aligned
153        })
154        .collect();
155
156    // ── Build bootstrap_means matrix ──
157    let mut bootstrap_means = FdMatrix::zeros(config.n_bootstrap, m);
158    for (b, bm) in boot_means.iter().enumerate() {
159        for j in 0..m {
160            bootstrap_means[(b, j)] = bm[j];
161        }
162    }
163
164    // ── Pointwise confidence bands ──
165    let alpha = 1.0 - config.confidence_level;
166    let mut lower_band = vec![0.0; m];
167    let mut upper_band = vec![0.0; m];
168
169    for j in 0..m {
170        let mut col_vals: Vec<f64> = (0..config.n_bootstrap)
171            .map(|b| bootstrap_means[(b, j)])
172            .collect();
173        col_vals.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
174
175        lower_band[j] = quantile_sorted(&col_vals, alpha / 2.0);
176        upper_band[j] = quantile_sorted(&col_vals, 1.0 - alpha / 2.0);
177    }
178
179    Ok(ShapeCiResult {
180        mean: full_karcher.mean,
181        lower_band,
182        upper_band,
183        bootstrap_means,
184    })
185}
186
187use crate::helpers::quantile_sorted;
188
189// ─── Tests ──────────────────────────────────────────────────────────────────
190
191#[cfg(test)]
192mod tests {
193    use super::*;
194    use crate::simulation::{sim_fundata, EFunType, EValType};
195    use crate::test_helpers::uniform_grid;
196
197    fn make_data(n: usize, m: usize) -> (FdMatrix, Vec<f64>) {
198        let t = uniform_grid(m);
199        let data = sim_fundata(n, &t, 3, EFunType::Fourier, EValType::Exponential, Some(99));
200        (data, t)
201    }
202
203    #[test]
204    fn shape_ci_band_contains_mean() {
205        let (data, t) = make_data(8, 20);
206        let config = ShapeCiConfig {
207            n_bootstrap: 30,
208            confidence_level: 0.95,
209            max_iter: 5,
210            tol: 1e-2,
211            ..Default::default()
212        };
213        let result = shape_confidence_interval(&data, &t, &config).unwrap();
214        let m = t.len();
215        for j in 0..m {
216            assert!(
217                result.lower_band[j] <= result.mean[j] + 1e-6
218                    && result.mean[j] <= result.upper_band[j] + 1e-6,
219                "mean[{j}]={} not in [{}, {}]",
220                result.mean[j],
221                result.lower_band[j],
222                result.upper_band[j],
223            );
224        }
225    }
226
227    #[test]
228    fn shape_ci_band_width_positive() {
229        let (data, t) = make_data(8, 20);
230        let config = ShapeCiConfig {
231            n_bootstrap: 30,
232            confidence_level: 0.95,
233            max_iter: 5,
234            tol: 1e-2,
235            ..Default::default()
236        };
237        let result = shape_confidence_interval(&data, &t, &config).unwrap();
238        let m = t.len();
239        let n_positive = (0..m)
240            .filter(|&j| result.upper_band[j] > result.lower_band[j] + 1e-12)
241            .count();
242        assert!(
243            n_positive > m / 2,
244            "upper > lower for only {n_positive}/{m} points, expected > {}/{}",
245            m / 2,
246            m
247        );
248    }
249
250    #[test]
251    fn shape_ci_bootstrap_means_shape() {
252        let (data, t) = make_data(6, 20);
253        let n_boot = 15;
254        let config = ShapeCiConfig {
255            n_bootstrap: n_boot,
256            confidence_level: 0.90,
257            max_iter: 3,
258            tol: 1e-2,
259            ..Default::default()
260        };
261        let result = shape_confidence_interval(&data, &t, &config).unwrap();
262        assert_eq!(result.bootstrap_means.shape(), (n_boot, t.len()));
263    }
264
265    #[test]
266    fn shape_ci_rejects_too_few_curves() {
267        let t = uniform_grid(20);
268        let data = FdMatrix::zeros(2, 20);
269        let config = ShapeCiConfig::default();
270        assert!(shape_confidence_interval(&data, &t, &config).is_err());
271    }
272}
fdars_core/alignment/shape_ci.rs

fdars_core/alignment/
shape_ci.rs
