stats_claw/resampling/
jackknife.rs

1//! Deterministic jackknife (leave-one-out) resampling.
2//!
3//! The jackknife recomputes a statistic on each of the `n` subsamples that omit
4//! one observation, then combines the resulting replicates into an estimate of
5//! the statistic's bias and standard error. Unlike the bootstrap it draws no
6//! random numbers, so [`jackknife_statistic`] is fully deterministic for a given
7//! input.
8//!
9//! # Examples
10//!
11//! ```
12//! use stats_claw::resampling::jackknife_indices;
13//!
14//! // Three leave-one-out index sets, each omitting one position in order.
15//! assert_eq!(jackknife_indices(3)?, vec![vec![1, 2], vec![0, 2], vec![0, 1]]);
16//! # Ok::<(), stats_claw::error::Error>(())
17//! ```
18
19use crate::error::{Error, Result};
20use crate::numeric::count_to_f64;
21use crate::resampling::JackknifeResampling;
22
23/// The jackknife bias and standard-error estimate for a statistic.
24///
25/// Bundles the statistic evaluated on the full sample together with the
26/// leave-one-out replicates and the two classical jackknife summaries computed
27/// from them (Efron & Tibshirani, *An Introduction to the Bootstrap*, 1993,
28/// §10.2). The fields are private (the struct owns a `Vec`, so it fully
29/// encapsulates its storage); read them through the [`estimate`](Self::estimate),
30/// [`bias`](Self::bias), [`std_error`](Self::std_error), and
31/// [`replicates`](Self::replicates) accessors.
32///
33/// # Examples
34///
35/// ```
36/// use stats_claw::resampling::jackknife_statistic;
37///
38/// let data = [2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0];
39/// let mean = |s: &[f64]| s.iter().sum::<f64>() / f64::from(u32::try_from(s.len()).unwrap_or(0));
40/// let est = jackknife_statistic(&data, mean)?;
41/// assert!((est.estimate() - 5.0).abs() < 1e-12, "estimate was {}", est.estimate());
42/// # Ok::<(), stats_claw::error::Error>(())
43/// ```
44#[derive(Debug, Clone, PartialEq)]
45pub struct JackknifeEstimate {
46    /// The statistic evaluated on the full sample.
47    estimate: f64,
48    /// Jackknife bias estimate, `(n - 1) * (mean(replicates) - estimate)`.
49    bias: f64,
50    /// Jackknife standard error,
51    /// `sqrt((n - 1) / n * Σ (replicate_i - mean(replicates))²)`.
52    std_error: f64,
53    /// The statistic evaluated on each leave-one-out subsample, in index order.
54    replicates: Vec<f64>,
55}
56
57impl JackknifeEstimate {
58    /// Returns the statistic evaluated on the full sample.
59    ///
60    /// # Returns
61    ///
62    /// The full-sample estimate `stat(data)`.
63    ///
64    /// # Examples
65    ///
66    /// ```
67    /// use stats_claw::resampling::jackknife_statistic;
68    ///
69    /// let data = [1.0, 2.0, 3.0];
70    /// let est = jackknife_statistic(&data, |s| s.iter().copied().fold(f64::MIN, f64::max))?;
71    /// assert!((est.estimate() - 3.0).abs() < 1e-12, "estimate was {}", est.estimate());
72    /// # Ok::<(), stats_claw::error::Error>(())
73    /// ```
74    #[must_use]
75    pub const fn estimate(&self) -> f64 {
76        self.estimate
77    }
78
79    /// Returns the jackknife bias estimate.
80    ///
81    /// # Returns
82    ///
83    /// `(n - 1) * (mean(replicates) - estimate)`; subtract it from
84    /// [`estimate`](Self::estimate) for the bias-corrected value.
85    ///
86    /// # Examples
87    ///
88    /// ```
89    /// use stats_claw::resampling::jackknife_statistic;
90    ///
91    /// let data = [2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0];
92    /// let mean = |s: &[f64]| s.iter().sum::<f64>() / f64::from(u32::try_from(s.len()).unwrap_or(0));
93    /// let est = jackknife_statistic(&data, mean)?;
94    /// // The jackknife bias of the mean is identically zero.
95    /// assert!(est.bias().abs() < 1e-12, "bias was {}", est.bias());
96    /// # Ok::<(), stats_claw::error::Error>(())
97    /// ```
98    #[must_use]
99    pub const fn bias(&self) -> f64 {
100        self.bias
101    }
102
103    /// Returns the jackknife standard error.
104    ///
105    /// # Returns
106    ///
107    /// `sqrt((n - 1) / n * Σ (replicate_i - mean(replicates))²)`.
108    ///
109    /// # Examples
110    ///
111    /// ```
112    /// use stats_claw::resampling::jackknife_statistic;
113    ///
114    /// let data = [2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0];
115    /// let mean = |s: &[f64]| s.iter().sum::<f64>() / f64::from(u32::try_from(s.len()).unwrap_or(0));
116    /// let est = jackknife_statistic(&data, mean)?;
117    /// // Jackknife SE of the mean equals the classic sd(ddof=1)/sqrt(n).
118    /// assert!((est.std_error() - 0.755_928_946_018_454_4).abs() < 1e-12, "se was {}", est.std_error());
119    /// # Ok::<(), stats_claw::error::Error>(())
120    /// ```
121    #[must_use]
122    pub const fn std_error(&self) -> f64 {
123        self.std_error
124    }
125
126    /// Returns the per-subsample replicates.
127    ///
128    /// # Returns
129    ///
130    /// A borrowed slice of the statistic evaluated on each leave-one-out
131    /// subsample, in index order; its length equals the input sample size.
132    ///
133    /// # Examples
134    ///
135    /// ```
136    /// use stats_claw::resampling::jackknife_statistic;
137    ///
138    /// let data = [1.0, 2.0, 3.0, 4.0];
139    /// let est = jackknife_statistic(&data, |s| s.iter().sum::<f64>())?;
140    /// assert_eq!(est.replicates().len(), data.len(), "one replicate per observation");
141    /// // Leaving out the largest value yields the smallest sum.
142    /// assert!((est.replicates()[3] - 6.0).abs() < 1e-12, "was {}", est.replicates()[3]);
143    /// # Ok::<(), stats_claw::error::Error>(())
144    /// ```
145    #[must_use]
146    pub fn replicates(&self) -> &[f64] {
147        &self.replicates
148    }
149}
150
151/// Builds the leave-one-out index sets for a sample of size `n`.
152///
153/// Returns `n` index vectors; set `i` contains every index in `0..n` except `i`,
154/// in ascending order. These are the subsamples the jackknife evaluates a
155/// statistic on.
156///
157/// # Arguments
158///
159/// * `n` — the sample size; must be `>= 2` (a jackknife needs at least two
160///   observations to leave one out and still have data).
161///
162/// # Returns
163///
164/// A length-`n` vector whose `i`-th entry lists the `n - 1` indices `0..n` with
165/// `i` removed.
166///
167/// # Errors
168///
169/// Returns [`Error::InsufficientData`] when `n < 2`.
170///
171/// # Examples
172///
173/// ```
174/// use stats_claw::resampling::jackknife_indices;
175///
176/// assert_eq!(jackknife_indices(3)?, vec![vec![1, 2], vec![0, 2], vec![0, 1]]);
177/// # Ok::<(), stats_claw::error::Error>(())
178/// ```
179pub fn jackknife_indices(n: usize) -> Result<Vec<Vec<usize>>> {
180    if n < 2 {
181        return Err(Error::InsufficientData);
182    }
183    Ok((0..n)
184        .map(|i| (0..n).filter(|&j| j != i).collect())
185        .collect())
186}
187
188/// Computes the jackknife bias and standard error of `stat` over `data`.
189///
190/// Evaluates `stat` on the full sample and on each leave-one-out subsample, then
191/// combines the replicates into the classical jackknife summaries (Efron &
192/// Tibshirani, 1993, §10.2). Fully deterministic — no random numbers are drawn.
193///
194/// # Arguments
195///
196/// * `data` — the observed sample; must contain at least two values.
197/// * `stat` — the statistic to jackknife, mapping a sample view to a scalar.
198///
199/// # Returns
200///
201/// A [`JackknifeEstimate`] holding the full-sample estimate, the bias and
202/// standard-error estimates, and the per-subsample replicates.
203///
204/// # Errors
205///
206/// Returns [`Error::InsufficientData`] when `data` has fewer than two elements.
207///
208/// # Examples
209///
210/// ```
211/// use stats_claw::resampling::jackknife_statistic;
212///
213/// let data = [2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0];
214/// let mean = |s: &[f64]| s.iter().sum::<f64>() / f64::from(u32::try_from(s.len()).unwrap_or(0));
215/// let est = jackknife_statistic(&data, mean)?;
216/// assert!((est.estimate() - 5.0).abs() < 1e-12, "estimate was {}", est.estimate());
217/// // Jackknife SE of the mean equals the classic sd(ddof=1)/sqrt(n).
218/// assert!((est.std_error() - 0.755_928_946_018_454_4).abs() < 1e-12, "se was {}", est.std_error());
219/// assert_eq!(est.replicates().len(), data.len(), "one replicate per observation");
220/// # Ok::<(), stats_claw::error::Error>(())
221/// ```
222pub fn jackknife_statistic(
223    data: &[f64],
224    stat: impl Fn(&[f64]) -> f64,
225) -> Result<JackknifeEstimate> {
226    let n = data.len();
227    if n < 2 {
228        return Err(Error::InsufficientData);
229    }
230    let estimate = stat(data);
231    let mut replicates = Vec::with_capacity(n);
232    for i in 0..n {
233        let sample: Vec<f64> = data
234            .iter()
235            .enumerate()
236            .filter(|&(j, _)| j != i)
237            .map(|(_, &value)| value)
238            .collect();
239        replicates.push(stat(&sample));
240    }
241    let n_f = count_to_f64(n);
242    // `replicates.len() == n`, so this shared mean matches the old inline
243    // `sum / n_f` bit-for-bit; the jackknife bias/SE below keep their own
244    // (non-`sample_variance`) `(n-1)/n` scaling.
245    let mean_rep = crate::numeric::mean(&replicates);
246    let bias = (n_f - 1.0) * (mean_rep - estimate);
247    let sum_sq: f64 = replicates
248        .iter()
249        .map(|&r| {
250            let deviation = r - mean_rep;
251            deviation * deviation
252        })
253        .sum();
254    let std_error = ((n_f - 1.0) / n_f * sum_sq).sqrt();
255    Ok(JackknifeEstimate {
256        estimate,
257        bias,
258        std_error,
259        replicates,
260    })
261}
262
263impl JackknifeResampling {
264    /// Computes the jackknife bias and standard error of `stat` over `data`.
265    ///
266    /// Convenience wrapper that lets the [`JackknifeResampling`] scheme
267    /// drive the numerics directly; it delegates to [`jackknife_statistic`] and
268    /// shares its contract exactly.
269    ///
270    /// # Arguments
271    ///
272    /// * `data` — the observed sample; must contain at least two values.
273    /// * `stat` — the statistic to jackknife, mapping a sample view to a scalar.
274    ///
275    /// # Returns
276    ///
277    /// A [`JackknifeEstimate`] with the full-sample estimate, bias,
278    /// standard error, and per-subsample replicates.
279    ///
280    /// # Errors
281    ///
282    /// Returns [`Error::InsufficientData`] when `data` has fewer than two
283    /// elements.
284    ///
285    /// # Examples
286    ///
287    /// ```
288    /// use stats_claw::resampling::JackknifeResampling;
289    ///
290    /// let scheme = JackknifeResampling::default();
291    /// let data = [1.0, 2.0, 3.0, 4.0];
292    /// let est = scheme.estimate(&data, |s| s.iter().copied().fold(f64::MIN, f64::max))?;
293    /// assert!((est.estimate() - 4.0).abs() < 1e-12, "max was {}", est.estimate());
294    /// # Ok::<(), stats_claw::error::Error>(())
295    /// ```
296    pub fn estimate(
297        &self,
298        data: &[f64],
299        stat: impl Fn(&[f64]) -> f64,
300    ) -> Result<JackknifeEstimate> {
301        jackknife_statistic(data, stat)
302    }
303}
304
305/// Kani formal-verification harnesses for the deterministic jackknife index
306/// construction.
307///
308/// [`jackknife_indices`] draws no random numbers, so these prove its
309/// input-validation and index-safety properties over symbolic and small-fixed
310/// sizes rather than the sampled fixtures the `#[cfg(test)]` suite uses. Compiled
311/// only under `cargo kani` (behind `#[cfg(kani)]`); invisible to normal
312/// build/test/clippy. Run e.g. with
313/// `cargo kani -Z stubbing -p stats-claw --harness resampling_jackknife_rejects_small_n`.
314#[cfg(kani)]
315mod verification {
316    use super::{Error, jackknife_indices};
317
318    /// Proves the input-validation path: for *every* symbolic `n < 2`,
319    /// [`jackknife_indices`] returns [`Error::InsufficientData`] and never panics —
320    /// a single observation has no held-out complement.
321    ///
322    /// The `#[kani::unwind(2)]` bound caps the (unreachable-on-feasible-paths)
323    /// index-building loop: with `n < 2` the guard returns before it, and CBMC
324    /// discharges the over-unwinding of the infeasible `n >= 2` branch vacuously.
325    #[kani::proof]
326    #[kani::unwind(2)]
327    fn resampling_jackknife_rejects_small_n() {
328        let n: usize = kani::any();
329        kani::assume(n < 2);
330        let result = jackknife_indices(n);
331        assert!(
332            matches!(result, Err(Error::InsufficientData)),
333            "n < 2 must be rejected with InsufficientData"
334        );
335    }
336
337    /// Proves the leave-one-out index sets are total and in bounds for `n = 3`:
338    /// exactly `n` sets, each of size `n - 1`, every index in `0..n`, and set `i`
339    /// omitting its own index `i`. This is the index-safety guarantee the jackknife
340    /// statistic loop relies on when it views each subsample.
341    #[kani::proof]
342    #[kani::unwind(5)]
343    fn resampling_jackknife_indices_in_bounds() {
344        const N: usize = 3;
345        let result = jackknife_indices(N);
346        assert!(result.is_ok(), "n >= 2 must produce leave-one-out sets");
347        if let Ok(sets) = result {
348            assert!(sets.len() == N, "expected one set per observation");
349            for (i, set) in sets.iter().enumerate() {
350                assert!(set.len() == N - 1, "each set must omit exactly one index");
351                for &j in set {
352                    assert!(j < N, "jackknife index {j} escaped 0..N");
353                    assert!(j != i, "set {i} must not contain its own index");
354                }
355            }
356        }
357    }
358}
359
360#[cfg(test)]
361mod tests;
stats_claw/resampling/jackknife.rs

stats_claw/resampling/
jackknife.rs
