qs_backtest/evaluation/
stats.rs

1use super::{BootstrapConfig, ConfidenceInterval, MetricValue};
2
3pub(crate) fn mean(values: &[f64]) -> Option<f64> {
4    (!values.is_empty()).then(|| values.iter().sum::<f64>() / values.len() as f64)
5}
6
7pub(crate) fn median(values: &[f64]) -> Option<f64> {
8    if values.is_empty() {
9        return None;
10    }
11
12    let mut sorted = values.to_vec();
13    sorted.sort_by(f64::total_cmp);
14    let midpoint = sorted.len() / 2;
15    Some(if sorted.len().is_multiple_of(2) {
16        (sorted[midpoint - 1] + sorted[midpoint]) / 2.0
17    } else {
18        sorted[midpoint]
19    })
20}
21
22pub(crate) fn sample_standard_deviation(values: &[f64]) -> Option<f64> {
23    if values.len() < 2 {
24        return None;
25    }
26
27    let average = mean(values)?;
28    let variance = values
29        .iter()
30        .map(|value| (value - average).powi(2))
31        .sum::<f64>()
32        / (values.len() - 1) as f64;
33    Some(variance.sqrt())
34}
35
36/// Computes a two-sided Wilson score interval for a binomial proportion.
37///
38/// `confidence_level` must be strictly between zero and one. A zero-trial input
39/// is reported as insufficient data rather than represented as `0 / 0`.
40pub fn wilson_interval(
41    successes: usize,
42    trials: usize,
43    confidence_level: f64,
44) -> MetricValue<ConfidenceInterval> {
45    if successes > trials {
46        return MetricValue::invalid_input("successes cannot exceed trials");
47    }
48    if !confidence_level.is_finite() || !(0.0..1.0).contains(&confidence_level) {
49        return MetricValue::invalid_input("confidence_level must be finite and between 0 and 1");
50    }
51    if trials == 0 {
52        return MetricValue::insufficient_data("at least one trial is required");
53    }
54
55    let n = trials as f64;
56    let estimate = successes as f64 / n;
57    let z = inverse_standard_normal(0.5 + confidence_level / 2.0);
58    let z_squared = z * z;
59    let denominator = 1.0 + z_squared / n;
60    let center = (estimate + z_squared / (2.0 * n)) / denominator;
61    let margin =
62        z * ((estimate * (1.0 - estimate) / n) + z_squared / (4.0 * n * n)).sqrt() / denominator;
63
64    MetricValue::available(ConfidenceInterval {
65        estimate,
66        lower: (center - margin).max(0.0),
67        upper: (center + margin).min(1.0),
68        confidence_level,
69    })
70}
71
72/// Computes a deterministic percentile-bootstrap confidence interval for a mean.
73///
74/// Resampling uses an internal SplitMix64 generator, making the output stable for
75/// the same values, ordering, configuration, and crate version without adding a
76/// random-number dependency.
77pub fn bootstrap_mean_confidence(
78    values: &[f64],
79    config: BootstrapConfig,
80) -> MetricValue<ConfidenceInterval> {
81    if values.iter().any(|value| !value.is_finite()) {
82        return MetricValue::invalid_input("bootstrap values must all be finite");
83    }
84    if !config.confidence_level.is_finite() || !(0.0..1.0).contains(&config.confidence_level) {
85        return MetricValue::invalid_input("confidence_level must be finite and between 0 and 1");
86    }
87    if config.samples == 0 {
88        return MetricValue::invalid_input("bootstrap samples must be greater than zero");
89    }
90    if values.len() < config.minimum_sample_size {
91        return MetricValue::insufficient_data(format!(
92            "at least {} observations are required for bootstrap confidence",
93            config.minimum_sample_size
94        ));
95    }
96    if values.is_empty() {
97        return MetricValue::insufficient_data("at least one observation is required");
98    }
99
100    let estimate = mean(values).expect("non-empty values checked above");
101    if !estimate.is_finite() {
102        return MetricValue::invalid_input("bootstrap mean exceeds the finite f64 range");
103    }
104
105    let mut rng = SplitMix64::new(config.seed);
106    let mut bootstrap_means = Vec::with_capacity(config.samples);
107    for _ in 0..config.samples {
108        let mut total = 0.0;
109        for _ in values {
110            total += values[rng.index(values.len())];
111        }
112        let bootstrap_mean = total / values.len() as f64;
113        if !bootstrap_mean.is_finite() {
114            return MetricValue::invalid_input(
115                "a bootstrap sample mean exceeds the finite f64 range",
116            );
117        }
118        bootstrap_means.push(bootstrap_mean);
119    }
120    bootstrap_means.sort_by(f64::total_cmp);
121
122    let alpha = 1.0 - config.confidence_level;
123    let lower = quantile_sorted(&bootstrap_means, alpha / 2.0);
124    let upper = quantile_sorted(&bootstrap_means, 1.0 - alpha / 2.0);
125
126    MetricValue::available(ConfidenceInterval {
127        estimate,
128        lower,
129        upper,
130        confidence_level: config.confidence_level,
131    })
132}
133
134pub(crate) fn quantile_sorted(sorted: &[f64], probability: f64) -> f64 {
135    debug_assert!(!sorted.is_empty());
136    let index = probability * (sorted.len() - 1) as f64;
137    let lower = index.floor() as usize;
138    let upper = index.ceil() as usize;
139    if lower == upper {
140        sorted[lower]
141    } else {
142        let weight = index - lower as f64;
143        sorted[lower] * (1.0 - weight) + sorted[upper] * weight
144    }
145}
146
147// Peter J. Acklam's inverse-normal approximation. Accuracy is more than
148// sufficient for confidence bounds and avoids a statistics dependency.
149fn inverse_standard_normal(probability: f64) -> f64 {
150    const A: [f64; 6] = [
151        -3.969_683_028_665_376e1,
152        2.209_460_984_245_205e2,
153        -2.759_285_104_469_687e2,
154        1.383_577_518_672_69e2,
155        -3.066_479_806_614_716e1,
156        2.506_628_277_459_239,
157    ];
158    const B: [f64; 5] = [
159        -5.447_609_879_822_406e1,
160        1.615_858_368_580_409e2,
161        -1.556_989_798_598_866e2,
162        6.680_131_188_771_972e1,
163        -1.328_068_155_288_572e1,
164    ];
165    const C: [f64; 6] = [
166        -7.784_894_002_430_293e-3,
167        -3.223_964_580_411_365e-1,
168        -2.400_758_277_161_838,
169        -2.549_732_539_343_734,
170        4.374_664_141_464_968,
171        2.938_163_982_698_783,
172    ];
173    const D: [f64; 4] = [
174        7.784_695_709_041_462e-3,
175        3.224_671_290_700_398e-1,
176        2.445_134_137_142_996,
177        3.754_408_661_907_416,
178    ];
179    const LOW: f64 = 0.024_25;
180    const HIGH: f64 = 1.0 - LOW;
181
182    if probability < LOW {
183        let q = (-2.0 * probability.ln()).sqrt();
184        (((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5])
185            / ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0)
186    } else if probability <= HIGH {
187        let q = probability - 0.5;
188        let r = q * q;
189        (((((A[0] * r + A[1]) * r + A[2]) * r + A[3]) * r + A[4]) * r + A[5]) * q
190            / (((((B[0] * r + B[1]) * r + B[2]) * r + B[3]) * r + B[4]) * r + 1.0)
191    } else {
192        let q = (-2.0 * (1.0 - probability).ln()).sqrt();
193        -(((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5])
194            / ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0)
195    }
196}
197
198struct SplitMix64 {
199    state: u64,
200}
201
202impl SplitMix64 {
203    fn new(seed: u64) -> Self {
204        Self { state: seed }
205    }
206
207    fn next_u64(&mut self) -> u64 {
208        self.state = self.state.wrapping_add(0x9E37_79B9_7F4A_7C15);
209        let mut value = self.state;
210        value = (value ^ (value >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
211        value = (value ^ (value >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
212        value ^ (value >> 31)
213    }
214
215    fn index(&mut self, upper_bound: usize) -> usize {
216        ((self.next_u64() as u128 * upper_bound as u128) >> 64) as usize
217    }
218}
qs_backtest/evaluation/stats.rs

qs_backtest/evaluation/
stats.rs
