use rand::Rng;
use rand::SeedableRng;
use rand::rngs::StdRng;
pub fn paired_bootstrap_ci(deltas: &[f32], n_resamples: usize, ci: f64, seed: u64) -> (f32, f32) {
if deltas.is_empty() || n_resamples == 0 {
return (0.0, 0.0);
}
let mut rng = StdRng::seed_from_u64(seed);
let n = deltas.len();
let mut means: Vec<f32> = Vec::with_capacity(n_resamples);
for _ in 0..n_resamples {
let mut sum = 0.0_f64;
for _ in 0..n {
let idx = rng.random_range(0..n);
sum += deltas[idx] as f64;
}
means.push((sum / n as f64) as f32);
}
means.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
let tail = (1.0 - ci) / 2.0;
let lo = percentile(&means, tail);
let hi = percentile(&means, 1.0 - tail);
(lo, hi)
}
fn percentile(sorted: &[f32], p: f64) -> f32 {
let last = sorted.len() - 1;
let rank = (p * last as f64).round() as usize;
sorted[rank.min(last)]
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn all_positive_deltas_collapse_to_one() {
let deltas = vec![1.0_f32; 50];
let (lo, hi) = paired_bootstrap_ci(&deltas, 2000, 0.95, 7);
assert!((lo - 1.0).abs() < 1e-6, "lo was {lo}");
assert!((hi - 1.0).abs() < 1e-6, "hi was {hi}");
}
#[test]
fn all_zero_deltas_give_zero_interval() {
let deltas = vec![0.0_f32; 50];
let (lo, hi) = paired_bootstrap_ci(&deltas, 2000, 0.95, 7);
assert_eq!(lo, 0.0);
assert_eq!(hi, 0.0);
}
#[test]
fn mixed_deltas_bracket_the_true_mean() {
let deltas: Vec<f32> = (0..40).map(|i| if i % 2 == 0 { 1.0 } else { -0.5 }).collect();
let mean = deltas.iter().sum::<f32>() / deltas.len() as f32;
let (lo, hi) = paired_bootstrap_ci(&deltas, 2000, 0.95, 7);
assert!(lo < hi, "expected lo < hi, got ({lo}, {hi})");
assert!(lo <= mean && mean <= hi, "mean {mean} not in [{lo}, {hi}]");
}
#[test]
fn is_deterministic_for_a_fixed_seed() {
let deltas: Vec<f32> = (0..30).map(|i| (i as f32) * 0.01 - 0.15).collect();
let a = paired_bootstrap_ci(&deltas, 2000, 0.95, 42);
let b = paired_bootstrap_ci(&deltas, 2000, 0.95, 42);
assert_eq!(a, b);
}
}