pub(crate) fn mean(samples: &[f64]) -> f64 {
if samples.is_empty() {
return 0.0;
}
samples.iter().sum::<f64>() / samples.len() as f64
}
pub(crate) fn std_dev(samples: &[f64]) -> f64 {
batuta_common::math::std_dev(samples)
}
pub(crate) fn cv(samples: &[f64]) -> f64 {
let m = mean(samples);
if m <= 0.0 || samples.len() < 2 {
return 0.0;
}
(std_dev(samples) / m) * 100.0
}
pub(crate) fn t_test(a: &[f64], b: &[f64]) -> f64 {
if a.len() < 2 || b.len() < 2 {
return 1.0; }
let mean_a = mean(a);
let mean_b = mean(b);
let var_a = std_dev(a).powi(2);
let var_b = std_dev(b).powi(2);
let n_a = a.len() as f64;
let n_b = b.len() as f64;
let se = ((var_a / n_a) + (var_b / n_b)).sqrt();
if se == 0.0 {
return 1.0;
}
let t = (mean_a - mean_b).abs() / se;
let num = ((var_a / n_a) + (var_b / n_b)).powi(2);
let denom = (var_a / n_a).powi(2) / (n_a - 1.0) + (var_b / n_b).powi(2) / (n_b - 1.0);
let df = if denom > 0.0 { num / denom } else { 1.0 };
if df > 30.0 {
2.0 * (1.0 - normal_cdf(t))
} else {
let adjusted_t = t * (1.0 + 0.5 / df).sqrt();
2.0 * (1.0 - normal_cdf(adjusted_t))
}
}
fn normal_cdf(x: f64) -> f64 {
let a1 = 0.254829592;
let a2 = -0.284496736;
let a3 = 1.421413741;
let a4 = -1.453152027;
let a5 = 1.061405429;
let p = 0.3275911;
let sign = if x < 0.0 { -1.0 } else { 1.0 };
let x = x.abs() / std::f64::consts::SQRT_2;
let t = 1.0 / (1.0 + p * x);
let y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * (-x * x).exp();
0.5 * (1.0 + sign * y)
}