use std::f64::consts::{LN_2, PI, SQRT_2};
use num_traits::{One, Signed, Zero};
use super::common::{
check_confidence, check_unit_open, chi_squared_sf, ex, f_sf, f_sf_rational, invalid, qi, qu,
};
use super::data::{self, Q};
use super::hypothesis::{self, Alternative, PValue, TestResult, p_value_accessors};
use super::numdist::norm::{pdf as norm_pdf, sf as norm_sf};
use super::regression::ols;
use crate::api::context::Context;
use crate::api::expr::Ex;
use crate::base::errors::SymplexError;
use crate::base::interval::{Bounds, Interval};
use crate::base::numeric::ratio_to_f64;
use crate::domains::exact_matrix::QMatrix;
use crate::domains::optimize::{RootOpts, brent_root, grow_bracket};
use crate::output::codegen::numeric_rt::{erf, erfc};
fn failed(op: &'static str, reason: impl Into<String>) -> SymplexError {
SymplexError::computation_failed(op, reason)
}
fn to_f64(op: &'static str, q: &Q) -> Result<f64, SymplexError> {
ratio_to_f64(q).ok_or_else(|| failed(op, format!("{q} does not fit in an f64")))
}
fn centred_ss(x: &[Q]) -> Q {
data::sum_of_squares(x).unwrap_or_else(|_| Q::zero())
}
pub(crate) struct GroupSums {
pub(crate) ss_between: Q,
pub(crate) ss_within: Q,
pub(crate) n: usize,
pub(crate) k: usize,
}
pub(crate) fn sums_of_squares(
op: &'static str,
groups: &[Vec<Q>],
) -> Result<GroupSums, SymplexError> {
let k = groups.len();
if k < 2 {
return Err(invalid(op, "at least two groups are needed"));
}
if groups.iter().any(Vec::is_empty) {
return Err(invalid(op, "every group must be non-empty"));
}
let all: Vec<Q> = groups.iter().flatten().cloned().collect();
let grand = data::mean(&all)?;
let mut ss_between = Q::zero();
let mut ss_within = Q::zero();
for g in groups {
let m = data::mean(g)?;
let d = &m - &grand;
ss_between += qu(g.len()) * &d * &d;
ss_within += data::sum_of_squares(g)?;
}
Ok(GroupSums {
ss_between,
ss_within,
n: all.len(),
k,
})
}
#[derive(Clone, Debug, PartialEq)]
pub struct AnovaResult {
pub f: Q,
pub df_between: usize,
pub df_within: usize,
pub p_value: Ex,
pub ss_between: Q,
pub ss_within: Q,
pub eta_squared: Q,
}
impl AnovaResult {
pub fn p_value_f64(&self) -> Result<f64, SymplexError> {
self.p_value.eval_f64()
}
}
impl PValue for AnovaResult {
fn p_value_ex(&self) -> &Ex {
&self.p_value
}
}
p_value_accessors!(AnovaResult);
pub fn anova_one_way(ctx: &Context, groups: &[Vec<Q>]) -> Result<AnovaResult, SymplexError> {
const OP: &str = "anova_one_way";
let GroupSums {
ss_between,
ss_within,
n,
k,
} = sums_of_squares(OP, groups)?;
if n <= k {
return Err(invalid(
OP,
"at least one group needs more than one observation",
));
}
if ss_within.is_zero() {
return Err(invalid(OP, "the within-group variance is zero"));
}
let (df_between, df_within) = (k - 1, n - k);
let f = (&ss_between / qu(df_between)) / (&ss_within / qu(df_within));
let total = &ss_between + &ss_within;
let eta_squared = &ss_between / &total;
Ok(AnovaResult {
p_value: f_sf(ctx, df_between, df_within, &f),
f,
df_between,
df_within,
ss_between,
ss_within,
eta_squared,
})
}
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum Source {
FactorA,
FactorB,
Interaction,
Conditions,
Subjects,
Residual,
Total,
}
#[derive(Clone, Debug, PartialEq)]
pub struct AnovaRow {
pub source: Source,
pub ss: Q,
pub df: usize,
pub ms: Option<Q>,
pub f: Option<Q>,
pub p_value: Option<Ex>,
pub eta_squared: Option<Q>,
pub partial_eta_squared: Option<Q>,
}
impl AnovaRow {
fn effect(
ctx: &Context,
source: Source,
ss: Q,
df: usize,
ss_resid: &Q,
df_resid: usize,
ss_total: &Q,
) -> Self {
let test = FTest::of(ctx, &ss, df, ss_resid, df_resid);
Self::with_test(source, ss, df, test, ss_resid, ss_total)
}
fn with_test(
source: Source,
ss: Q,
df: usize,
test: FTest,
ss_resid: &Q,
ss_total: &Q,
) -> Self {
Self {
source,
eta_squared: Some(&ss / ss_total),
partial_eta_squared: Some(&ss / (&ss + ss_resid)),
ms: Some(&ss / qu(df)),
ss,
df,
f: Some(test.f),
p_value: Some(test.p_value),
}
}
fn untested(source: Source, ss: Q, df: usize) -> Self {
Self {
source,
ms: Some(&ss / qu(df)),
ss,
df,
f: None,
p_value: None,
eta_squared: None,
partial_eta_squared: None,
}
}
fn total(ss: Q, df: usize) -> Self {
Self {
source: Source::Total,
ss,
df,
ms: None,
f: None,
p_value: None,
eta_squared: None,
partial_eta_squared: None,
}
}
pub fn p_value_f64(&self) -> Result<f64, SymplexError> {
self.tested_p_value("AnovaRow::p_value_f64")?.eval_f64()
}
pub fn p_value_log10(&self) -> Result<f64, SymplexError> {
hypothesis::p_value_log10_of(self.tested_p_value("AnovaRow::p_value_log10")?)
}
pub fn p_value_ln(&self) -> Result<f64, SymplexError> {
hypothesis::p_value_ln_of(self.tested_p_value("AnovaRow::p_value_ln")?)
}
pub fn p_value_decimal(&self, digits: u32) -> Result<String, SymplexError> {
self.tested_p_value("AnovaRow::p_value_decimal")?
.eval_decimal(digits)
}
fn tested_p_value(&self, op: &'static str) -> Result<&Ex, SymplexError> {
self.p_value
.as_ref()
.ok_or_else(|| invalid(op, format!("the {:?} row has no F test", self.source)))
}
}
struct FTest {
f: Q,
p_value: Ex,
}
impl FTest {
fn of(ctx: &Context, ss: &Q, df: usize, ss_resid: &Q, df_resid: usize) -> Self {
let f = (ss / qu(df)) / (ss_resid / qu(df_resid));
let p_value = f_sf(ctx, df, df_resid, &f);
Self { f, p_value }
}
}
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct Observation {
pub a: usize,
pub b: usize,
pub y: Q,
}
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct TwoWayData {
a_levels: usize,
b_levels: usize,
cells: Vec<Vec<Vec<Q>>>,
}
impl TwoWayData {
pub fn from_cells(cells: Vec<Vec<Vec<Q>>>) -> Result<Self, SymplexError> {
const OP: &str = "TwoWayData::from_cells";
let a_levels = cells.len();
if a_levels < 2 {
return Err(invalid(OP, "factor A needs at least two levels"));
}
let b_levels = cells.first().map_or(0, Vec::len);
if b_levels < 2 {
return Err(invalid(OP, "factor B needs at least two levels"));
}
for (a, row) in cells.iter().enumerate() {
if row.len() != b_levels {
return Err(invalid(
OP,
format!(
"level {a} of A has {} cells, expected {b_levels}",
row.len()
),
));
}
if let Some((b, _)) = row.iter().enumerate().find(|(_, c)| c.is_empty()) {
return Err(invalid(
OP,
format!("cell (A = {a}, B = {b}) has no observations"),
));
}
}
Ok(Self {
a_levels,
b_levels,
cells,
})
}
pub fn from_i64(cells: &[&[&[i64]]]) -> Result<Self, SymplexError> {
Self::from_cells(
cells
.iter()
.map(|row| row.iter().map(|c| data::from_i64(c)).collect())
.collect(),
)
}
pub fn from_long(rows: &[Observation]) -> Result<Self, SymplexError> {
const OP: &str = "TwoWayData::from_long";
if rows.is_empty() {
return Err(invalid(OP, "no observations"));
}
let a_levels = rows.iter().map(|o| o.a).max().unwrap_or(0).checked_add(1);
let b_levels = rows.iter().map(|o| o.b).max().unwrap_or(0).checked_add(1);
let (a_levels, b_levels) = match (a_levels, b_levels) {
(Some(a), Some(b)) if a.checked_mul(b).is_some_and(|ab| ab <= rows.len()) => (a, b),
_ => {
return Err(invalid(
OP,
"more combinations of levels than observations: some combination has none",
));
}
};
let mut cells = vec![vec![Vec::new(); b_levels]; a_levels];
for o in rows {
if let Some(cell) = cells.get_mut(o.a).and_then(|row| row.get_mut(o.b)) {
cell.push(o.y.clone());
}
}
Self::from_cells(cells).map_err(|e| match e {
SymplexError::InvalidArgument { reason, .. } => invalid(OP, reason),
other => other,
})
}
#[must_use]
pub fn a_levels(&self) -> usize {
self.a_levels
}
#[must_use]
pub fn b_levels(&self) -> usize {
self.b_levels
}
#[must_use]
pub fn n_obs(&self) -> usize {
self.cells.iter().flatten().map(Vec::len).sum()
}
#[must_use]
pub fn cell(&self, a: usize, b: usize) -> Option<&[Q]> {
self.cells.get(a)?.get(b).map(Vec::as_slice)
}
#[must_use]
pub fn cells(&self) -> &[Vec<Vec<Q>>] {
&self.cells
}
#[must_use]
pub fn is_balanced(&self) -> bool {
let mut sizes = self.cells.iter().flatten().map(Vec::len);
match sizes.next() {
Some(first) => sizes.all(|s| s == first),
None => true,
}
}
fn flatten(&self) -> Vec<Q> {
self.cells.iter().flatten().flatten().cloned().collect()
}
fn a_slice(&self, a: usize) -> Vec<Q> {
self.cells
.get(a)
.map(|row| row.iter().flatten().cloned().collect())
.unwrap_or_default()
}
fn b_slice(&self, b: usize) -> Vec<Q> {
self.cells
.iter()
.filter_map(|row| row.get(b))
.flatten()
.cloned()
.collect()
}
}
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum SsType {
TypeI,
TypeII,
TypeIII,
}
#[derive(Clone, Debug, PartialEq)]
pub struct TwoWayAnova {
pub ss_type: SsType,
pub factor_a: AnovaRow,
pub factor_b: AnovaRow,
pub interaction: AnovaRow,
pub residual: AnovaRow,
pub total: AnovaRow,
pub grand_mean: Q,
pub cell_means: Vec<Vec<Q>>,
}
impl TwoWayAnova {
#[must_use]
pub fn rows(&self) -> Vec<&AnovaRow> {
vec![
&self.factor_a,
&self.factor_b,
&self.interaction,
&self.residual,
&self.total,
]
}
}
#[derive(Clone, Copy)]
struct Terms {
a: bool,
b: bool,
ab: bool,
}
#[derive(Clone, Copy)]
enum Coding {
Treatment,
Sum,
}
fn codes(level: usize, n_levels: usize, coding: Coding) -> Vec<Q> {
(1..n_levels)
.map(|i| match coding {
Coding::Treatment => {
if level == i {
Q::one()
} else {
Q::zero()
}
}
Coding::Sum => {
if level == i {
Q::one()
} else if level == 0 {
-Q::one()
} else {
Q::zero()
}
}
})
.collect()
}
fn reduced_rss(
op: &'static str,
data: &TwoWayData,
y: &[Q],
terms: Terms,
coding: Coding,
) -> Result<Q, SymplexError> {
let mut rows: Vec<Vec<Q>> = Vec::with_capacity(y.len());
for (a, row) in data.cells.iter().enumerate() {
let ca = codes(a, data.a_levels, coding);
for (b, cell) in row.iter().enumerate() {
let cb = codes(b, data.b_levels, coding);
let mut r = Vec::new();
if terms.a {
r.extend(ca.iter().cloned());
}
if terms.b {
r.extend(cb.iter().cloned());
}
if terms.ab {
for x in &ca {
for z in &cb {
r.push(x * z);
}
}
}
rows.extend(std::iter::repeat_n(r, cell.len()));
}
}
ols(y, &rows, true).map(|fit| fit.ssr).map_err(|e| {
failed(
op,
format!("least-squares fit of a reduced model failed: {e}"),
)
})
}
pub fn anova_two_way(ctx: &Context, data: &TwoWayData) -> Result<TwoWayAnova, SymplexError> {
anova_two_way_with(ctx, data, SsType::TypeII)
}
pub fn anova_two_way_with(
ctx: &Context,
data: &TwoWayData,
ss_type: SsType,
) -> Result<TwoWayAnova, SymplexError> {
const OP: &str = "anova_two_way";
let (a, b) = (data.a_levels, data.b_levels);
let n = data.n_obs();
if n <= a * b {
return Err(invalid(
OP,
"every cell has a single observation: no residual degrees of freedom",
));
}
let df_resid = n - a * b;
let y = data.flatten();
let grand_mean = data::mean(&y)?;
let ss_total = centred_ss(&y);
if ss_total.is_zero() {
return Err(invalid(OP, "the response is constant"));
}
let cell_means = data
.cells
.iter()
.map(|row| {
row.iter()
.map(|c| data::mean(c))
.collect::<Result<Vec<_>, _>>()
})
.collect::<Result<Vec<_>, _>>()?;
let rss_full = data
.cells
.iter()
.flatten()
.fold(Q::zero(), |acc, c| acc + centred_ss(c));
if rss_full.is_zero() {
return Err(invalid(
OP,
"the within-cell variance is zero: F is undefined",
));
}
let rss_a = (0..a).fold(Q::zero(), |acc, i| acc + centred_ss(&data.a_slice(i)));
let rss_b = (0..b).fold(Q::zero(), |acc, j| acc + centred_ss(&data.b_slice(j)));
let additive = Terms {
a: true,
b: true,
ab: false,
};
let rss_ab = reduced_rss(OP, data, &y, additive, Coding::Treatment)?;
let (ss_a, ss_b) = match ss_type {
SsType::TypeI => (&ss_total - &rss_a, &rss_a - &rss_ab),
SsType::TypeII => (&rss_b - &rss_ab, &rss_a - &rss_ab),
SsType::TypeIII => {
let without_a = Terms {
a: false,
b: true,
ab: true,
};
let without_b = Terms {
a: true,
b: false,
ab: true,
};
let rss_no_a = reduced_rss(OP, data, &y, without_a, Coding::Sum)?;
let rss_no_b = reduced_rss(OP, data, &y, without_b, Coding::Sum)?;
(rss_no_a - &rss_full, rss_no_b - &rss_full)
}
};
let ss_ab = &rss_ab - &rss_full;
Ok(TwoWayAnova {
ss_type,
factor_a: AnovaRow::effect(
ctx,
Source::FactorA,
ss_a,
a - 1,
&rss_full,
df_resid,
&ss_total,
),
factor_b: AnovaRow::effect(
ctx,
Source::FactorB,
ss_b,
b - 1,
&rss_full,
df_resid,
&ss_total,
),
interaction: AnovaRow::effect(
ctx,
Source::Interaction,
ss_ab,
(a - 1) * (b - 1),
&rss_full,
df_resid,
&ss_total,
),
residual: AnovaRow::untested(Source::Residual, rss_full, df_resid),
total: AnovaRow::total(ss_total, n - 1),
grand_mean,
cell_means,
})
}
#[derive(Clone, Debug, PartialEq)]
pub struct Mauchly {
pub w: Q,
pub chi_squared: Ex,
pub df: usize,
pub p_value: Ex,
}
impl Mauchly {
pub fn chi_squared_f64(&self) -> Result<f64, SymplexError> {
self.chi_squared.eval_f64()
}
pub fn p_value_f64(&self) -> Result<f64, SymplexError> {
self.p_value.eval_f64()
}
}
impl PValue for Mauchly {
fn p_value_ex(&self) -> &Ex {
&self.p_value
}
}
p_value_accessors!(Mauchly);
#[derive(Clone, Debug, PartialEq)]
pub struct RepeatedMeasuresAnova {
pub n_subjects: usize,
pub n_conditions: usize,
pub conditions: AnovaRow,
pub subjects: AnovaRow,
pub error: AnovaRow,
pub total: AnovaRow,
pub f: Q,
pub p_value: Ex,
pub epsilon_gg: Q,
pub epsilon_hf: Option<Q>,
pub p_value_gg: Ex,
pub p_value_hf: Option<Ex>,
pub mauchly: Option<Mauchly>,
pub grand_mean: Q,
pub condition_means: Vec<Q>,
pub subject_means: Vec<Q>,
}
impl RepeatedMeasuresAnova {
pub fn p_value_f64(&self) -> Result<f64, SymplexError> {
self.p_value.eval_f64()
}
pub fn p_value_gg_f64(&self) -> Result<f64, SymplexError> {
self.p_value_gg.eval_f64()
}
pub fn p_value_hf_f64(&self) -> Result<Option<f64>, SymplexError> {
self.p_value_hf.as_ref().map(Ex::eval_f64).transpose()
}
pub fn p_value_gg_log10(&self) -> Result<f64, SymplexError> {
hypothesis::p_value_log10_of(&self.p_value_gg)
}
pub fn p_value_hf_log10(&self) -> Result<Option<f64>, SymplexError> {
self.p_value_hf
.as_ref()
.map(hypothesis::p_value_log10_of)
.transpose()
}
#[must_use]
pub fn rows(&self) -> Vec<&AnovaRow> {
vec![&self.conditions, &self.subjects, &self.error, &self.total]
}
}
impl PValue for RepeatedMeasuresAnova {
fn p_value_ex(&self) -> &Ex {
&self.p_value
}
}
p_value_accessors!(RepeatedMeasuresAnova);
fn mauchly(
ctx: &Context,
s: &[Vec<Q>],
trace: &Q,
n: usize,
k: usize,
) -> Result<Option<Mauchly>, SymplexError> {
const OP: &str = "anova_repeated_measures";
let d = k - 1;
let m = QMatrix::new(
(0..k)
.map(|row| {
(0..d)
.map(|j| {
if row == j {
Q::one()
} else if row == j + 1 {
-Q::one()
} else {
Q::zero()
}
})
.collect()
})
.collect(),
)
.map_err(|e| failed(OP, e.to_string()))?;
let s_mat = QMatrix::new(s.to_vec()).map_err(|e| failed(OP, e.to_string()))?;
let mt = m.transpose();
let mtsm = mt.matmul(&s_mat)?.matmul(&m)?;
let mtm = mt.matmul(&m)?;
let product = mtsm.det()? / mtm.det()?;
if !product.is_positive() {
return Ok(None);
}
let mean_eigenvalue = trace / qu(d);
let denominator = (0..d).fold(Q::one(), |acc, _| acc * &mean_eigenvalue);
let w = product / denominator;
let (dq, n1) = (qu(d), qu(n - 1));
let rho = Q::one() - (qi(2) * &dq * &dq + &dq + qi(2)) / (qi(6) * &dq * &n1);
if !rho.is_positive() {
return Ok(None);
}
let scale = &n1 * &dq * ρ
let omega2 = (&dq + qi(2))
* (&dq - qi(1))
* (&dq - qi(2))
* (qi(2) * &dq * &dq * &dq + qi(6) * &dq * &dq + qi(3) * &dq + qi(2))
/ (qi(288) * &scale * &scale);
let df = d * (d + 1) / 2 - 1;
let chi_squared = ex(ctx, &(-(&n1 * &rho))) * ex(ctx, &w).ln();
let p1 = chi_squared_sf(ctx, df, &chi_squared);
let p2 = chi_squared_sf(ctx, df + 4, &chi_squared);
let p_value = (&p1 + ex(ctx, &omega2) * (p2 - &p1)).min_with(&ctx.one());
Ok(Some(Mauchly {
w,
chi_squared,
df,
p_value,
}))
}
pub fn anova_repeated_measures(
ctx: &Context,
subjects_by_condition: &[Vec<Q>],
) -> Result<RepeatedMeasuresAnova, SymplexError> {
const OP: &str = "anova_repeated_measures";
let n = subjects_by_condition.len();
if n < 2 {
return Err(invalid(OP, "at least two subjects are needed"));
}
let k = subjects_by_condition.first().map_or(0, Vec::len);
if k < 2 {
return Err(invalid(OP, "at least two conditions are needed"));
}
if let Some((i, row)) = subjects_by_condition
.iter()
.enumerate()
.find(|(_, r)| r.len() != k)
{
return Err(invalid(
OP,
format!(
"subject {i} has {} observations, expected {k} (one per condition)",
row.len()
),
));
}
let all: Vec<Q> = subjects_by_condition.iter().flatten().cloned().collect();
let grand_mean = data::mean(&all)?;
let ss_total = centred_ss(&all);
let condition_means: Vec<Q> = (0..k)
.map(|j| {
subjects_by_condition
.iter()
.fold(Q::zero(), |acc, r| acc + &r[j])
/ qu(n)
})
.collect();
let subject_means: Vec<Q> = subjects_by_condition
.iter()
.map(|r| data::sum(r) / qu(k))
.collect();
let sq = |x: &Q| x * x;
let ss_conditions = qu(n)
* condition_means
.iter()
.fold(Q::zero(), |acc, m| acc + sq(&(m - &grand_mean)));
let ss_subjects = qu(k)
* subject_means
.iter()
.fold(Q::zero(), |acc, m| acc + sq(&(m - &grand_mean)));
let ss_error = &ss_total - &ss_conditions - &ss_subjects;
if !ss_error.is_positive() {
return Err(invalid(
OP,
"the error variance is zero (every subject's profile is a shift of the condition means): F is undefined",
));
}
let (df_c, df_s, df_e) = (k - 1, n - 1, (k - 1) * (n - 1));
let test = FTest::of(ctx, &ss_conditions, df_c, &ss_error, df_e);
let (f, p_value) = (test.f.clone(), test.p_value.clone());
let conditions = AnovaRow::with_test(
Source::Conditions,
ss_conditions,
df_c,
test,
&ss_error,
&ss_total,
);
let s: Vec<Vec<Q>> = (0..k)
.map(|a| {
(0..k)
.map(|b| {
subjects_by_condition.iter().fold(Q::zero(), |acc, r| {
acc + (&r[a] - &condition_means[a]) * (&r[b] - &condition_means[b])
}) / qu(n - 1)
})
.collect()
})
.collect();
let row_means: Vec<Q> = s.iter().map(|r| data::sum(r) / qu(k)).collect();
let all_mean = data::sum(&row_means) / qu(k);
let s_tilde: Vec<Vec<Q>> = (0..k)
.map(|a| {
(0..k)
.map(|b| &s[a][b] - &row_means[a] - &row_means[b] + &all_mean)
.collect()
})
.collect();
let trace = (0..k).fold(Q::zero(), |acc, a| acc + &s_tilde[a][a]);
let trace_sq = s_tilde
.iter()
.flatten()
.fold(Q::zero(), |acc, v| acc + sq(v));
if !trace_sq.is_positive() {
return Err(failed(
OP,
"the double-centred covariance vanished although SS_error > 0",
));
}
let epsilon_gg = sq(&trace) / (qu(df_c) * &trace_sq);
let hf_denominator = qu(df_c) * (qu(df_s) - qu(df_c) * &epsilon_gg);
let epsilon_hf = hf_denominator
.is_positive()
.then(|| (qu(n) * qu(df_c) * &epsilon_gg - qi(2)) / hf_denominator);
let corrected = |eps: &Q| f_sf_rational(ctx, &(qu(df_c) * eps), &(qu(df_e) * eps), &f);
let p_value_gg = corrected(&epsilon_gg);
let one = Q::one();
let p_value_hf = epsilon_hf
.as_ref()
.map(|e| corrected(if *e > one { &one } else { e }));
let mauchly = if k >= 3 {
mauchly(ctx, &s, &trace, n, k)?
} else {
None
};
Ok(RepeatedMeasuresAnova {
n_subjects: n,
n_conditions: k,
conditions,
subjects: AnovaRow::untested(Source::Subjects, ss_subjects, df_s),
error: AnovaRow::untested(Source::Residual, ss_error, df_e),
total: AnovaRow::total(ss_total, n * k - 1),
f,
p_value,
epsilon_gg,
epsilon_hf,
p_value_gg,
p_value_hf,
mauchly,
grand_mean,
condition_means,
subject_means,
})
}
struct GaussLegendre {
nodes: Vec<f64>,
weights: Vec<f64>,
}
impl GaussLegendre {
fn new(m: usize) -> Self {
let mut nodes = Vec::with_capacity(m);
let mut weights = Vec::with_capacity(m);
for i in 0..m {
let mut x = (PI * (i as f64 + 0.75) / (m as f64 + 0.5)).cos();
for _ in 0..100 {
let (p, d) = legendre(m, x);
let step = p / d;
x -= step;
if step.abs() < 1e-15 {
break;
}
}
let (_, dp) = legendre(m, x);
nodes.push(x);
weights.push(2.0 / ((1.0 - x * x) * dp * dp));
}
Self { nodes, weights }
}
fn integrate(&self, f: &dyn Fn(f64) -> f64, a: f64, b: f64, panels: usize) -> f64 {
let panels = panels.max(1);
let h = (b - a) / panels as f64;
let half = 0.5 * h;
let mut total = 0.0;
for p in 0..panels {
let mid = a + (p as f64 + 0.5) * h;
let mut acc = 0.0;
for (x, w) in self.nodes.iter().zip(&self.weights) {
acc += w * f(mid + half * x);
}
total += half * acc;
}
total
}
}
fn legendre(m: usize, x: f64) -> (f64, f64) {
let mut p0 = 1.0;
let mut p1 = x;
for j in 2..=m {
let jf = j as f64;
let p2 = ((2.0 * jf - 1.0) * x * p1 - (jf - 1.0) * p0) / jf;
p0 = p1;
p1 = p2;
}
let dp = m as f64 * (x * p1 - p0) / (x * x - 1.0);
(p1, dp)
}
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum Tail {
Lower,
Upper,
}
impl Tail {
fn at_zero(self) -> f64 {
match self {
Tail::Lower => 0.0,
Tail::Upper => 1.0,
}
}
fn at_infinity(self) -> f64 {
1.0 - self.at_zero()
}
}
fn norm_interval(a: f64, w: f64) -> f64 {
let h = 0.5 * w;
let m = a + h;
if h * (m.abs() + 1.0) <= 0.1 {
let (mut even, mut odd) = (1.0, m); let (mut sum, mut hpow, mut fact) = (1.0, 1.0, 1.0);
for j in 1..=12 {
let n = (2 * j) as f64;
even = m * odd - (n - 1.0) * even; odd = m * even - n * odd; hpow *= h * h;
fact *= (n - 1.0) * n; let term = even * hpow / ((n + 1.0) * fact);
sum += term;
if term.abs() <= 1e-17 * sum.abs() {
break;
}
}
return norm_pdf(m) * w * sum;
}
let b = a + w;
if a >= 0.0 {
0.5 * (erfc(a / SQRT_2) - erfc(b / SQRT_2))
} else if b <= 0.0 {
0.5 * (erfc(-b / SQRT_2) - erfc(-a / SQRT_2))
} else {
0.5 * (erf(b / SQRT_2) + erf(-a / SQRT_2))
}
}
fn min_normal_cut(k: usize) -> f64 {
(2.0 * ((k as f64).ln() + 39.2)).sqrt()
}
fn normal_range_tail(w: f64, k: usize, tail: Tail, rule: &GaussLegendre) -> f64 {
if w.is_nan() || w <= 0.0 {
return tail.at_zero();
}
if w == f64::INFINITY {
return tail.at_infinity();
}
let kf = k as f64;
let km1 = kf - 1.0;
let cut = min_normal_cut(k);
let integral = match tail {
Tail::Lower => {
let f = |z: f64| norm_pdf(z) * norm_interval(z, w).powf(km1);
let (a, b, width) = (-cut, 8.5, (2.8 / kf.sqrt()).min(2.0));
rule.integrate(&f, a, b, ((b - a) / width).ceil() as usize)
}
Tail::Upper => {
let f = |z: f64| {
let above = norm_sf(z);
if above <= 0.0 {
return 0.0;
}
let r = (norm_sf(z + w) / above).min(1.0);
let one_exceeds = -(km1 * (-r).ln_1p()).exp_m1();
norm_pdf(z) * above.powf(km1) * one_exceeds
};
let centre = -0.5 * w;
let (a, b) = ((centre - 7.5).min(-cut), centre + 7.5);
let width = (2.8 / (2.0 * kf.ln()).max(1.0).sqrt()).min(2.0);
rule.integrate(&f, a, b, ((b - a) / width).ceil() as usize)
}
};
(kf * integral).clamp(0.0, 1.0)
}
const TABLE_DEGREE: usize = 20;
const TABLE_PIECE: f64 = 1.5;
const TABLE_EXP_CAP: f64 = 700.0;
struct RangeUpperTable<'r> {
k: usize,
rule: &'r GaussLegendre,
width: f64,
zero_from: f64,
points: Vec<f64>,
weights: Vec<f64>,
pieces: std::cell::RefCell<Vec<Option<Option<Vec<f64>>>>>,
}
impl<'r> RangeUpperTable<'r> {
fn new(k: usize, rule: &'r GaussLegendre) -> Self {
let n = TABLE_DEGREE;
let kf = k as f64;
let width = TABLE_PIECE / (2.0 * kf.ln()).max(1.0).sqrt();
let zero_from = 2.0 * ((0.5 * kf * (kf - 1.0)).ln() + 1075.0 * LN_2).sqrt();
let points = (0..=n).map(|j| (PI * j as f64 / n as f64).cos()).collect();
let weights = (0..=n)
.map(|j| {
let sign = if j % 2 == 0 { 1.0 } else { -1.0 };
if j == 0 || j == n { 0.5 * sign } else { sign }
})
.collect();
let count = (zero_from / width).ceil() as usize;
Self {
k,
rule,
width,
zero_from,
points,
weights,
pieces: std::cell::RefCell::new(vec![None; count]),
}
}
fn square(w: f64) -> (f64, f64) {
let hi = w * w;
(hi, w.mul_add(w, -hi))
}
fn sample(&self, v: f64) -> f64 {
let (hi, lo) = Self::square(v);
let r = self.direct(v);
if 0.25 * hi < TABLE_EXP_CAP {
(r * (0.25 * hi).exp()).ln() + 0.25 * lo
} else {
r.ln() + 0.25 * hi + 0.25 * lo
}
}
fn direct(&self, w: f64) -> f64 {
normal_range_tail(w, self.k, Tail::Upper, self.rule)
}
fn upper(&self, w: f64) -> f64 {
if w >= self.zero_from {
return 0.0;
}
if w.is_nan() || w <= 0.0 {
return self.direct(w);
}
let pieces_len = self.pieces.borrow().len();
let index = ((w / self.width) as usize).min(pieces_len - 1);
let a = index as f64 * self.width;
if self.pieces.borrow()[index].is_none() {
let values: Vec<f64> = self
.points
.iter()
.map(|&x| self.sample(a + 0.5 * self.width * (1.0 + x)))
.collect();
let usable = values.iter().all(|f| f.is_finite());
self.pieces.borrow_mut()[index] = Some(usable.then_some(values));
}
let pieces = self.pieces.borrow();
let Some(Some(values)) = &pieces[index] else {
return self.direct(w);
};
let x = 2.0 * (w - a) / self.width - 1.0;
let (mut num, mut den) = (0.0, 0.0);
for ((&p, &c), &f) in self.points.iter().zip(&self.weights).zip(values) {
let d = x - p;
if d == 0.0 {
(num, den) = (f, 1.0);
break;
}
num += c * f / d;
den += c / d;
}
let (hi, lo) = Self::square(w);
let f = num / den - 0.25 * lo;
let r = if 0.25 * hi < TABLE_EXP_CAP {
f.exp() * (-0.25 * hi).exp()
} else {
(f - 0.25 * hi).exp()
};
r.clamp(0.0, 1.0)
}
}
fn stirling_remainder(x: f64) -> f64 {
const B: [f64; 7] = [
1.0 / 12.0,
-1.0 / 360.0,
1.0 / 1260.0,
-1.0 / 1680.0,
1.0 / 1188.0,
-691.0 / 360_360.0,
1.0 / 156.0,
];
let mut x = x;
let mut shift = 0.0;
while x < 10.0 {
let v2 = (1.0 / (2.0 * x + 1.0)).powi(2);
let (mut pow, mut sum) = (v2, 0.0);
for j in 1..=40 {
let term = pow / (2 * j + 1) as f64;
sum += term;
if term <= 1e-17 * sum {
break;
}
pow *= v2;
}
shift += sum;
x += 1.0;
}
let inv2 = 1.0 / (x * x);
shift + B.iter().rev().fold(0.0, |acc, &b| acc * inv2 + b) / x
}
fn ln1p_minus_u(u: f64) -> f64 {
if u.abs() >= 0.25 {
return u.ln_1p() - u;
}
let mut term = u * u;
let mut sum = 0.0;
for n in 2..=40 {
let contribution = term / n as f64;
sum += if n % 2 == 0 {
-contribution
} else {
contribution
};
term *= u;
if contribution.abs() <= 1e-18 * sum.abs() {
break;
}
}
sum
}
const NU_NORMAL_LIMIT: f64 = 1e30;
const OUTER_LOG_DROP: f64 = 45.0;
const OUTER_PANEL_LOG_CHANGE: f64 = 12.0;
const UNDERFLOW_LOG_MODE: f64 = -760.0;
const OUTER_SEARCH_CAP: usize = 4096;
fn studentized_range_tail(q: f64, k: usize, nu: f64, tail: Tail) -> Result<f64, SymplexError> {
let rule = GaussLegendre::new(16);
studentized_range_tail_by(q, k, nu, tail, &|w| normal_range_tail(w, k, tail, &rule))
}
fn studentized_range_tail_by(
q: f64,
k: usize,
nu: f64,
tail: Tail,
range_tail_at: &dyn Fn(f64) -> f64,
) -> Result<f64, SymplexError> {
const OP: &str = "studentized_range";
if q.is_nan() || q <= 0.0 {
return Ok(tail.at_zero());
}
if q == f64::INFINITY {
return Ok(tail.at_infinity());
}
let rule = GaussLegendre::new(16);
if nu >= NU_NORMAL_LIMIT {
return Ok(range_tail_at(q));
}
let x = 0.5 * nu;
let log_prefactor = LN_2 + 0.5 * (x / (2.0 * PI)).ln() - stirling_remainder(x);
let log_weight = |t: f64| {
let u = t.exp_m1();
let t_minus_u = if u.abs() < 0.25 {
ln1p_minus_u(u)
} else {
t - u
};
log_prefactor + 2.0 * x * t_minus_u - x * u * u
};
let range_tail = |t: f64| range_tail_at(q * t.exp());
let log_h = |t: f64| log_weight(t) + range_tail(t).ln();
let not_found = || {
failed(
OP,
format!(
"the studentized range quadrature found no window (q = {q}, k = {k}, df = {nu})"
),
)
};
let scale = 1.0 / (2.0 * nu + 1.0).sqrt();
let (t0, default_dir) = match tail {
Tail::Upper => {
let r = q / (2.0 * nu).sqrt();
let t0 = if r > 1e100 {
-r.ln()
} else {
-0.5 * (r * r).ln_1p()
};
(t0, -1.0)
}
Tail::Lower => (0.5 * ((k - 1) as f64 / nu).ln_1p(), 1.0),
};
let l0 = log_h(t0);
let (l_right, l_left) = (log_h(t0 + scale), log_h(t0 - scale));
let dir = if l_right > l0 && l_right >= l_left {
1.0
} else if l_left > l0 {
-1.0
} else if l0 == f64::NEG_INFINITY {
default_dir
} else {
0.0
};
let (mut lo, mut mode, mut hi, mut l_mode) = (t0 - scale, t0, t0 + scale, l0);
if dir != 0.0 {
let l_first = if dir > 0.0 { l_right } else { l_left };
let (mut prev, mut cur, mut l_cur) = (t0, t0 + dir * scale, l_first);
let mut step = scale;
let mut bracketed = false;
for _ in 0..OUTER_SEARCH_CAP {
step *= 2.0;
let next = cur + dir * step;
if next.is_nan() || next.abs() >= 1500.0 {
if l_cur == f64::NEG_INFINITY {
return Ok(0.0);
}
break;
}
let l_next = log_h(next);
if l_next > l_cur || (l_next == f64::NEG_INFINITY && l_cur == f64::NEG_INFINITY) {
(prev, cur, l_cur) = (cur, next, l_next);
} else {
(lo, hi) = if dir > 0.0 {
(prev, next)
} else {
(next, prev)
};
(mode, l_mode) = (cur, l_cur);
bracketed = true;
break;
}
}
if !bracketed || !l_mode.is_finite() {
return Err(not_found());
}
}
if l_mode < UNDERFLOW_LOG_MODE {
return Ok(0.0);
}
let tol = 0.05 * scale / (k as f64).sqrt();
for _ in 0..OUTER_SEARCH_CAP {
if hi - lo <= tol {
break;
}
const GOLDEN: f64 = 0.381_966_011_250_105_1;
let left = mode - lo > hi - mode;
let probe = if left {
mode - GOLDEN * (mode - lo)
} else {
mode + GOLDEN * (hi - mode)
};
if probe == mode {
break;
}
let l_probe = log_h(probe);
match (l_probe > l_mode, left) {
(true, true) => (hi, mode, l_mode) = (mode, probe, l_probe),
(true, false) => (lo, mode, l_mode) = (mode, probe, l_probe),
(false, true) => lo = probe,
(false, false) => hi = probe,
}
}
let h = |t: f64| log_weight(t).exp() * range_tail(t);
let floor = l_mode - OUTER_LOG_DROP;
let first = 0.25 * scale / (k as f64).sqrt();
let mut total = 0.0;
for dir in [-1.0, 1.0] {
let (mut pos, mut l_pos, mut step) = (mode, l_mode, first);
let mut reached = false;
'march: for _ in 0..OUTER_SEARCH_CAP {
let (mut next, mut l_next) = (pos + dir * step, f64::NAN);
for _ in 0..OUTER_SEARCH_CAP {
next = pos + dir * step;
if !next.is_finite() || next == pos {
break 'march;
}
l_next = log_h(next);
if (l_next - l_pos).abs() <= OUTER_PANEL_LOG_CHANGE || l_next == f64::NEG_INFINITY {
break;
}
step *= 0.5;
}
total += rule.integrate(&h, pos.min(next), pos.max(next), 1);
if l_next < floor {
reached = true;
break;
}
(pos, l_pos, step) = (next, l_next, 2.0 * step);
}
if !reached {
return Err(not_found());
}
}
Ok(total.clamp(0.0, 1.0))
}
fn check_studentized_range_args(op: &'static str, k: usize, df: f64) -> Result<(), SymplexError> {
if k < 2 {
return Err(invalid(op, format!("k must be at least 2, got {k}")));
}
if df.is_nan() || df < 1.0 {
return Err(invalid(
op,
format!("the degrees of freedom must be at least 1 (or +inf), got {df}"),
));
}
Ok(())
}
pub fn studentized_range_cdf(q: f64, k: usize, df: f64) -> Result<f64, SymplexError> {
const OP: &str = "studentized_range_cdf";
check_studentized_range_args(OP, k, df)?;
if q.is_nan() {
return Err(invalid(OP, "q must not be NaN"));
}
studentized_range_tail(q, k, df, Tail::Lower)
}
pub fn studentized_range_sf(q: f64, k: usize, df: f64) -> Result<f64, SymplexError> {
const OP: &str = "studentized_range_sf";
check_studentized_range_args(OP, k, df)?;
if q.is_nan() {
return Err(invalid(OP, "q must not be NaN"));
}
studentized_range_tail(q, k, df, Tail::Upper)
}
pub fn studentized_range_quantile(p: f64, k: usize, df: f64) -> Result<f64, SymplexError> {
check_studentized_range_args("studentized_range_quantile", k, df)?;
check_unit_open("studentized_range_quantile", "p", p)?;
studentized_range_root(p, k, df, None)
}
fn studentized_range_root(
p: f64,
k: usize,
df: f64,
table: Option<&RangeUpperTable<'_>>,
) -> Result<f64, SymplexError> {
const OP: &str = "studentized_range_quantile";
let (tail, level) = if p <= 0.5 {
(Tail::Lower, p)
} else {
(Tail::Upper, 1.0 - p)
};
let rule = GaussLegendre::new(16);
let range_tail_at = |w: f64| match table {
Some(t) if tail == Tail::Upper => t.upper(w),
_ => normal_range_tail(w, k, tail, &rule),
};
let ln_level = level.ln();
let g = |x: f64| match studentized_range_tail_by(x.exp(), k, df, tail, &range_tail_at) {
Ok(v) => v.ln().max(-1e4) - ln_level,
Err(_) => f64::NAN,
};
let bracket = grow_bracket(g, 0.0, 1.0, Bounds::free(), 64).map_err(|e| {
failed(
OP,
format!("no bracket for the studentized range quantile: {e}"),
)
})?;
let opts = RootOpts {
xtol: 1e-14,
..RootOpts::default()
};
let root = brent_root(g, bracket.lower, bracket.upper, &opts)
.map_err(|e| failed(OP, e.to_string()))?
.exp();
if root.is_finite() && root > 0.0 {
Ok(root)
} else {
Err(failed(
OP,
"the studentized range quantile did not converge",
))
}
}
#[derive(Clone, Debug, PartialEq)]
pub struct PairwiseComparison {
pub i: usize,
pub j: usize,
pub diff: Q,
pub se: Ex,
pub statistic: Ex,
pub p_adj: f64,
pub ci: Interval<f64>,
}
pub fn tukey_hsd(
ctx: &Context,
groups: &[Vec<Q>],
confidence: f64,
) -> Result<Vec<PairwiseComparison>, SymplexError> {
const OP: &str = "tukey_hsd";
check_confidence(OP, confidence)?;
let k = groups.len();
if k < 2 {
return Err(invalid(OP, "at least two groups are needed"));
}
if groups.iter().any(Vec::is_empty) {
return Err(invalid(OP, "every group must be non-empty"));
}
let n: usize = groups.iter().map(Vec::len).sum();
if n <= k {
return Err(invalid(
OP,
"at least one group needs more than one observation",
));
}
let df = n - k;
let means = groups
.iter()
.map(|g| data::mean(g))
.collect::<Result<Vec<_>, _>>()?;
let ss_within = groups.iter().fold(Q::zero(), |acc, g| acc + centred_ss(g));
if ss_within.is_zero() {
return Err(invalid(OP, "the within-group variance is zero"));
}
let mse = ss_within / qu(df);
let df_f = df as f64;
let rule = GaussLegendre::new(16);
let table = RangeUpperTable::new(k, &rule);
let upper = |w: f64| table.upper(w);
let q_crit = studentized_range_root(confidence, k, df_f, Some(&table))?;
let mut out = Vec::with_capacity(k * (k - 1) / 2);
for i in 0..k {
for j in i + 1..k {
let diff = &means[i] - &means[j];
let var = &mse / qi(2) * (qu(groups[i].len()).recip() + qu(groups[j].len()).recip());
let se = ex(ctx, &var).sqrt();
let statistic = ex(ctx, &diff.abs()) / &se;
let (diff_f, se_f) = (to_f64(OP, &diff)?, to_f64(OP, &var)?.sqrt());
let stat_f = to_f64(OP, &(&diff * &diff / &var))?.sqrt();
let p_adj = studentized_range_tail_by(stat_f, k, df_f, Tail::Upper, &upper)?;
out.push(PairwiseComparison {
i,
j,
diff,
se,
statistic,
p_adj,
ci: Interval::closed(diff_f - q_crit * se_f, diff_f + q_crit * se_f),
});
}
}
Ok(out)
}
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum Adjustment {
Bonferroni,
Holm,
}
#[derive(Clone, Debug, PartialEq)]
pub struct PairwiseTTest {
pub i: usize,
pub j: usize,
pub diff: Q,
pub test: TestResult,
pub p_adj: f64,
pub reject: bool,
}
pub fn pairwise_t_tests(
ctx: &Context,
groups: &[Vec<Q>],
adjustment: Adjustment,
alpha: f64,
) -> Result<Vec<PairwiseTTest>, SymplexError> {
const OP: &str = "pairwise_t_tests";
let k = groups.len();
if k < 2 {
return Err(invalid(OP, "at least two groups are needed"));
}
let mut pairs = Vec::with_capacity(k * (k - 1) / 2);
let mut p_values = Vec::with_capacity(k * (k - 1) / 2);
for i in 0..k {
for j in i + 1..k {
let test = hypothesis::t_test_two_sample(
ctx,
&groups[i],
&groups[j],
false,
Alternative::TwoSided,
)?;
p_values.push(test.p_value_f64()?);
let diff = data::mean(&groups[i])? - data::mean(&groups[j])?;
pairs.push((i, j, diff, test));
}
}
let adjusted = match adjustment {
Adjustment::Bonferroni => hypothesis::bonferroni(&p_values, alpha)?,
Adjustment::Holm => hypothesis::holm(&p_values, alpha)?,
};
Ok(pairs
.into_iter()
.zip(adjusted.p_adjusted)
.zip(adjusted.reject)
.map(|(((i, j, diff, test), p_adj), reject)| PairwiseTTest {
i,
j,
diff,
test,
p_adj,
reject,
})
.collect())
}