symplex 0.14.0

//! Hypothesis tests, effect sizes, multiple-comparison corrections,
//! resampling and power / sample-size utilities.
//!
//! The tests take exact observations ([`Q`], `Ratio<BigInt>`, see
//! [`super::data`]) and follow one principle: **exact where possible,
//! honest where not**.
//!
//! * A test **statistic** is exact: a rational (`U`, `H`, `χ²`, `F`, an odds
//!   ratio) or a rational times a square root (`t`, `z`, Cohen's `d`), as an
//!   [`Ex`].
//! * A **p-value** is an exact *expression* of the statistic — the Student-t
//!   tail through `betainc_regularized`, a χ² tail through `uppergamma`,
//!   a normal tail through `erfc`, or an exact rational for the discrete
//!   exact tests (binomial, Fisher, McNemar, sign, exact Mann–Whitney /
//!   Wilcoxon / Kendall) — that the caller evaluates with
//!   [`Ex::eval_f64`] (or [`TestResult::p_value_f64`]).  Nothing is rounded
//!   before the caller asks.
//! * The few quantities that are inherently numerical (confidence limits
//!   through a quantile, the Kolmogorov–Smirnov `D` of a transcendental
//!   CDF, bootstrap / permutation p-values, power) are `f64` and say so.
//!
//! Every function names the `scipy.stats` / `statsmodels` routine whose
//! conventions it follows (alternative hypotheses, tie corrections,
//! continuity corrections, two-sided definitions of the discrete tests);
//! the tests in `tests/v13/v13_hypothesis.rs` pin the agreement to ≥ 1e-9.
//!
//! ```
//! use symplex::prelude::*;
//! use symplex::stats::data::from_i64;
//! use symplex::stats::hypothesis::{t_test_one_sample, Alternative};
//!
//! let ctx = Context::new();
//! let x = from_i64(&[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]);
//! let r = t_test_one_sample(&ctx, &x, &symplex::linprog::qi(5), Alternative::TwoSided)?;
//! // scipy: ttest_1samp(range(1, 11), 5) → statistic 0.5222329678670935, pvalue 0.6141172548083939
//! assert!((r.statistic_f64()? - 0.522_232_967_867_093_5).abs() < 1e-12);
//! assert!((r.p_value_f64()? - 0.614_117_254_808_393_9).abs() < 1e-12);
//! assert_eq!(r.df, Some(ctx.int(9)));
//! # Ok::<(), SymplexError>(())
//! ```

use std::cmp::Ordering;

use num_bigint::BigInt;
use num_traits::{One, Signed, ToPrimitive, Zero};

use super::data::{self, Ddof, Q};
use super::family::Distribution;
use super::sample::Rng;
use crate::api::context::Context;
use crate::api::expr::Ex;
use crate::base::errors::SymplexError;
use crate::calculus::definite::{QuadOpts, quadrature};
use crate::domains::optimize::{RootOpts, brent_root};
use crate::output::codegen::numeric_rt::{erfc, erfcinv, lgamma};

// ═══════════════════════════════════════════════════════════════════════════
// Small helpers
// ═══════════════════════════════════════════════════════════════════════════

fn invalid(op: &'static str, reason: impl Into<String>) -> SymplexError {
    SymplexError::invalid_argument(op, reason)
}

fn qi(n: i64) -> Q {
    Q::from_integer(BigInt::from(n))
}

fn qu(n: usize) -> Q {
    Q::from_integer(BigInt::from(n))
}

fn ex(ctx: &Context, q: &Q) -> Ex {
    ctx.from_ratio(q.clone())
}

fn q_to_f64(q: &Q) -> f64 {
    data::to_f64(std::slice::from_ref(q))
        .first()
        .copied()
        .unwrap_or(f64::NAN)
}

fn usize_to_i64(op: &'static str, n: usize) -> Result<i64, SymplexError> {
    i64::try_from(n).map_err(|_| invalid(op, format!("{n} does not fit in an i64")))
}

fn factorial_big(n: usize) -> BigInt {
    (1..=n).fold(BigInt::one(), |acc, k| acc * BigInt::from(k))
}

fn sum_big(v: &[BigInt]) -> BigInt {
    v.iter().fold(BigInt::zero(), |acc, x| acc + x)
}

fn check_finite(op: &'static str, name: &str, v: f64) -> Result<(), SymplexError> {
    if v.is_finite() {
        Ok(())
    } else {
        Err(invalid(op, format!("{name} must be finite, got {v}")))
    }
}

fn check_unit_open(op: &'static str, name: &str, v: f64) -> Result<(), SymplexError> {
    if v > 0.0 && v < 1.0 {
        Ok(())
    } else {
        Err(invalid(
            op,
            format!("{name} must lie strictly between 0 and 1, got {v}"),
        ))
    }
}

fn check_sample(op: &'static str, name: &str, x: &[Q], min: usize) -> Result<(), SymplexError> {
    if x.len() < min {
        return Err(invalid(
            op,
            format!(
                "{name} needs at least {min} observation{}, got {}",
                if min == 1 { "" } else { "s" },
                x.len()
            ),
        ));
    }
    Ok(())
}

fn check_same_len(op: &'static str, x: &[Q], y: &[Q]) -> Result<(), SymplexError> {
    if x.len() != y.len() {
        return Err(invalid(
            op,
            format!(
                "paired samples must have the same size ({} and {})",
                x.len(),
                y.len()
            ),
        ));
    }
    Ok(())
}

/// Φ(x) in `f64`.
fn norm_cdf(x: f64) -> f64 {
    0.5 * erfc(-x / std::f64::consts::SQRT_2)
}

/// 1 − Φ(x) in `f64`.
fn norm_sf(x: f64) -> f64 {
    0.5 * erfc(x / std::f64::consts::SQRT_2)
}

/// Φ⁻¹(1 − α) in `f64`.
fn norm_isf(alpha: f64) -> f64 {
    std::f64::consts::SQRT_2 * erfcinv(2.0 * alpha)
}

/// `P(T ≤ t)` for Student's t with `ν` degrees of freedom, as an `f64`:
/// the exact expression `½ I_{ν/(t²+ν)}(ν/2, ½)` evaluated by `eval_f64`
/// (arbitrary precision, then rounded).
fn student_t_cdf_f64(ctx: &Context, df: f64, t: f64) -> Result<f64, SymplexError> {
    let nu = ctx.from_f64(df)?;
    let z = &nu / (ctx.from_f64(t * t)? + &nu);
    let tail = ctx.rational(1, 2)
        * z.betainc_regularized(&(&nu / ctx.int(2)), &ctx.rational(1, 2), &ctx.zero());
    let tail = tail.eval_f64()?;
    Ok(if t < 0.0 { tail } else { 1.0 - tail })
}

/// The Student-t quantile `t_{p, ν}` by Brent's method on
/// [`student_t_cdf_f64`] over a bracket grown from `0` (the distribution
/// is symmetric, so `p < ½` is mirrored).
fn student_t_quantile_f64(
    op: &'static str,
    ctx: &Context,
    df: f64,
    p: f64,
) -> Result<f64, SymplexError> {
    if p < 0.5 {
        return student_t_quantile_f64(op, ctx, df, 1.0 - p).map(|t| -t);
    }
    if p == 0.5 {
        return Ok(0.0);
    }
    let g = |t: f64| student_t_cdf_f64(ctx, df, t).unwrap_or(f64::NAN) - p;
    let mut hi = 1.0;
    for _ in 0..64 {
        let v = g(hi);
        if v.is_nan() {
            return Err(SymplexError::computation_failed(
                op,
                "the Student-t distribution function could not be evaluated",
            ));
        }
        if v >= 0.0 {
            break;
        }
        hi *= 2.0;
    }
    let root = brent_root(g, 0.0, hi, &RootOpts::default())
        .map_err(|e| SymplexError::computation_failed(op, e.to_string()))?;
    if root.is_finite() {
        Ok(root)
    } else {
        Err(SymplexError::computation_failed(
            op,
            "the Student-t quantile did not converge",
        ))
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// Result types
// ═══════════════════════════════════════════════════════════════════════════

/// The alternative hypothesis of a test (`scipy`'s `alternative=`).
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum Alternative {
    /// `'two-sided'`: the parameter differs from its null value.
    TwoSided,
    /// `'less'`: the parameter is smaller than its null value.
    Less,
    /// `'greater'`: the parameter is larger than its null value.
    Greater,
}

/// The outcome of a test: an exact statistic, an exact p-value expression
/// and (when the reference distribution has one) the degrees of freedom.
#[derive(Clone, Debug, PartialEq)]
pub struct TestResult {
    /// The test statistic, exact (a rational, or a rational times a root).
    pub statistic: Ex,
    /// The p-value as an exact expression; evaluate with
    /// [`p_value_f64`](Self::p_value_f64).
    pub p_value: Ex,
    /// Degrees of freedom of the reference distribution, if any (exact; the
    /// Welch–Satterthwaite `ν` is a rational).
    pub df: Option<Ex>,
    /// The alternative the p-value refers to.
    pub alternative: Alternative,
}

impl TestResult {
    /// The p-value as an `f64`.
    ///
    /// # Errors
    ///
    /// Propagates the evaluation error of the expression (not expected for
    /// the expressions this module builds).
    pub fn p_value_f64(&self) -> Result<f64, SymplexError> {
        self.p_value.eval_f64()
    }

    /// The statistic as an `f64`.
    ///
    /// # Errors
    ///
    /// As [`p_value_f64`](Self::p_value_f64).
    pub fn statistic_f64(&self) -> Result<f64, SymplexError> {
        self.statistic.eval_f64()
    }

    /// The p-value as an exact rational when it is one (the discrete exact
    /// tests: binomial, Fisher, McNemar, sign, exact rank tests).
    #[must_use]
    pub fn p_value_exact(&self) -> Option<Q> {
        self.p_value.as_rational()
    }

    /// The statistic as an exact rational when it is one (`U`, `H`, `χ²`,
    /// `F`, counts and proportions).
    #[must_use]
    pub fn statistic_exact(&self) -> Option<Q> {
        self.statistic.as_rational()
    }
}

/// How the p-value of a rank test is computed.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum RankMethod {
    /// The exact permutation distribution of the statistic (an exact
    /// rational p-value).  Only valid without ties.
    Exact,
    /// The normal approximation with the tie correction, optionally with the
    /// continuity correction (`scipy`'s `use_continuity` / `correction`).
    Asymptotic {
        /// Move the statistic half a unit towards its null mean.
        continuity: bool,
    },
}

/// The outcome of a Pearson χ² test.
#[derive(Clone, Debug, PartialEq)]
pub struct ChiSquareResult {
    /// The χ² statistic, exact.
    pub statistic: Q,
    /// Degrees of freedom.
    pub df: usize,
    /// `P(χ²_df ≥ statistic)` as an exact expression (`uppergamma(df/2,
    /// statistic/2) / Γ(df/2)`).
    pub p_value: Ex,
    /// The expected counts under the null, exact (one row for a
    /// goodness-of-fit test).
    pub expected: Vec<Vec<Q>>,
}

impl ChiSquareResult {
    /// The p-value as an `f64`.
    ///
    /// # Errors
    ///
    /// Propagates the evaluation error of the expression.
    pub fn p_value_f64(&self) -> Result<f64, SymplexError> {
        self.p_value.eval_f64()
    }
}

/// The outcome of a one-way analysis of variance.
#[derive(Clone, Debug, PartialEq)]
pub struct AnovaResult {
    /// The `F` statistic `(SS_between/df_between) / (SS_within/df_within)`,
    /// exact.
    pub f: Q,
    /// `k − 1`.
    pub df_between: usize,
    /// `N − k`.
    pub df_within: usize,
    /// `P(F_{df_between, df_within} ≥ f)` as an exact expression.
    pub p_value: Ex,
    /// `Σ nᵢ (x̄ᵢ − x̄)²`.
    pub ss_between: Q,
    /// `Σᵢ Σⱼ (xᵢⱼ − x̄ᵢ)²`.
    pub ss_within: Q,
    /// `SS_between / SS_total`, the proportion of variance explained.
    pub eta_squared: Q,
}

impl AnovaResult {
    /// The p-value as an `f64`.
    ///
    /// # Errors
    ///
    /// Propagates the evaluation error of the expression.
    pub fn p_value_f64(&self) -> Result<f64, SymplexError> {
        self.p_value.eval_f64()
    }
}

/// A ratio estimate (odds ratio, relative risk) with a Wald confidence
/// interval on the log scale.
#[derive(Clone, Debug, PartialEq)]
pub struct RatioEstimate {
    /// The point estimate, exact.
    pub estimate: Q,
    /// `(lower, upper)` of the `confidence` Wald interval, computed on the
    /// log scale and exponentiated.
    pub ci: (f64, f64),
    /// The confidence level of `ci`.
    pub confidence: f64,
}

/// The outcome of a Kolmogorov–Smirnov test (numerical: the statistic
/// compares an empirical CDF against a transcendental one).
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct KsResult {
    /// `D` (two-sided), `D⁺` (`Greater`) or `D⁻` (`Less`).
    pub statistic: f64,
    /// The p-value: Kolmogorov's asymptotic distribution two-sided,
    /// Smirnov's exact one-sided distribution otherwise.
    pub p_value: f64,
    /// The alternative the p-value refers to.
    pub alternative: Alternative,
}

/// Adjusted p-values and rejection flags of a multiple-comparison
/// procedure, in the order of the input.
#[derive(Clone, Debug, PartialEq)]
pub struct Adjusted {
    /// The adjusted p-values (clipped to `1`).
    pub p_adjusted: Vec<f64>,
    /// Whether each hypothesis is rejected at the family-wise (or false
    /// discovery) level `alpha`.
    pub reject: Vec<bool>,
}

/// How a bootstrap confidence interval is read off the resampled
/// statistics.
#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
pub enum BootstrapMethod {
    /// The `α/2` and `1 − α/2` quantiles of the bootstrap distribution.
    Percentile,
    /// The reflected ("basic", "reverse percentile") interval
    /// `(2θ̂ − q_{1−α/2}, 2θ̂ − q_{α/2})`.
    Basic,
}

/// The outcome of a randomised permutation test.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct PermutationResult {
    /// The observed statistic.
    pub statistic: f64,
    /// `(#{permuted statistics at least as extreme} + 1) / (n_permutations + 1)`.
    pub p_value: f64,
}

// ═══════════════════════════════════════════════════════════════════════════
// Exact tails of the reference distributions
// ═══════════════════════════════════════════════════════════════════════════

/// A statistic `num / √var` with `num`, `var` rational and `var > 0`: a
/// `t` or `z` whose square is rational and whose sign is known exactly.
struct RootRatio {
    num: Q,
    var: Q,
}

impl RootRatio {
    fn to_ex(&self, ctx: &Context) -> Ex {
        (ex(ctx, &self.num) / ex(ctx, &self.var).sqrt()).simplify()
    }

    fn square(&self) -> Q {
        &self.num * &self.num / &self.var
    }

    fn in_tail(&self, alt: Alternative) -> bool {
        match alt {
            Alternative::Greater | Alternative::TwoSided => !self.num.is_negative(),
            Alternative::Less => !self.num.is_positive(),
        }
    }
}

/// The p-value of a Student-t statistic with `df` degrees of freedom.
/// With `z = ν/(t² + ν)` (rational), `P(|T| ≥ |t|) = I_z(ν/2, ½)` and
/// `P(T ≥ t) = ½ I_z(ν/2, ½)` for `t ≥ 0`.
fn student_p_value(ctx: &Context, df: &Q, stat: &RootRatio, alt: Alternative) -> Ex {
    let t2 = stat.square();
    let z = df / (&t2 + df);
    let two_sided =
        ex(ctx, &z).betainc_regularized(&ex(ctx, &(df / qi(2))), &ctx.rational(1, 2), &ctx.zero());
    one_sided_from_symmetric(ctx, two_sided, stat.in_tail(alt), alt)
}

/// The p-value of a standard-normal statistic: `P(|Z| ≥ |z|) = erfc(|z|/√2)`,
/// `P(Z ≥ z) = ½ erfc(z/√2)`.
fn normal_p_value(ctx: &Context, stat: &RootRatio, alt: Alternative) -> Ex {
    let two_sided = ex(ctx, &(stat.square() / qi(2))).sqrt().erfc();
    one_sided_from_symmetric(ctx, two_sided, stat.in_tail(alt), alt)
}

/// Turn the two-sided tail `P(|X| ≥ |x|)` of a symmetric distribution into
/// the requested one: half of it when `x` lies in the alternative's tail,
/// its complement otherwise.
fn one_sided_from_symmetric(ctx: &Context, two_sided: Ex, in_tail: bool, alt: Alternative) -> Ex {
    match alt {
        Alternative::TwoSided => two_sided,
        Alternative::Greater | Alternative::Less => {
            let half = ctx.rational(1, 2) * two_sided;
            if in_tail { half } else { ctx.one() - half }
        }
    }
}

/// `P(χ²_df ≥ x) = Γ(df/2, x/2) / Γ(df/2)` as an expression.
fn chi_squared_sf(ctx: &Context, df: usize, x: &Ex) -> Ex {
    let half_df = ctx.rational(df as i64, 2);
    (x / ctx.int(2)).uppergamma(&half_df) / half_df.gamma()
}

fn chi_squared_sf_q(ctx: &Context, df: usize, x: &Q) -> Ex {
    if !x.is_positive() {
        return ctx.one();
    }
    chi_squared_sf(ctx, df, &ex(ctx, x))
}

/// `P(F_{d₁,d₂} ≥ f) = I_{d₂/(d₂ + d₁f)}(d₂/2, d₁/2)` as an expression.
fn f_sf(ctx: &Context, d1: usize, d2: usize, f: &Q) -> Ex {
    if !f.is_positive() {
        return ctx.one();
    }
    let z = qu(d2) / (qu(d2) + qu(d1) * f);
    ex(ctx, &z).betainc_regularized(
        &ctx.rational(d2 as i64, 2),
        &ctx.rational(d1 as i64, 2),
        &ctx.zero(),
    )
}

fn result(
    ctx: &Context,
    statistic: Ex,
    p_value: Q,
    df: Option<Ex>,
    alt: Alternative,
) -> TestResult {
    TestResult {
        statistic,
        p_value: ex(ctx, &p_value),
        df,
        alternative: alt,
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// 1. Exact discrete tests
// ═══════════════════════════════════════════════════════════════════════════

/// `P(X = i)` for `i = 0..=n`, `X ~ Binomial(n, p)`, exact.
fn binomial_pmf_table(n: usize, p: &Q) -> Vec<Q> {
    let q = Q::one() - p;
    let mut pp = vec![Q::one(); n + 1];
    let mut qq = vec![Q::one(); n + 1];
    for i in 1..=n {
        pp[i] = &pp[i - 1] * p;
        qq[i] = &qq[i - 1] * &q;
    }
    let mut coef = Q::one();
    (0..=n)
        .map(|i| {
            if i > 0 {
                coef = &coef * qu(n + 1 - i) / qu(i);
            }
            &coef * &pp[i] * &qq[n - i]
        })
        .collect()
}

/// The sum of `pmf` over the indices selected by `keep`.
fn mass_where(pmf: &[Q], keep: impl Fn(usize, &Q) -> bool) -> Q {
    pmf.iter()
        .enumerate()
        .filter(|(i, m)| keep(*i, m))
        .fold(Q::zero(), |acc, (_, m)| acc + m)
}

/// scipy's two-sided p-value of a discrete exact test: the total mass of
/// the outcomes no more likely than the observed one, `Σ_{i : P(i) ≤ P(k)}
/// P(i)`.  (scipy compares with a relative slack of `1e-7`; here the
/// comparison is exact.)
fn two_sided_mass(pmf: &[Q], observed: usize) -> Q {
    let at = &pmf[observed];
    mass_where(pmf, |_, m| m <= at).min(Q::one())
}

/// Exact binomial test of `P(success) = p₀` from `k` successes in `n`
/// trials: the statistic is the proportion `k/n`; the p-value is
/// `P(X ≤ k)` (`Less`), `P(X ≥ k)` (`Greater`) or, two-sided, the total
/// probability of the outcomes no more likely than `k` (`Σ_{i: P(i) ≤ P(k)}
/// P(i)`), an exact rational.  `scipy.stats.binomtest(k, n, p, alternative)`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::linprog::q;
/// use symplex::stats::hypothesis::{binomial_test, Alternative};
///
/// let ctx = Context::new();
/// // scipy: binomtest(7, 10, 0.5).pvalue = 0.34375 = 11/32
/// let r = binomial_test(&ctx, 7, 10, &q(1, 2), Alternative::TwoSided)?;
/// assert_eq!(r.p_value_exact(), Some(q(11, 32)));
/// assert_eq!(r.statistic_exact(), Some(q(7, 10)));
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for `n = 0`, `k > n` or `p₀ ∉ [0, 1]`.
pub fn binomial_test(
    ctx: &Context,
    k: usize,
    n: usize,
    p0: &Q,
    alt: Alternative,
) -> Result<TestResult, SymplexError> {
    const OP: &str = "binomial_test";
    if n == 0 {
        return Err(invalid(OP, "the number of trials must be positive"));
    }
    if k > n {
        return Err(invalid(OP, format!("{k} successes exceed {n} trials")));
    }
    if p0.is_negative() || *p0 > Q::one() {
        return Err(invalid(OP, "the null proportion must lie in [0, 1]"));
    }
    let pmf = binomial_pmf_table(n, p0);
    let p = match alt {
        Alternative::Less => mass_where(&pmf, |i, _| i <= k),
        Alternative::Greater => mass_where(&pmf, |i, _| i >= k),
        Alternative::TwoSided => two_sided_mass(&pmf, k),
    };
    Ok(result(ctx, ex(ctx, &(qu(k) / qu(n))), p, None, alt))
}

/// `P(X = x)` of `Hypergeometric(N = n₁ + n₂, n₁, n)` for every `x` in the
/// support `max(0, n − n₂) ..= min(n, n₁)`; returns `(lowest x, masses)`.
fn hypergeometric_pmf_table(n1: usize, n2: usize, n: usize) -> (usize, Vec<Q>) {
    let lo = n.saturating_sub(n2);
    let hi = n.min(n1);
    let total = data::binomial_q(n1 + n2, n);
    let pmf = (lo..=hi)
        .map(|x| data::binomial_q(n1, x) * data::binomial_q(n2, n - x) / &total)
        .collect();
    (lo, pmf)
}

/// Fisher's exact test on a 2×2 table `[[a, b], [c, d]]`: the statistic is
/// the sample odds ratio `ad/(bc)` (`+∞` when `bc = 0`); the p-value is
/// exact from the hypergeometric distribution of `a` given the margins —
/// `P(X ≤ a)` (`Less`), `P(X ≥ a)` (`Greater`), or two-sided the total mass
/// of the tables no more likely than the observed one.
/// `scipy.stats.fisher_exact(table, alternative)`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::linprog::{q, qi};
/// use symplex::stats::hypothesis::{fisher_exact, Alternative};
///
/// let ctx = Context::new();
/// // scipy: fisher_exact([[8, 2], [1, 5]]) → statistic 20.0, pvalue 0.034965034965034975 = 5/143
/// let r = fisher_exact(&ctx, [[8, 2], [1, 5]], Alternative::TwoSided)?;
/// assert_eq!(r.statistic_exact(), Some(qi(20)));
/// assert_eq!(r.p_value_exact(), Some(q(5, 143)));
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] if a row or a column of the table is
/// empty (the test is undefined; scipy reports `p = 1`, odds ratio NaN).
pub fn fisher_exact(
    ctx: &Context,
    table: [[usize; 2]; 2],
    alt: Alternative,
) -> Result<TestResult, SymplexError> {
    const OP: &str = "fisher_exact";
    let [[a, b], [c, d]] = table;
    let (n1, n2, n) = (a + b, c + d, a + c);
    if n1 == 0 || n2 == 0 || n == 0 || b + d == 0 {
        return Err(invalid(OP, "a row or a column of the table is empty"));
    }
    let (lo, pmf) = hypergeometric_pmf_table(n1, n2, n);
    let idx = a - lo;
    let p = match alt {
        Alternative::Less => mass_where(&pmf, |i, _| i <= idx),
        Alternative::Greater => mass_where(&pmf, |i, _| i >= idx),
        Alternative::TwoSided => two_sided_mass(&pmf, idx),
    };
    let odds = if b * c == 0 {
        ctx.infinity()
    } else {
        ex(ctx, &(qu(a * d) / qu(b * c)))
    };
    Ok(result(ctx, odds, p, None, alt))
}

/// McNemar's test on the discordant counts `b` (row 1 / column 2) and `c`
/// (row 2 / column 1) of a paired 2×2 table.  `exact`: statistic
/// `min(b, c)`, p-value `min(1, 2·P(Binomial(b + c, ½) ≤ min(b, c)))`, an
/// exact rational.  Otherwise the χ² statistic `(|b − c| − 1)²/(b + c)`
/// (`correction`) or `(b − c)²/(b + c)` with `P(χ²₁ ≥ ·)`.
/// `statsmodels.stats.contingency_tables.mcnemar(table, exact, correction)`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::linprog::q;
/// use symplex::stats::hypothesis::mcnemar_test;
///
/// let ctx = Context::new();
/// // statsmodels: mcnemar([[100, 5], [15, 100]], exact=True) → statistic 5.0, pvalue 0.04138946533203125 = 5425/131072
/// let r = mcnemar_test(&ctx, 5, 15, true, true)?;
/// assert_eq!(r.p_value_exact(), Some(q(5425, 131072)));
/// // statsmodels: mcnemar(..., exact=False, correction=True) → statistic 4.05, pvalue 0.0441713449084427
/// let r = mcnemar_test(&ctx, 5, 15, false, true)?;
/// assert_eq!(r.statistic_exact(), Some(q(81, 20)));
/// assert!((r.p_value_f64()? - 0.044_171_344_908_442_7).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] when `b + c = 0` (no discordant pairs).
pub fn mcnemar_test(
    ctx: &Context,
    b: usize,
    c: usize,
    exact: bool,
    correction: bool,
) -> Result<TestResult, SymplexError> {
    const OP: &str = "mcnemar_test";
    let n = b + c;
    if n == 0 {
        return Err(invalid(OP, "there are no discordant pairs"));
    }
    if exact {
        let k = b.min(c);
        let pmf = binomial_pmf_table(n, &Q::new(BigInt::one(), BigInt::from(2)));
        let p = (mass_where(&pmf, |i, _| i <= k) * qi(2)).min(Q::one());
        return Ok(result(
            ctx,
            ctx.int(k as i64),
            p,
            None,
            Alternative::TwoSided,
        ));
    }
    let diff = qu(b.max(c) - b.min(c)) - if correction { Q::one() } else { Q::zero() };
    let stat = &diff * &diff / qu(n);
    Ok(TestResult {
        statistic: ex(ctx, &stat),
        p_value: chi_squared_sf_q(ctx, 1, &stat),
        df: Some(ctx.one()),
        alternative: Alternative::TwoSided,
    })
}

/// The sign test of `median = μ₀`: with `n₊` observations above and `n₋`
/// below `μ₀` (ties dropped), the statistic is `M = (n₊ − n₋)/2` and the
/// p-value is the exact binomial test of `n₊` out of `n₊ + n₋` with
/// `p = ½` (`statsmodels.stats.descriptivestats.sign_test(x, mu0)`, which
/// is two-sided; `Greater` is `P(X ≥ n₊)`, `Less` is `P(X ≤ n₊)`).
///
/// ```
/// use symplex::prelude::*;
/// use symplex::linprog::{q, qi};
/// use symplex::stats::data::from_i64;
/// use symplex::stats::hypothesis::{sign_test, Alternative};
///
/// let ctx = Context::new();
/// let x = from_i64(&[3, 5, 7, 8, 9, 11, 12, 15, 2, 6]);
/// // statsmodels: sign_test(x, mu0=5) = (2.5, 0.1796875)   (7 above, 2 below; 0.1796875 = 23/128)
/// let r = sign_test(&ctx, &x, &qi(5), Alternative::TwoSided)?;
/// assert_eq!(r.statistic_exact(), Some(q(5, 2)));
/// assert_eq!(r.p_value_exact(), Some(q(23, 128)));
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] when every observation equals `μ₀`.
pub fn sign_test(
    ctx: &Context,
    x: &[Q],
    mu0: &Q,
    alt: Alternative,
) -> Result<TestResult, SymplexError> {
    const OP: &str = "sign_test";
    let pos = x.iter().filter(|v| *v > mu0).count();
    let neg = x.iter().filter(|v| *v < mu0).count();
    if pos + neg == 0 {
        return Err(invalid(OP, "every observation equals the null median"));
    }
    let half = Q::new(BigInt::one(), BigInt::from(2));
    let inner = binomial_test(ctx, pos, pos + neg, &half, alt)?;
    let m = (qu(pos) - qu(neg)) / qi(2);
    Ok(TestResult {
        statistic: ex(ctx, &m),
        ..inner
    })
}

// ═══════════════════════════════════════════════════════════════════════════
// 2. Categorical
// ═══════════════════════════════════════════════════════════════════════════

/// A table of counts from integer rows (a convenience for the table-valued
/// functions of this module).
///
/// ```
/// use symplex::stats::hypothesis::counts;
/// use symplex::linprog::qi;
/// let t = counts(&[&[10, 20], &[30, 40]]);
/// assert_eq!(t[1][0], qi(30));
/// ```
#[must_use]
pub fn counts(rows: &[&[i64]]) -> Vec<Vec<Q>> {
    rows.iter().map(|r| data::from_i64(r)).collect()
}

/// Check a rectangular table of non-negative entries; returns `(rows, cols)`.
fn check_table(op: &'static str, table: &[Vec<Q>]) -> Result<(usize, usize), SymplexError> {
    let r = table.len();
    let c = table.first().map_or(0, Vec::len);
    if r == 0 || c == 0 {
        return Err(invalid(op, "the table is empty"));
    }
    if table.iter().any(|row| row.len() != c) {
        return Err(invalid(op, "the table is not rectangular"));
    }
    if table.iter().flatten().any(Signed::is_negative) {
        return Err(invalid(op, "counts must be non-negative"));
    }
    Ok((r, c))
}

/// Row sums, column sums and the grand total.
fn margins(table: &[Vec<Q>]) -> (Vec<Q>, Vec<Q>, Q) {
    let cols = table[0].len();
    let rows: Vec<Q> = table.iter().map(|r| data::sum(r)).collect();
    let col_sums: Vec<Q> = (0..cols)
        .map(|j| table.iter().fold(Q::zero(), |acc, r| acc + &r[j]))
        .collect();
    let total = data::sum(&rows);
    (rows, col_sums, total)
}

/// `Eᵢⱼ = rowᵢ · colⱼ / N`, erroring on an empty row or column.
fn expected_counts(op: &'static str, table: &[Vec<Q>]) -> Result<Vec<Vec<Q>>, SymplexError> {
    let (rows, cols, total) = margins(table);
    if rows.iter().any(Zero::is_zero) || cols.iter().any(Zero::is_zero) {
        return Err(invalid(
            op,
            "a row or a column of the table is empty, so an expected count is zero",
        ));
    }
    Ok(rows
        .iter()
        .map(|r| cols.iter().map(|c| r * c / &total).collect())
        .collect())
}

/// `Σ (|O − E| − shift)² / E` over a table and its expected counts.
fn pearson_statistic(observed: &[Vec<Q>], expected: &[Vec<Q>], shift: &Q) -> Q {
    observed
        .iter()
        .zip(expected)
        .flat_map(|(o, e)| o.iter().zip(e))
        .fold(Q::zero(), |acc, (o, e)| {
            let d = (o - e).abs() - shift;
            acc + &d * &d / e
        })
}

/// Pearson's χ² test of independence on an `r × c` table of counts:
/// `χ² = Σ (Oᵢⱼ − Eᵢⱼ)²/Eᵢⱼ` with `Eᵢⱼ = rowᵢ · colⱼ / N`, `df = (r−1)(c−1)`;
/// with `correction` and `df = 1` Yates' `(|O − E| − ½)²` (exactly as
/// scipy, without clamping).  The statistic and the expected counts are
/// exact rationals; the p-value is `P(χ²_df ≥ χ²)`.
/// `scipy.stats.chi2_contingency(table, correction)`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::linprog::q;
/// use symplex::stats::hypothesis::{chi_square_independence, counts};
///
/// let ctx = Context::new();
/// let t = counts(&[&[10, 20, 30], &[6, 9, 17]]);
/// // scipy: chi2_contingency(t) → statistic 0.27157465150403504, dof 2, pvalue 0.873028283380073
/// let r = chi_square_independence(&ctx, &t, true)?;
/// assert_eq!(r.df, 2);
/// assert_eq!(r.expected[0][0], q(240, 23));
/// assert!((r.p_value_f64()? - 0.873_028_283_380_073).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for a table with fewer than two rows
/// or columns, a ragged or negative table, or an empty row / column.
pub fn chi_square_independence(
    ctx: &Context,
    table: &[Vec<Q>],
    correction: bool,
) -> Result<ChiSquareResult, SymplexError> {
    const OP: &str = "chi_square_independence";
    let (r, c) = check_table(OP, table)?;
    if r < 2 || c < 2 {
        return Err(invalid(
            OP,
            "the table needs at least two rows and two columns",
        ));
    }
    let expected = expected_counts(OP, table)?;
    let df = (r - 1) * (c - 1);
    let shift = if correction && df == 1 {
        Q::new(BigInt::one(), BigInt::from(2))
    } else {
        Q::zero()
    };
    let statistic = pearson_statistic(table, &expected, &shift);
    Ok(ChiSquareResult {
        p_value: chi_squared_sf_q(ctx, df, &statistic),
        statistic,
        df,
        expected,
    })
}

/// Pearson's χ² goodness-of-fit test of observed counts against expected
/// ones (`None`: all equal to the mean count): `χ² = Σ (O − E)²/E`, `df = k −
/// 1 − ddof`.  The expected counts must sum to the observed total (exactly;
/// scipy tolerates `1e-8`).  `scipy.stats.chisquare(f_obs, f_exp, ddof)`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::linprog::qi;
/// use symplex::stats::data::from_i64;
/// use symplex::stats::hypothesis::chi_square_goodness_of_fit;
///
/// let ctx = Context::new();
/// // scipy: chisquare([16, 18, 16, 14, 12, 12]) → statistic 2.0, pvalue 0.8491450360846096
/// let r = chi_square_goodness_of_fit(&ctx, &from_i64(&[16, 18, 16, 14, 12, 12]), None, 0)?;
/// assert_eq!(r.statistic, qi(2));
/// assert_eq!(r.df, 5);
/// assert!((r.p_value_f64()? - 0.849_145_036_084_609_6).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for fewer than two categories, a
/// non-positive expected count, mismatched lengths or totals, or `ddof ≥ k − 1`.
pub fn chi_square_goodness_of_fit(
    ctx: &Context,
    observed: &[Q],
    expected: Option<&[Q]>,
    ddof: usize,
) -> Result<ChiSquareResult, SymplexError> {
    const OP: &str = "chi_square_goodness_of_fit";
    let k = observed.len();
    if k < 2 {
        return Err(invalid(OP, "at least two categories are needed"));
    }
    if observed.iter().any(Signed::is_negative) {
        return Err(invalid(OP, "observed counts must be non-negative"));
    }
    let expected: Vec<Q> = match expected {
        Some(e) => {
            if e.len() != k {
                return Err(invalid(OP, "observed and expected have different lengths"));
            }
            if e.iter().any(|v| !v.is_positive()) {
                return Err(invalid(OP, "expected counts must be positive"));
            }
            if data::sum(e) != data::sum(observed) {
                return Err(invalid(
                    OP,
                    "observed and expected counts have different totals",
                ));
            }
            e.to_vec()
        }
        None => {
            let m = data::mean(observed)?;
            if !m.is_positive() {
                return Err(invalid(OP, "the observed counts are all zero"));
            }
            vec![m; k]
        }
    };
    if ddof + 1 >= k {
        return Err(invalid(
            OP,
            format!("ddof = {ddof} leaves no degrees of freedom"),
        ));
    }
    let df = k - 1 - ddof;
    let statistic = observed
        .iter()
        .zip(&expected)
        .fold(Q::zero(), |acc, (o, e)| {
            let d = o - e;
            acc + &d * &d / e
        });
    Ok(ChiSquareResult {
        p_value: chi_squared_sf_q(ctx, df, &statistic),
        statistic,
        df,
        expected: vec![expected],
    })
}

/// The likelihood-ratio (G) test of independence on an `r × c` table:
/// `G = 2 Σ Oᵢⱼ ln(Oᵢⱼ/Eᵢⱼ)` (cells with `O = 0` contribute `0`), an exact
/// expression, with `P(χ²_{(r−1)(c−1)} ≥ G)`.  No continuity correction:
/// `scipy.stats.chi2_contingency(table, correction=False, lambda_='log-likelihood')`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::stats::hypothesis::{counts, g_test};
///
/// let ctx = Context::new();
/// let t = counts(&[&[10, 20, 30], &[6, 9, 17]]);
/// // scipy: chi2_contingency(t, correction=False, lambda_='log-likelihood')
/// //        → statistic 0.2740265420246606, pvalue 0.8719586542812721
/// let r = g_test(&ctx, &t)?;
/// assert!((r.statistic_f64()? - 0.274_026_542_024_660_6).abs() < 1e-12);
/// assert!((r.p_value_f64()? - 0.871_958_654_281_272_1).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// As [`chi_square_independence`].
pub fn g_test(ctx: &Context, table: &[Vec<Q>]) -> Result<TestResult, SymplexError> {
    const OP: &str = "g_test";
    let (r, c) = check_table(OP, table)?;
    if r < 2 || c < 2 {
        return Err(invalid(
            OP,
            "the table needs at least two rows and two columns",
        ));
    }
    let expected = expected_counts(OP, table)?;
    let df = (r - 1) * (c - 1);
    let mut terms: Vec<Ex> = Vec::new();
    for (o_row, e_row) in table.iter().zip(&expected) {
        for (o, e) in o_row.iter().zip(e_row) {
            if o.is_zero() || o == e {
                continue;
            }
            terms.push(ex(ctx, o) * (ex(ctx, o) / ex(ctx, e)).ln());
        }
    }
    if terms.is_empty() {
        return Ok(TestResult {
            statistic: ctx.zero(),
            p_value: ctx.one(),
            df: Some(ctx.int(df as i64)),
            alternative: Alternative::TwoSided,
        });
    }
    let sum = terms.into_iter().fold(ctx.zero(), |acc, t| acc + t);
    let statistic = ctx.int(2) * sum;
    Ok(TestResult {
        p_value: chi_squared_sf(ctx, df, &statistic),
        statistic,
        df: Some(ctx.int(df as i64)),
        alternative: Alternative::TwoSided,
    })
}

/// Cramér's `V = √(χ² / (N · min(r − 1, c − 1)))` of an `r × c` table
/// (without Yates' correction), as an exact expression.
/// `scipy.stats.contingency.association(table, method='cramer', correction=False)`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::stats::hypothesis::{counts, cramers_v};
///
/// let ctx = Context::new();
/// // scipy: association([[10, 20, 30], [6, 9, 17]], method='cramer', correction=False) = 0.05433137570422292
/// let v = cramers_v(&ctx, &counts(&[&[10, 20, 30], &[6, 9, 17]]))?;
/// assert!((v.eval_f64()? - 0.054_331_375_704_222_92).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// As [`chi_square_independence`].
pub fn cramers_v(ctx: &Context, table: &[Vec<Q>]) -> Result<Ex, SymplexError> {
    let r = chi_square_independence(ctx, table, false)?;
    let (rows, cols) = (table.len(), table[0].len());
    let (_, _, total) = margins(table);
    let k = qu((rows - 1).min(cols - 1));
    Ok(ex(ctx, &(r.statistic / (total * k))).sqrt().simplify())
}

fn check_2x2_margins(op: &'static str, table: [[usize; 2]; 2]) -> Result<(), SymplexError> {
    let [[a, b], [c, d]] = table;
    if a + b == 0 || c + d == 0 || a + c == 0 || b + d == 0 {
        return Err(invalid(op, "a row or a column of the table is empty"));
    }
    Ok(())
}

/// The φ coefficient of a 2×2 table `[[a, b], [c, d]]`:
/// `(ad − bc) / √((a+b)(c+d)(a+c)(b+d))`, as an exact expression (Pearson's
/// `r` of the two indicator variables; `φ² = χ²/N` without correction).
///
/// ```
/// use symplex::prelude::*;
/// use symplex::stats::hypothesis::phi_coefficient;
///
/// let ctx = Context::new();
/// // numpy: corrcoef of the indicator vectors of [[20, 10], [5, 15]] = 0.40824829046386296
/// let phi = phi_coefficient(&ctx, [[20, 10], [5, 15]])?;
/// assert!((phi.eval_f64()? - 0.408_248_290_463_862_96).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for an empty row or column.
pub fn phi_coefficient(ctx: &Context, table: [[usize; 2]; 2]) -> Result<Ex, SymplexError> {
    const OP: &str = "phi_coefficient";
    check_2x2_margins(OP, table)?;
    let [[a, b], [c, d]] = table;
    let num = qu(a * d) - qu(b * c);
    let den = qu((a + b) * (c + d) * (a + c) * (b + d));
    Ok((ex(ctx, &num) / ex(ctx, &den).sqrt()).simplify())
}

/// Wald interval `exp(ln θ̂ ± z_{1−α/2} · se)`.
fn log_wald_ci(estimate: &Q, se: f64, confidence: f64) -> (f64, f64) {
    let z = norm_isf((1.0 - confidence) / 2.0);
    let log = q_to_f64(estimate).ln();
    ((log - z * se).exp(), (log + z * se).exp())
}

/// The sample odds ratio `ad/(bc)` of a 2×2 table `[[a, b], [c, d]]` with
/// the log-scale Wald interval `exp(ln OR ± z √(1/a + 1/b + 1/c + 1/d))`.
/// `statsmodels.stats.contingency_tables.Table2x2(t).oddsratio` /
/// `.oddsratio_confint(alpha)`.
///
/// ```
/// use symplex::linprog::q;
/// use symplex::stats::hypothesis::odds_ratio;
///
/// // statsmodels: Table2x2([[20, 10], [5, 15]]).oddsratio = 6.0,
/// //              .oddsratio_confint(0.05) = (1.6931795592741443, 21.261773332199304)
/// let r = odds_ratio([[20, 10], [5, 15]], 0.95)?;
/// assert_eq!(r.estimate, q(6, 1));
/// assert!((r.ci.0 - 1.693_179_559_274_144_3).abs() < 1e-9);
/// assert!((r.ci.1 - 21.261_773_332_199_304).abs() < 1e-9);
/// # Ok::<(), symplex::prelude::SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] if any cell is zero (the estimate or
/// its standard error is infinite) or `confidence ∉ (0, 1)`.
pub fn odds_ratio(table: [[usize; 2]; 2], confidence: f64) -> Result<RatioEstimate, SymplexError> {
    const OP: &str = "odds_ratio";
    check_unit_open(OP, "confidence", confidence)?;
    let [[a, b], [c, d]] = table;
    if a == 0 || b == 0 || c == 0 || d == 0 {
        return Err(invalid(
            OP,
            "every cell must be positive for a finite odds ratio",
        ));
    }
    let estimate = qu(a * d) / qu(b * c);
    let se = (1.0 / a as f64 + 1.0 / b as f64 + 1.0 / c as f64 + 1.0 / d as f64).sqrt();
    Ok(RatioEstimate {
        ci: log_wald_ci(&estimate, se, confidence),
        estimate,
        confidence,
    })
}

/// The relative risk of a 2×2 table `[[exposed cases, exposed non-cases],
/// [control cases, control non-cases]]`: `(a/(a+b)) / (c/(c+d))` with the
/// log-scale Wald interval `exp(ln RR ± z √(1/a − 1/(a+b) + 1/c − 1/(c+d)))`.
/// `scipy.stats.contingency.relative_risk(a, a+b, c, c+d)` and its
/// `.confidence_interval(confidence)`.
///
/// ```
/// use symplex::linprog::q;
/// use symplex::stats::hypothesis::relative_risk;
///
/// // scipy: relative_risk(20, 30, 5, 20).relative_risk = 2.6666666666666665,
/// //        .confidence_interval(0.95) = (1.198028521436089, 5.935677643623166)
/// let r = relative_risk([[20, 10], [5, 15]], 0.95)?;
/// assert_eq!(r.estimate, q(8, 3));
/// assert!((r.ci.0 - 1.198_028_521_436_089).abs() < 1e-9);
/// assert!((r.ci.1 - 5.935_677_643_623_166).abs() < 1e-9);
/// # Ok::<(), symplex::prelude::SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] if either case count is zero or a
/// row is empty, or `confidence ∉ (0, 1)`.
pub fn relative_risk(
    table: [[usize; 2]; 2],
    confidence: f64,
) -> Result<RatioEstimate, SymplexError> {
    const OP: &str = "relative_risk";
    check_unit_open(OP, "confidence", confidence)?;
    let [[a, b], [c, d]] = table;
    let (n1, n2) = (a + b, c + d);
    if a == 0 || c == 0 {
        return Err(invalid(
            OP,
            "both case counts must be positive for a finite relative risk",
        ));
    }
    let estimate = (qu(a) / qu(n1)) / (qu(c) / qu(n2));
    let se = (1.0 / a as f64 - 1.0 / n1 as f64 + 1.0 / c as f64 - 1.0 / n2 as f64).sqrt();
    Ok(RatioEstimate {
        ci: log_wald_ci(&estimate, se, confidence),
        estimate,
        confidence,
    })
}

/// Cohen's `h = 2 asin √p₁ − 2 asin √p₂`, the effect size of a difference
/// of proportions, as an exact expression.
/// `statsmodels.stats.proportion.proportion_effectsize(p1, p2)`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::linprog::q;
/// use symplex::stats::hypothesis::cohens_h;
///
/// let ctx = Context::new();
/// // statsmodels: proportion_effectsize(0.5, 0.4) = 0.20135792079033088
/// let h = cohens_h(&ctx, &q(1, 2), &q(2, 5))?;
/// assert!((h.eval_f64()? - 0.201_357_920_790_330_88).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for a proportion outside `[0, 1]`.
pub fn cohens_h(ctx: &Context, p1: &Q, p2: &Q) -> Result<Ex, SymplexError> {
    const OP: &str = "cohens_h";
    for p in [p1, p2] {
        if p.is_negative() || *p > Q::one() {
            return Err(invalid(OP, "proportions must lie in [0, 1]"));
        }
    }
    let two = ctx.int(2);
    Ok(&two * ex(ctx, p1).sqrt().asin() - &two * ex(ctx, p2).sqrt().asin())
}

// ═══════════════════════════════════════════════════════════════════════════
// 3. Means and proportions
// ═══════════════════════════════════════════════════════════════════════════

fn student_result(ctx: &Context, df: &Q, stat: &RootRatio, alt: Alternative) -> TestResult {
    TestResult {
        statistic: stat.to_ex(ctx),
        p_value: student_p_value(ctx, df, stat, alt),
        df: Some(ex(ctx, df)),
        alternative: alt,
    }
}

fn normal_result(ctx: &Context, stat: &RootRatio, alt: Alternative) -> TestResult {
    TestResult {
        statistic: stat.to_ex(ctx),
        p_value: normal_p_value(ctx, stat, alt),
        df: None,
        alternative: alt,
    }
}

/// One-sample Student t-test of `mean = μ₀`: `t = (x̄ − μ₀) / (s/√n)` with
/// the sample standard deviation `s`, `df = n − 1`.  The statistic is exact
/// (`(x̄ − μ₀)√n / √s²`); the p-value is the exact Student-t tail
/// `I_{ν/(t²+ν)}(ν/2, ½)` (two-sided) or half of it / its complement.
/// `scipy.stats.ttest_1samp(x, popmean, alternative)`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::linprog::qi;
/// use symplex::stats::data::from_i64;
/// use symplex::stats::hypothesis::{t_test_one_sample, Alternative};
///
/// let ctx = Context::new();
/// let x = from_i64(&[5, 7, 8, 9, 10, 12]);
/// // scipy: ttest_1samp(x, 6, alternative='greater') → statistic 2.521097420448054, pvalue 0.026551745875898917
/// let r = t_test_one_sample(&ctx, &x, &qi(6), Alternative::Greater)?;
/// assert!((r.statistic_f64()? - 2.521_097_420_448_054).abs() < 1e-12);
/// assert!((r.p_value_f64()? - 0.026_551_745_875_898_917).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for fewer than two observations or a
/// constant sample (zero variance).
pub fn t_test_one_sample(
    ctx: &Context,
    x: &[Q],
    mu0: &Q,
    alt: Alternative,
) -> Result<TestResult, SymplexError> {
    const OP: &str = "t_test_one_sample";
    check_sample(OP, "the sample", x, 2)?;
    let n = x.len();
    let var = data::variance(x, Ddof::Sample)?;
    if var.is_zero() {
        return Err(invalid(OP, "the sample is constant (zero variance)"));
    }
    let stat = RootRatio {
        num: data::mean(x)? - mu0,
        var: var / qu(n),
    };
    Ok(student_result(ctx, &qu(n - 1), &stat, alt))
}

/// Two-sample t-test of `mean(x) = mean(y)`.  `equal_var`: Student's test
/// with the pooled variance `s²ₚ = ((n₁−1)s₁² + (n₂−1)s₂²)/(n₁+n₂−2)`,
/// `t = (x̄ − ȳ)/√(s²ₚ(1/n₁ + 1/n₂))`, `df = n₁ + n₂ − 2`.  Otherwise Welch's
/// test `t = (x̄ − ȳ)/√(s₁²/n₁ + s₂²/n₂)` with the Welch–Satterthwaite
/// `ν = (s₁²/n₁ + s₂²/n₂)² / ((s₁²/n₁)²/(n₁−1) + (s₂²/n₂)²/(n₂−1))`, an exact
/// rational.  `scipy.stats.ttest_ind(x, y, equal_var, alternative)`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::linprog::q;
/// use symplex::stats::data::from_i64;
/// use symplex::stats::hypothesis::{t_test_two_sample, Alternative};
///
/// let ctx = Context::new();
/// let x = from_i64(&[20, 22, 19, 20, 22, 20, 21]);
/// let y = from_i64(&[28, 32, 36, 24, 29, 32]);
/// // scipy: ttest_ind(x, y, equal_var=False) → statistic -5.529270507777645,
/// //        pvalue 0.0017938799124542811, df 5.650760744239527 (exactly 3527433605/624240481)
/// let r = t_test_two_sample(&ctx, &x, &y, false, Alternative::TwoSided)?;
/// assert!((r.statistic_f64()? - -5.529_270_507_777_645).abs() < 1e-12);
/// assert!((r.p_value_f64()? - 0.001_793_879_912_454_281_1).abs() < 1e-12);
/// assert_eq!(r.df, Some(ctx.from_ratio(q(3_527_433_605, 624_240_481))));
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] if either sample has fewer than two
/// observations or both are constant.
pub fn t_test_two_sample(
    ctx: &Context,
    x: &[Q],
    y: &[Q],
    equal_var: bool,
    alt: Alternative,
) -> Result<TestResult, SymplexError> {
    const OP: &str = "t_test_two_sample";
    check_sample(OP, "the first sample", x, 2)?;
    check_sample(OP, "the second sample", y, 2)?;
    let (n1, n2) = (x.len(), y.len());
    let (v1, v2) = (
        data::variance(x, Ddof::Sample)?,
        data::variance(y, Ddof::Sample)?,
    );
    let num = data::mean(x)? - data::mean(y)?;
    if v1.is_zero() && v2.is_zero() {
        return Err(invalid(OP, "both samples are constant (zero variance)"));
    }
    if equal_var {
        let df = qu(n1 + n2 - 2);
        let pooled = (qu(n1 - 1) * &v1 + qu(n2 - 1) * &v2) / &df;
        let var = pooled * (qu(n1).recip() + qu(n2).recip());
        Ok(student_result(ctx, &df, &RootRatio { num, var }, alt))
    } else {
        let (a, b) = (&v1 / qu(n1), &v2 / qu(n2));
        let var = &a + &b;
        let df = &var * &var / (&a * &a / qu(n1 - 1) + &b * &b / qu(n2 - 1));
        Ok(student_result(ctx, &df, &RootRatio { num, var }, alt))
    }
}

/// Paired t-test: the one-sample test of the differences `xᵢ − yᵢ` against
/// `0`.  `scipy.stats.ttest_rel(x, y, alternative)`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::stats::data::from_i64;
/// use symplex::stats::hypothesis::{t_test_paired, Alternative};
///
/// let ctx = Context::new();
/// let before = from_i64(&[200, 190, 210, 180, 195, 205]);
/// let after = from_i64(&[190, 185, 200, 182, 190, 195]);
/// // scipy: ttest_rel(before, after) → statistic 3.258473117707668, pvalue 0.022483670687634263
/// let r = t_test_paired(&ctx, &before, &after, Alternative::TwoSided)?;
/// assert!((r.statistic_f64()? - 3.258_473_117_707_668).abs() < 1e-12);
/// assert!((r.p_value_f64()? - 0.022_483_670_687_634_263).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for unequal lengths, fewer than two
/// pairs, or constant differences.
pub fn t_test_paired(
    ctx: &Context,
    x: &[Q],
    y: &[Q],
    alt: Alternative,
) -> Result<TestResult, SymplexError> {
    const OP: &str = "t_test_paired";
    check_same_len(OP, x, y)?;
    let d: Vec<Q> = x.iter().zip(y).map(|(a, b)| a - b).collect();
    check_sample(OP, "the paired sample", &d, 2)?;
    if data::variance(&d, Ddof::Sample)?.is_zero() {
        return Err(invalid(OP, "the differences are constant (zero variance)"));
    }
    t_test_one_sample(ctx, &d, &Q::zero(), alt)
}

/// One-proportion z-test of `p = p₀` from `k` successes in `n` trials:
/// `z = (k/n − p₀) / √(p₀(1 − p₀)/n)` (the null variance), with the normal
/// tail `½ erfc(z/√2)`.
/// `statsmodels.stats.proportion.proportions_ztest(k, n, value=p0, prop_var=p0)`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::linprog::q;
/// use symplex::stats::hypothesis::{z_test_proportion, Alternative};
///
/// let ctx = Context::new();
/// // statsmodels: proportions_ztest(60, 100, value=0.5, prop_var=0.5) = (2.0, 0.04550026389635844)
/// let r = z_test_proportion(&ctx, 60, 100, &q(1, 2), Alternative::TwoSided)?;
/// assert_eq!(r.statistic, ctx.int(2));
/// assert!((r.p_value_f64()? - 0.045_500_263_896_358_44).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for `n = 0`, `k > n` or `p₀ ∉ (0, 1)`.
pub fn z_test_proportion(
    ctx: &Context,
    k: usize,
    n: usize,
    p0: &Q,
    alt: Alternative,
) -> Result<TestResult, SymplexError> {
    const OP: &str = "z_test_proportion";
    if n == 0 {
        return Err(invalid(OP, "the number of trials must be positive"));
    }
    if k > n {
        return Err(invalid(OP, format!("{k} successes exceed {n} trials")));
    }
    if !p0.is_positive() || *p0 >= Q::one() {
        return Err(invalid(
            OP,
            "the null proportion must lie strictly between 0 and 1",
        ));
    }
    let stat = RootRatio {
        num: qu(k) / qu(n) - p0,
        var: p0 * (Q::one() - p0) / qu(n),
    };
    Ok(normal_result(ctx, &stat, alt))
}

/// Two-proportion z-test of `p₁ = p₂` with the pooled estimate
/// `p̂ = (k₁ + k₂)/(n₁ + n₂)`: `z = (k₁/n₁ − k₂/n₂) / √(p̂(1 − p̂)(1/n₁ + 1/n₂))`.
/// `statsmodels.stats.proportion.proportions_ztest([k1, k2], [n1, n2])`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::stats::hypothesis::{two_proportion_z_test, Alternative};
///
/// let ctx = Context::new();
/// // statsmodels: proportions_ztest([45, 30], [100, 100]) = (2.1908902300206647, 0.028459736916310555)
/// let r = two_proportion_z_test(&ctx, 45, 100, 30, 100, Alternative::TwoSided)?;
/// assert!((r.statistic_f64()? - 2.190_890_230_020_664_7).abs() < 1e-12);
/// assert!((r.p_value_f64()? - 0.028_459_736_916_310_555).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for an empty sample, `k > n`, or a
/// pooled proportion of `0` or `1` (zero variance).
pub fn two_proportion_z_test(
    ctx: &Context,
    k1: usize,
    n1: usize,
    k2: usize,
    n2: usize,
    alt: Alternative,
) -> Result<TestResult, SymplexError> {
    const OP: &str = "two_proportion_z_test";
    if n1 == 0 || n2 == 0 {
        return Err(invalid(OP, "both samples must be non-empty"));
    }
    if k1 > n1 || k2 > n2 {
        return Err(invalid(OP, "successes exceed trials"));
    }
    let pooled = qu(k1 + k2) / qu(n1 + n2);
    if pooled.is_zero() || pooled.is_one() {
        return Err(invalid(
            OP,
            "the pooled proportion is 0 or 1 (zero variance)",
        ));
    }
    let stat = RootRatio {
        num: qu(k1) / qu(n1) - qu(k2) / qu(n2),
        var: &pooled * (Q::one() - &pooled) * (qu(n1).recip() + qu(n2).recip()),
    };
    Ok(normal_result(ctx, &stat, alt))
}

/// `(SS_between, SS_within, N, k)` of a collection of groups.
fn sums_of_squares(
    op: &'static str,
    groups: &[Vec<Q>],
) -> Result<(Q, Q, usize, usize), 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((ss_between, ss_within, all.len(), k))
}

/// One-way analysis of variance of `k` independent groups:
/// `F = (SS_between/(k−1)) / (SS_within/(N−k))` with the sums of squares and
/// `F` exact, `η² = SS_between/SS_total`, and `P(F_{k−1, N−k} ≥ F)` as an
/// exact expression.  `scipy.stats.f_oneway(*groups)`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::linprog::q;
/// use symplex::stats::data::from_i64;
/// use symplex::stats::hypothesis::anova_one_way;
///
/// let ctx = Context::new();
/// let g = [from_i64(&[6, 8, 4, 5, 3, 4]), from_i64(&[8, 12, 9, 11, 6, 8]), from_i64(&[13, 9, 11, 8, 7, 12])];
/// // scipy: f_oneway(*g) → statistic 9.264705882352942 (= 315/34), pvalue 0.0023987773293929083
/// let r = anova_one_way(&ctx, &g)?;
/// assert_eq!(r.f, q(315, 34));
/// assert_eq!((r.df_between, r.df_within), (2, 15));
/// assert!((r.p_value_f64()? - 0.002_398_777_329_392_908_3).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for fewer than two groups, an empty
/// group, `N ≤ k`, or zero within-group variance.
pub fn anova_one_way(ctx: &Context, groups: &[Vec<Q>]) -> Result<AnovaResult, SymplexError> {
    const OP: &str = "anova_one_way";
    let (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,
    })
}

/// The `confidence` interval for the mean, `x̄ ± t_{(1+c)/2, n−1} · s/√n`,
/// with the Student-t quantile found by Brent's method on the exact CDF
/// expression (evaluated numerically).  `scipy.stats.t.interval(c, n−1,
/// loc=mean, scale=sem)`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::stats::data::from_i64;
/// use symplex::stats::hypothesis::confidence_interval_mean;
///
/// let ctx = Context::new();
/// let x = from_i64(&[5, 7, 8, 9, 10, 12]);
/// // scipy: t.interval(0.95, 5, loc=mean(x), scale=sem(x)) = (5.9509296876164886, 11.049070312383511)
/// let (lo, hi) = confidence_interval_mean(&ctx, &x, 0.95)?;
/// assert!((lo - 5.950_929_687_616_488_6).abs() < 1e-9);
/// assert!((hi - 11.049_070_312_383_511).abs() < 1e-9);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for fewer than two observations or
/// `confidence ∉ (0, 1)`; the quantile's error if it fails to converge.
pub fn confidence_interval_mean(
    ctx: &Context,
    x: &[Q],
    confidence: f64,
) -> Result<(f64, f64), SymplexError> {
    const OP: &str = "confidence_interval_mean";
    check_sample(OP, "the sample", x, 2)?;
    check_unit_open(OP, "confidence", confidence)?;
    let n = x.len();
    let mean = q_to_f64(&data::mean(x)?);
    let sem = q_to_f64(&(data::variance(x, Ddof::Sample)? / qu(n))).sqrt();
    let t = student_t_quantile_f64(OP, ctx, (n - 1) as f64, (1.0 + confidence) / 2.0)?;
    Ok((mean - t * sem, mean + t * sem))
}

// ═══════════════════════════════════════════════════════════════════════════
// 4. Rank tests
// ═══════════════════════════════════════════════════════════════════════════

/// `Σ (t³ − t)` over the tie groups of a sample.
fn tie_term(x: &[Q]) -> Q {
    data::tie_sizes(x)
        .into_iter()
        .fold(Q::zero(), |acc, t| acc + qu(t * t * t - t))
}

/// Frequencies of `U = u`, `u = 0..=mn`, over the `C(m+n, m)` equally likely
/// arrangements of `m` and `n` distinct values: the coefficients of the
/// Gaussian binomial `[m+n choose m]_q = Π_{i=1}^{m} (1 − q^{n+i})/(1 − q^i)`.
fn mann_whitney_frequencies(m: usize, n: usize) -> Vec<BigInt> {
    let size = m * n + 1;
    let mut c = vec![BigInt::zero(); size];
    c[0] = BigInt::one();
    for i in 1..=m {
        let shift = n + i;
        for t in (shift..size).rev() {
            let v = c[t - shift].clone();
            c[t] -= v;
        }
        for t in i..size {
            let v = c[t - i].clone();
            c[t] += v;
        }
    }
    c
}

/// Frequencies of `T⁺ = s`, `s = 0..=n(n+1)/2`, over the `2ⁿ` equally likely
/// sign patterns of the ranks `1..=n` (the number of subsets with sum `s`).
fn signed_rank_frequencies(n: usize) -> Vec<BigInt> {
    let size = n * (n + 1) / 2 + 1;
    let mut c = vec![BigInt::zero(); size];
    c[0] = BigInt::one();
    for k in 1..=n {
        for s in (k..size).rev() {
            let v = c[s - k].clone();
            c[s] += v;
        }
    }
    c
}

/// Mahonian numbers `M(n, k)`, `k = 0..=cmax`: permutations of `n` with `k`
/// inversions (Kendall, "Rank Correlation Methods", ch. 5).
fn inversion_counts(n: usize, cmax: usize) -> Vec<BigInt> {
    let mut c = vec![BigInt::zero(); cmax + 1];
    c[0] = BigInt::one();
    for j in 2..=n {
        let mut s = c.clone();
        for t in 1..=cmax {
            let v = s[t - 1].clone();
            s[t] += v;
        }
        for k in 0..=cmax {
            c[k] = if k >= j {
                &s[k] - &s[k - j]
            } else {
                s[k].clone()
            };
        }
    }
    c
}

/// `Σ_{k ≤ upto} freq[k] / total`.
fn cdf_of(freq: &[BigInt], total: &BigInt, upto: usize) -> Q {
    let upto = upto.min(freq.len().saturating_sub(1));
    Q::new(sum_big(&freq[..=upto]), total.clone())
}

/// `Σ_{k ≥ from} freq[k] / total`.
fn sf_of(freq: &[BigInt], total: &BigInt, from: usize) -> Q {
    if from >= freq.len() {
        return Q::zero();
    }
    Q::new(sum_big(&freq[from..]), total.clone())
}

fn rank_statistic_to_index(op: &'static str, v: &Q) -> Result<usize, SymplexError> {
    if !v.is_integer() || v.is_negative() {
        return Err(invalid(
            op,
            "the exact distribution needs an integer statistic",
        ));
    }
    v.to_integer()
        .to_usize()
        .ok_or_else(|| invalid(op, "the statistic is too large"))
}

/// Continuity correction of a discrete statistic compared with a normal:
/// move `num` half a unit towards zero in the direction of the
/// alternative (`sign(num)` two-sided, `+1` greater, `−1` less).
fn continuity_shift(num: &Q, alt: Alternative) -> Q {
    let half = Q::new(BigInt::one(), BigInt::from(2));
    match alt {
        Alternative::Greater => num - half,
        Alternative::Less => num + half,
        Alternative::TwoSided => match num.cmp(&Q::zero()) {
            Ordering::Greater => num - half,
            Ordering::Less => num + half,
            Ordering::Equal => num.clone(),
        },
    }
}

/// Mann–Whitney U test of two independent samples.  With the ranks of the
/// pooled sample, `U₁ = R₁ − n₁(n₁+1)/2` is the statistic (for `x`; `U₂ =
/// n₁n₂ − U₁`).  `Exact`: the p-value is `P(U ≥ U₁)` (`Greater`), `P(U ≥
/// U₂)` (`Less`) or `min(1, 2 P(U ≥ max(U₁, U₂)))` from the exact
/// distribution of `U` (no ties allowed), an exact rational.  `Asymptotic`:
/// `z = (U₁ − n₁n₂/2 ∓ ½) / σ` with `σ² = n₁n₂/12 · ((N+1) − Σ(t³−t)/(N(N−1)))`
/// (tie correction) and the normal tail.
/// `scipy.stats.mannwhitneyu(x, y, use_continuity, alternative, method)`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::linprog::{q, qi};
/// use symplex::stats::data::from_i64;
/// use symplex::stats::hypothesis::{mann_whitney_u, Alternative, RankMethod};
///
/// let ctx = Context::new();
/// let males = from_i64(&[19, 22, 16, 29, 24]);
/// let females = from_i64(&[20, 11, 17, 12]);
/// // scipy: mannwhitneyu(males, females, method='exact') = (17.0, 0.1111111111111111)
/// let r = mann_whitney_u(&ctx, &males, &females, Alternative::TwoSided, RankMethod::Exact)?;
/// assert_eq!(r.statistic_exact(), Some(qi(17)));
/// assert_eq!(r.p_value_exact(), Some(q(1, 9)));
/// // scipy: mannwhitneyu(males, females, method='asymptotic') → pvalue 0.11134688653314039
/// let r = mann_whitney_u(&ctx, &males, &females, Alternative::TwoSided, RankMethod::Asymptotic { continuity: true })?;
/// assert!((r.p_value_f64()? - 0.111_346_886_533_140_39).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for an empty sample, ties with
/// `Exact`, or a pooled sample that is constant (`Asymptotic`).
pub fn mann_whitney_u(
    ctx: &Context,
    x: &[Q],
    y: &[Q],
    alt: Alternative,
    method: RankMethod,
) -> Result<TestResult, SymplexError> {
    const OP: &str = "mann_whitney_u";
    check_sample(OP, "the first sample", x, 1)?;
    check_sample(OP, "the second sample", y, 1)?;
    let (n1, n2) = (x.len(), y.len());
    let pooled: Vec<Q> = x.iter().chain(y).cloned().collect();
    let ranks = data::ranks(&pooled);
    let r1 = data::sum(&ranks[..n1]);
    let u1 = r1 - qu(n1 * (n1 + 1) / 2);
    let u2 = qu(n1 * n2) - &u1;
    let statistic = ex(ctx, &u1);
    match method {
        RankMethod::Exact => {
            if !data::tie_sizes(&pooled).is_empty() {
                return Err(invalid(
                    OP,
                    "the exact method needs a pooled sample without ties",
                ));
            }
            let u = match alt {
                Alternative::Greater => &u1,
                Alternative::Less => &u2,
                Alternative::TwoSided => (&u1).max(&u2),
            };
            let k = rank_statistic_to_index(OP, u)?;
            let freq = mann_whitney_frequencies(n1, n2);
            let total = sum_big(&freq);
            let sf = sf_of(&freq, &total, k);
            let p = match alt {
                Alternative::TwoSided => (sf * qi(2)).min(Q::one()),
                _ => sf,
            };
            Ok(result(ctx, statistic, p, None, alt))
        }
        RankMethod::Asymptotic { continuity } => {
            let n = n1 + n2;
            let var = qu(n1 * n2) / qi(12) * (qu(n + 1) - tie_term(&pooled) / qu(n * (n - 1)));
            if !var.is_positive() {
                return Err(invalid(OP, "every observation is identical"));
            }
            let mut num = &u1 - qu(n1 * n2) / qi(2);
            if continuity {
                num = continuity_shift(&num, alt);
            }
            Ok(TestResult {
                statistic,
                p_value: normal_p_value(ctx, &RootRatio { num, var }, alt),
                df: None,
                alternative: alt,
            })
        }
    }
}

/// Wilcoxon signed-rank test of the differences `dᵢ = xᵢ − yᵢ` (or of `x`
/// alone) against a symmetric distribution about `0`.  Zero differences
/// are dropped (`zero_method='wilcox'`); `|d|` is ranked with average
/// ranks; `T⁺`, `T⁻` are the rank sums of the positive and negative
/// differences.  The statistic is `min(T⁺, T⁻)` two-sided and `T⁺`
/// one-sided.  `Exact` (no ties in `|d|`): `P(T⁺ ≤ T⁺)` (`Less`), `P(T⁺ ≥
/// T⁺)` (`Greater`), `min(1, 2 min(cdf, sf))` two-sided, from the `2ⁿ` sign
/// patterns, an exact rational.  `Asymptotic`: `z = (T⁺ − n(n+1)/4 ∓ ½)/σ`
/// with `σ² = (n(n+1)(2n+1) − Σ(t³−t)/2)/24` and the normal tail.
/// `scipy.stats.wilcoxon(x, y, correction, alternative, method)`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::linprog::{q, qi};
/// use symplex::stats::data::from_i64;
/// use symplex::stats::hypothesis::{wilcoxon_signed_rank, Alternative, RankMethod};
///
/// let ctx = Context::new();
/// let x = from_i64(&[125, 115, 130, 140, 140, 115, 140, 125, 140, 135]);
/// let y = from_i64(&[110, 122, 125, 120, 140, 124, 123, 137, 134, 145]);
/// // scipy: wilcoxon(x, y, method='exact') = (18.0, 0.65234375)   (one zero difference dropped; 0.65234375 = 167/256)
/// let r = wilcoxon_signed_rank(&ctx, &x, Some(&y), Alternative::TwoSided, RankMethod::Exact)?;
/// assert_eq!(r.statistic_exact(), Some(qi(18)));
/// assert_eq!(r.p_value_exact(), Some(q(167, 256)));
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for unequal lengths, no non-zero
/// differences, or ties in `|d|` with `Exact`.
pub fn wilcoxon_signed_rank(
    ctx: &Context,
    x: &[Q],
    y: Option<&[Q]>,
    alt: Alternative,
    method: RankMethod,
) -> Result<TestResult, SymplexError> {
    const OP: &str = "wilcoxon_signed_rank";
    let d: Vec<Q> = match y {
        Some(y) => {
            check_same_len(OP, x, y)?;
            x.iter().zip(y).map(|(a, b)| a - b).collect()
        }
        None => x.to_vec(),
    };
    let d: Vec<Q> = d.into_iter().filter(|v| !v.is_zero()).collect();
    let n = d.len();
    if n == 0 {
        return Err(invalid(OP, "every difference is zero"));
    }
    let abs: Vec<Q> = d.iter().map(Signed::abs).collect();
    let ranks = data::ranks(&abs);
    let (mut r_plus, mut r_minus) = (Q::zero(), Q::zero());
    for (v, r) in d.iter().zip(&ranks) {
        if v.is_positive() {
            r_plus += r;
        } else {
            r_minus += r;
        }
    }
    let statistic = match alt {
        Alternative::TwoSided => ex(ctx, (&r_plus).min(&r_minus)),
        _ => ex(ctx, &r_plus),
    };
    match method {
        RankMethod::Exact => {
            if !data::tie_sizes(&abs).is_empty() {
                return Err(invalid(
                    OP,
                    "the exact method needs |differences| without ties",
                ));
            }
            let k = rank_statistic_to_index(OP, &r_plus)?;
            let freq = signed_rank_frequencies(n);
            let total = BigInt::one() << n;
            let (cdf, sf) = (cdf_of(&freq, &total, k), sf_of(&freq, &total, k));
            let p = match alt {
                Alternative::Less => cdf,
                Alternative::Greater => sf,
                Alternative::TwoSided => (cdf.min(sf) * qi(2)).min(Q::one()),
            };
            Ok(result(ctx, statistic, p, None, alt))
        }
        RankMethod::Asymptotic { continuity } => {
            let var = (qu(n * (n + 1) * (2 * n + 1)) - tie_term(&abs) / qi(2)) / qi(24);
            if !var.is_positive() {
                return Err(invalid(OP, "the variance of the rank sum is zero"));
            }
            let mut num = &r_plus - qu(n * (n + 1)) / qi(4);
            if continuity {
                num = continuity_shift(&num, alt);
            }
            Ok(TestResult {
                statistic,
                p_value: normal_p_value(ctx, &RootRatio { num, var }, alt),
                df: None,
                alternative: alt,
            })
        }
    }
}

/// Kruskal–Wallis H test of `k` independent groups: with the pooled ranks,
/// `H = 12/(N(N+1)) Σ Rᵢ²/nᵢ − 3(N+1)`, divided by the tie correction
/// `1 − Σ(t³−t)/(N³−N)`; `H` is exact, `df = k − 1`, and the p-value is
/// `P(χ²_{k−1} ≥ H)`.  `scipy.stats.kruskal(*groups)`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::linprog::q;
/// use symplex::stats::data::from_i64;
/// use symplex::stats::hypothesis::kruskal_wallis;
///
/// let ctx = Context::new();
/// let g = [from_i64(&[1, 3, 5, 7, 9]), from_i64(&[2, 4, 6, 8, 10]), from_i64(&[11, 12, 13, 14, 15])];
/// // scipy: kruskal(*g) → statistic 9.5, pvalue 0.008651695203120634
/// let r = kruskal_wallis(&ctx, &g)?;
/// assert_eq!(r.statistic_exact(), Some(q(19, 2)));
/// assert!((r.p_value_f64()? - 0.008_651_695_203_120_634).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for fewer than two groups, an empty
/// group, or identical observations throughout.
pub fn kruskal_wallis(ctx: &Context, groups: &[Vec<Q>]) -> Result<TestResult, SymplexError> {
    const OP: &str = "kruskal_wallis";
    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 pooled: Vec<Q> = groups.iter().flatten().cloned().collect();
    let n = pooled.len();
    let ranks = data::ranks(&pooled);
    let mut ssbn = Q::zero();
    let mut start = 0;
    for g in groups {
        let r = data::sum(&ranks[start..start + g.len()]);
        ssbn += &r * &r / qu(g.len());
        start += g.len();
    }
    let h = qi(12) / qu(n * (n + 1)) * ssbn - qu(3 * (n + 1));
    let correction = Q::one() - tie_term(&pooled) / qu(n * n * n - n);
    if correction.is_zero() {
        return Err(invalid(OP, "every observation is identical"));
    }
    let h = h / correction;
    Ok(TestResult {
        p_value: chi_squared_sf_q(ctx, k - 1, &h),
        statistic: ex(ctx, &h),
        df: Some(ctx.int(usize_to_i64(OP, k - 1)?)),
        alternative: Alternative::TwoSided,
    })
}

/// Friedman's test of `k ≥ 3` treatments measured on `n` blocks (`blocks[i]`
/// holds the `k` measurements of block `i`).  Ranking within each block,
/// `χ²_F = (12/(nk(k+1)) Σ Rⱼ² − 3n(k+1)) / (1 − Σ_blocks Σ(t³−t)/(n(k³−k)))`
/// is exact, `df = k − 1`, and the p-value is `P(χ²_{k−1} ≥ χ²_F)`.
/// `scipy.stats.friedmanchisquare(*columns)` (which takes the treatments as
/// separate arrays — the transpose of `blocks`).
///
/// ```
/// use symplex::prelude::*;
/// use symplex::linprog::q;
/// use symplex::stats::data::from_i64;
/// use symplex::stats::hypothesis::friedman;
///
/// let ctx = Context::new();
/// let blocks = [from_i64(&[10, 9, 8]), from_i64(&[9, 8, 6]), from_i64(&[7, 5, 4]), from_i64(&[8, 6, 3])];
/// // scipy: friedmanchisquare([10, 9, 7, 8], [9, 8, 5, 6], [8, 6, 4, 3]) → statistic 8.0, pvalue 0.018315638888734182
/// let r = friedman(&ctx, &blocks)?;
/// assert_eq!(r.statistic_exact(), Some(q(8, 1)));
/// assert!((r.p_value_f64()? - 0.018_315_638_888_734_182).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for fewer than three treatments, no
/// blocks, ragged blocks, or ties within every block throughout.
pub fn friedman(ctx: &Context, blocks: &[Vec<Q>]) -> Result<TestResult, SymplexError> {
    const OP: &str = "friedman";
    let n = blocks.len();
    if n == 0 {
        return Err(invalid(OP, "at least one block is needed"));
    }
    let k = blocks[0].len();
    if k < 3 {
        return Err(invalid(OP, "at least three treatments are needed"));
    }
    if blocks.iter().any(|b| b.len() != k) {
        return Err(invalid(
            OP,
            "every block must hold the same number of treatments",
        ));
    }
    let mut column_sums = vec![Q::zero(); k];
    let mut ties = Q::zero();
    for b in blocks {
        let r = data::ranks(b);
        for (s, v) in column_sums.iter_mut().zip(&r) {
            *s += v;
        }
        ties += tie_term(b);
    }
    let correction = Q::one() - ties / qu(k * (k * k - 1) * n);
    if correction.is_zero() {
        return Err(invalid(OP, "every block is constant"));
    }
    let ssbn = column_sums.iter().fold(Q::zero(), |acc, r| acc + r * r);
    let stat = (qi(12) / qu(k * n * (k + 1)) * ssbn - qu(3 * n * (k + 1))) / correction;
    Ok(TestResult {
        p_value: chi_squared_sf_q(ctx, k - 1, &stat),
        statistic: ex(ctx, &stat),
        df: Some(ctx.int(usize_to_i64(OP, k - 1)?)),
        alternative: Alternative::TwoSided,
    })
}

/// Spearman's rank correlation test: the statistic is `ρ`
/// ([`data::spearman`], exact); the p-value uses `t = ρ √((n−2)/(1−ρ²))`
/// with `n − 2` degrees of freedom (exact Student-t tail; `t²` is rational
/// even when `ρ` is not).  `|ρ| = 1` gives `p = 0` in the alternative's
/// direction.  `scipy.stats.spearmanr(x, y, alternative)`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::linprog::q;
/// use symplex::stats::data::from_i64;
/// use symplex::stats::hypothesis::{spearman_test, Alternative};
///
/// let ctx = Context::new();
/// let x = from_i64(&[1, 2, 3, 4, 5, 6, 7, 8]);
/// let y = from_i64(&[2, 1, 4, 3, 7, 8, 5, 6]);
/// // scipy: spearmanr(x, y) → statistic 0.7619047619047621 (= 16/21), pvalue 0.028004939153071815
/// let r = spearman_test(&ctx, &x, &y, Alternative::TwoSided)?;
/// assert_eq!(r.statistic, ctx.from_ratio(q(16, 21)));
/// assert!((r.p_value_f64()? - 0.028_004_939_153_071_815).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for unequal lengths, fewer than three
/// pairs, or a constant sample.
pub fn spearman_test(
    ctx: &Context,
    x: &[Q],
    y: &[Q],
    alt: Alternative,
) -> Result<TestResult, SymplexError> {
    const OP: &str = "spearman_test";
    check_same_len(OP, x, y)?;
    check_sample(OP, "the paired sample", x, 3)?;
    let n = x.len();
    let (rx, ry) = (data::ranks(x), data::ranks(y));
    let sxy = data::covariance(&rx, &ry, Ddof::Population)?;
    let sxx = data::variance(&rx, Ddof::Population)?;
    let syy = data::variance(&ry, Ddof::Population)?;
    if sxx.is_zero() || syy.is_zero() {
        return Err(invalid(OP, "a constant sample has no rank correlation"));
    }
    let rho = data::spearman(ctx, x, y)?;
    let df = qu(n - 2);
    let r2 = &sxy * &sxy / (&sxx * &syy);
    let one_minus = Q::one() - &r2;
    let p_value = if one_minus.is_zero() {
        // |ρ| = 1: t is infinite in the direction of sign(ρ).
        let extreme = match alt {
            Alternative::TwoSided => true,
            Alternative::Greater => sxy.is_positive(),
            Alternative::Less => sxy.is_negative(),
        };
        if extreme { ctx.zero() } else { ctx.one() }
    } else {
        // t = sxy / √(sxx·syy·(1 − ρ²)/(n − 2)) has t² = ρ²(n−2)/(1−ρ²).
        let stat = RootRatio {
            num: sxy,
            var: &sxx * &syy * &one_minus / &df,
        };
        student_p_value(ctx, &df, &stat, alt)
    };
    Ok(TestResult {
        statistic: rho,
        p_value,
        df: Some(ex(ctx, &df)),
        alternative: alt,
    })
}

/// Kendall's τ-b test.  The statistic is `τ_b = (C − D)/√((n₀−n₁)(n₀−n₂))`
/// ([`data::kendall_tau`], exact).  `exact` (no ties): the p-value from the
/// exact distribution of the number of inversions (Kendall, "Rank
/// Correlation Methods"), an exact rational — scipy's `method='exact'`.
/// Otherwise `z = (C − D)/√var` with the tie-corrected variance
/// `var = (m(2n+5) − Σt(t−1)(2t+5) − Σu(u−1)(2u+5))/18 + 2n₁n₂/m + x₀y₀/(9m(n−2))`,
/// `m = n(n−1)`, and the normal tail — scipy's `method='asymptotic'`.
/// `scipy.stats.kendalltau(x, y, method, alternative)`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::linprog::q;
/// use symplex::stats::data::from_i64;
/// use symplex::stats::hypothesis::{kendall_test, Alternative};
///
/// let ctx = Context::new();
/// let x = from_i64(&[1, 2, 3, 4, 5, 6, 7, 8]);
/// let y = from_i64(&[2, 1, 4, 3, 7, 8, 5, 6]);
/// // scipy: kendalltau(x, y, method='exact') → statistic 0.5714285714285714 (= 4/7), pvalue 0.06101190476190476 (= 41/672)
/// let r = kendall_test(&ctx, &x, &y, Alternative::TwoSided, true)?;
/// assert_eq!(r.statistic, ctx.from_ratio(q(4, 7)));
/// assert_eq!(r.p_value_exact(), Some(q(41, 672)));
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for unequal lengths, fewer than two
/// pairs, a constant sample, or ties with `exact`.
pub fn kendall_test(
    ctx: &Context,
    x: &[Q],
    y: &[Q],
    alt: Alternative,
    exact: bool,
) -> Result<TestResult, SymplexError> {
    const OP: &str = "kendall_test";
    check_same_len(OP, x, y)?;
    check_sample(OP, "the paired sample", x, 2)?;
    let n = x.len();
    let (mut con, mut dis) = (0usize, 0usize);
    for i in 0..n {
        for j in i + 1..n {
            let (sx, sy) = (x[i].cmp(&x[j]), y[i].cmp(&y[j]));
            if sx == Ordering::Equal || sy == Ordering::Equal {
                continue;
            }
            if sx == sy {
                con += 1;
            } else {
                dis += 1;
            }
        }
    }
    let tot = n * (n - 1) / 2;
    let tie_stats = |t: &[usize]| -> (usize, Q, Q) {
        // (pairs tied, Σ t(t−1)(t−2), Σ t(t−1)(2t+5))
        t.iter().fold((0, Q::zero(), Q::zero()), |(p, a, b), &t| {
            (
                p + t * (t - 1) / 2,
                a + qu(t * (t - 1) * (t - 2)),
                b + qu(t * (t - 1) * (2 * t + 5)),
            )
        })
    };
    let (xtie, x0, x1) = tie_stats(&data::tie_sizes(x));
    let (ytie, y0, y1) = tie_stats(&data::tie_sizes(y));
    if xtie == tot || ytie == tot {
        return Err(invalid(OP, "a constant sample has no rank correlation"));
    }
    let tau = data::kendall_tau(ctx, x, y)?;
    if !exact {
        let m = qu(n * (n - 1));
        let mut var = (&m * qu(2 * n + 5) - &x1 - &y1) / qi(18) + qi(2) * qu(xtie * ytie) / &m;
        if !x0.is_zero() && !y0.is_zero() {
            var += &x0 * &y0 / (qi(9) * &m * qu(n - 2));
        }
        if !var.is_positive() {
            return Err(invalid(OP, "the variance of the statistic is zero"));
        }
        let stat = RootRatio {
            num: qu(con) - qu(dis),
            var,
        };
        return Ok(TestResult {
            statistic: tau,
            p_value: normal_p_value(ctx, &stat, alt),
            df: None,
            alternative: alt,
        });
    }
    if xtie > 0 || ytie > 0 {
        return Err(invalid(OP, "the exact method needs samples without ties"));
    }
    // scipy: c = concordant count; work in the left tail of the symmetric
    // distribution of the inversion count.
    let c = tot - dis;
    let in_right_tail = c >= tot - c;
    let cmin = c.min(tot - c);
    let freq = inversion_counts(n, cmin);
    let total = factorial_big(n);
    let left = Q::new(sum_big(&freq), total.clone());
    let at = Q::new(freq[cmin].clone(), total);
    let p = match alt {
        Alternative::TwoSided => (left * qi(2)).min(Q::one()),
        Alternative::Greater | Alternative::Less => {
            if in_right_tail == (alt == Alternative::Greater) {
                left
            } else {
                Q::one() - left + at
            }
        }
    };
    Ok(result(ctx, tau, p, None, alt))
}

/// Kolmogorov's distribution: `P(K > x) = 2 Σ_{k≥1} (−1)^{k−1} e^{−2k²x²}`.
fn kolmogorov_sf(x: f64) -> f64 {
    if x <= 0.0 {
        return 1.0;
    }
    if x < 1.0 {
        // Jacobi form: K(x) = √(2π)/x · Σ_{k≥1} exp(−(2k−1)²π²/(8x²)).
        let mut s = 0.0;
        for k in 1..=20 {
            let m = f64::from(2 * k - 1);
            s += (-m * m * std::f64::consts::PI * std::f64::consts::PI / (8.0 * x * x)).exp();
        }
        1.0 - (2.0 * std::f64::consts::PI).sqrt() / x * s
    } else {
        let mut s = 0.0;
        for k in 1..=200 {
            let kf = f64::from(k);
            let term = (-2.0 * kf * kf * x * x).exp();
            s += if k % 2 == 1 { term } else { -term };
            if term < 1e-18 {
                break;
            }
        }
        (2.0 * s).clamp(0.0, 1.0)
    }
}

/// Smirnov's exact one-sided distribution (Birnbaum–Tingey):
/// `P(D⁺ₙ ≥ d) = d Σ_{j=0}^{⌊n(1−d)⌋} C(n, j) (1 − d − j/n)^{n−j} (d + j/n)^{j−1}`,
/// summed in log space.  `scipy.stats.ksone.sf(d, n)`.
fn smirnov_sf(n: usize, d: f64) -> f64 {
    if d <= 0.0 {
        return 1.0;
    }
    if d >= 1.0 {
        return 0.0;
    }
    let nf = n as f64;
    let log_n_fact = lgamma(nf + 1.0);
    let mut sum = 0.0;
    for j in 0..=n {
        let jf = j as f64;
        let a = 1.0 - d - jf / nf;
        if a <= 0.0 {
            break;
        }
        let log_term = if j == 0 {
            nf * a.ln() - d.ln()
        } else {
            log_n_fact - lgamma(jf + 1.0) - lgamma(nf - jf + 1.0)
                + (nf - jf) * a.ln()
                + (jf - 1.0) * (d + jf / nf).ln()
        };
        sum += log_term.exp();
    }
    (d * sum).clamp(0.0, 1.0)
}

/// One-sample Kolmogorov–Smirnov test of `x` against a continuous
/// distribution.  `D⁺ = max(i/n − F(x₍ᵢ₎))`, `D⁻ = max(F(x₍ᵢ₎) − (i−1)/n)`,
/// `D = max(D⁺, D⁻)`, with `F` evaluated numerically (the exact CDF
/// expression of `dist` at each observation, then `eval_f64`).  The
/// two-sided p-value is asymptotic — Kolmogorov's distribution at `√n D`;
/// the one-sided p-values are Smirnov's exact `P(D⁺ₙ ≥ d)` (as scipy does
/// for every `method`).  Numerical throughout, because `F` is
/// transcendental for the distributions this test is used with.
/// `scipy.stats.ks_1samp(x, cdf, alternative, method='asymp')`.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::stats::Distribution;
/// use symplex::stats::data::from_f64;
/// use symplex::stats::hypothesis::{ks_one_sample, Alternative};
///
/// let ctx = Context::new();
/// let x = from_f64(&[-1.2, -0.3, 0.1, 0.4, 0.9, 1.5, 2.2, -0.7])?;
/// let normal = Distribution::normal(ctx.int(0), ctx.int(1));
/// // scipy: ks_1samp(x, norm.cdf, method='asymp') → statistic 0.19093987465324047, pvalue 0.9324475218943081
/// let r = ks_one_sample(&x, &normal, Alternative::TwoSided)?;
/// assert!((r.statistic - 0.190_939_874_653_240_47).abs() < 1e-12);
/// assert!((r.p_value - 0.932_447_521_894_308_1).abs() < 1e-9);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for an empty sample or a discrete
/// distribution; the CDF's evaluation error if a parameter is symbolic.
pub fn ks_one_sample(
    x: &[Q],
    dist: &Distribution,
    alt: Alternative,
) -> Result<KsResult, SymplexError> {
    const OP: &str = "ks_one_sample";
    check_sample(OP, "the sample", x, 1)?;
    if !dist.is_continuous() {
        return Err(invalid(OP, "the reference distribution must be continuous"));
    }
    let ctx = dist.context();
    let sorted = data::sorted(x);
    let n = sorted.len() as f64;
    let (mut d_plus, mut d_minus) = (0.0f64, 0.0f64);
    for (i, v) in sorted.iter().enumerate() {
        let f = dist.cdf(&ex(&ctx, v)).eval_f64()?;
        d_plus = d_plus.max((i + 1) as f64 / n - f);
        d_minus = d_minus.max(f - i as f64 / n);
    }
    let (statistic, p_value) = match alt {
        Alternative::TwoSided => {
            let d = d_plus.max(d_minus);
            (d, kolmogorov_sf(n.sqrt() * d))
        }
        Alternative::Greater => (d_plus, smirnov_sf(sorted.len(), d_plus)),
        Alternative::Less => (d_minus, smirnov_sf(sorted.len(), d_minus)),
    };
    Ok(KsResult {
        statistic,
        p_value,
        alternative: alt,
    })
}

// ═══════════════════════════════════════════════════════════════════════════
// 5. Effect sizes
// ═══════════════════════════════════════════════════════════════════════════

/// Cohen's `d = (x̄ − ȳ) / s` as an exact expression, with `s²` the pooled
/// variance `((n₁−1)s₁² + (n₂−1)s₂²)/(n₁+n₂−2)` (`pooled`) or the plain
/// average `(s₁² + s₂²)/2` (Cohen's original definition for equal sizes).
///
/// ```
/// use symplex::prelude::*;
/// use symplex::stats::data::from_i64;
/// use symplex::stats::hypothesis::cohens_d;
///
/// let ctx = Context::new();
/// let x = from_i64(&[20, 22, 19, 20, 22, 20, 21]);
/// let y = from_i64(&[28, 32, 36, 24, 29, 32]);
/// // numpy: (mean(x) − mean(y)) / sqrt(((n1−1)var(x, ddof=1) + (n2−1)var(y, ddof=1))/(n1+n2−2)) = -3.3080302571461795
/// let d = cohens_d(&ctx, &x, &y, true)?;
/// assert!((d.eval_f64()? - -3.308_030_257_146_179_5).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for a sample with fewer than two
/// observations or two constant samples.
pub fn cohens_d(ctx: &Context, x: &[Q], y: &[Q], pooled: bool) -> Result<Ex, SymplexError> {
    const OP: &str = "cohens_d";
    check_sample(OP, "the first sample", x, 2)?;
    check_sample(OP, "the second sample", y, 2)?;
    let (n1, n2) = (x.len(), y.len());
    let (v1, v2) = (
        data::variance(x, Ddof::Sample)?,
        data::variance(y, Ddof::Sample)?,
    );
    let var = if pooled {
        (qu(n1 - 1) * &v1 + qu(n2 - 1) * &v2) / qu(n1 + n2 - 2)
    } else {
        (&v1 + &v2) / qi(2)
    };
    if var.is_zero() {
        return Err(invalid(OP, "both samples are constant (zero variance)"));
    }
    let num = data::mean(x)? - data::mean(y)?;
    Ok(RootRatio { num, var }.to_ex(ctx))
}

/// Hedges' `g = d · J` with the pooled Cohen's `d` and the small-sample
/// correction `J = 1 − 3/(4(n₁ + n₂) − 9)` (Hedges & Olkin's approximation
/// of `Γ(ν/2)/(√(ν/2) Γ((ν−1)/2))`), as an exact expression.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::stats::data::from_i64;
/// use symplex::stats::hypothesis::hedges_g;
///
/// let ctx = Context::new();
/// let x = from_i64(&[20, 22, 19, 20, 22, 20, 21]);
/// let y = from_i64(&[28, 32, 36, 24, 29, 32]);
/// // numpy: d · (1 − 3/(4·13 − 9)) = -3.3080302571461795 · (40/43) = -3.077237448508074
/// let g = hedges_g(&ctx, &x, &y)?;
/// assert!((g.eval_f64()? - -3.077_237_448_508_074).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// As [`cohens_d`].
pub fn hedges_g(ctx: &Context, x: &[Q], y: &[Q]) -> Result<Ex, SymplexError> {
    let d = cohens_d(ctx, x, y, true)?;
    let n = usize_to_i64("hedges_g", x.len() + y.len())?;
    let j = ctx.one() - ctx.rational(3, 4 * n - 9);
    Ok((d * j).simplify())
}

/// Glass's `Δ = (x̄ − ȳ) / s_y`, standardised by the second (control)
/// sample's standard deviation alone, as an exact expression.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::stats::data::from_i64;
/// use symplex::stats::hypothesis::glass_delta;
///
/// let ctx = Context::new();
/// let x = from_i64(&[20, 22, 19, 20, 22, 20, 21]);
/// let y = from_i64(&[28, 32, 36, 24, 29, 32]);
/// // numpy: (mean(x) − mean(y)) / std(y, ddof=1) = -2.329471985471971
/// let g = glass_delta(&ctx, &x, &y)?;
/// assert!((g.eval_f64()? - -2.329_471_985_471_971).abs() < 1e-12);
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for an empty `x`, fewer than two
/// observations in `y`, or a constant `y`.
pub fn glass_delta(ctx: &Context, x: &[Q], y: &[Q]) -> Result<Ex, SymplexError> {
    const OP: &str = "glass_delta";
    check_sample(OP, "the first sample", x, 1)?;
    check_sample(OP, "the control sample", y, 2)?;
    let var = data::variance(y, Ddof::Sample)?;
    if var.is_zero() {
        return Err(invalid(
            OP,
            "the control sample is constant (zero variance)",
        ));
    }
    let num = data::mean(x)? - data::mean(y)?;
    Ok(RootRatio { num, var }.to_ex(ctx))
}

/// The rank-biserial correlation of a Mann–Whitney `U₁` (the statistic for
/// the first sample of sizes `n₁`, `n₂`): `r = 2U₁/(n₁n₂) − 1 = P(X > Y) −
/// P(X < Y)`, which equals Cliff's δ.  (Some sources report `1 − 2U/(n₁n₂)`
/// with the smaller `U`, i.e. `|r|`.)
///
/// ```
/// use symplex::linprog::{q, qi};
/// use symplex::stats::hypothesis::rank_biserial;
///
/// // U₁ = 17 for sizes 5 and 4: r = 34/20 − 1 = 7/10
/// assert_eq!(rank_biserial(&qi(17), 5, 4)?, q(7, 10));
/// # Ok::<(), symplex::prelude::SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for an empty sample or `U ∉ [0, n₁n₂]`.
pub fn rank_biserial(u1: &Q, n1: usize, n2: usize) -> Result<Q, SymplexError> {
    const OP: &str = "rank_biserial";
    if n1 == 0 || n2 == 0 {
        return Err(invalid(OP, "both samples must be non-empty"));
    }
    if u1.is_negative() || *u1 > qu(n1 * n2) {
        return Err(invalid(OP, "U must lie in [0, n₁n₂]"));
    }
    Ok(qi(2) * u1 / qu(n1 * n2) - Q::one())
}

/// `η² = SS_between / SS_total` of `k` groups, exact (the same quantity as
/// [`AnovaResult::eta_squared`]).
///
/// ```
/// use symplex::linprog::q;
/// use symplex::stats::data::from_i64;
/// use symplex::stats::hypothesis::eta_squared;
///
/// let g = [from_i64(&[6, 8, 4, 5, 3, 4]), from_i64(&[8, 12, 9, 11, 6, 8]), from_i64(&[13, 9, 11, 8, 7, 12])];
/// // SS_between = 84, SS_within = 68
/// assert_eq!(eta_squared(&g)?, q(84, 152));
/// # Ok::<(), symplex::prelude::SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for fewer than two groups, an empty
/// group, or identical observations throughout.
pub fn eta_squared(groups: &[Vec<Q>]) -> Result<Q, SymplexError> {
    const OP: &str = "eta_squared";
    let (ss_between, ss_within, _, _) = sums_of_squares(OP, groups)?;
    let total = &ss_between + &ss_within;
    if total.is_zero() {
        return Err(invalid(OP, "every observation is identical"));
    }
    Ok(ss_between / total)
}

/// Cliff's `δ = (#{xᵢ > yⱼ} − #{xᵢ < yⱼ}) / (n₁n₂)`, exact.
///
/// ```
/// use symplex::linprog::q;
/// use symplex::stats::data::from_i64;
/// use symplex::stats::hypothesis::cliffs_delta;
///
/// let x = from_i64(&[19, 22, 16, 29, 24]);
/// let y = from_i64(&[20, 11, 17, 12]);
/// assert_eq!(cliffs_delta(&x, &y)?, q(7, 10));
/// # Ok::<(), symplex::prelude::SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for an empty sample.
pub fn cliffs_delta(x: &[Q], y: &[Q]) -> Result<Q, SymplexError> {
    const OP: &str = "cliffs_delta";
    check_sample(OP, "the first sample", x, 1)?;
    check_sample(OP, "the second sample", y, 1)?;
    let mut diff = 0i64;
    for a in x {
        for b in y {
            diff += match a.cmp(b) {
                Ordering::Greater => 1,
                Ordering::Less => -1,
                Ordering::Equal => 0,
            };
        }
    }
    Ok(qi(diff) / qu(x.len() * y.len()))
}

// ═══════════════════════════════════════════════════════════════════════════
// 6. Multiple comparisons
// ═══════════════════════════════════════════════════════════════════════════

fn check_pvalues(op: &'static str, p: &[f64], alpha: f64) -> Result<(), SymplexError> {
    if p.is_empty() {
        return Err(invalid(op, "no p-values"));
    }
    if let Some(bad) = p.iter().find(|&&v| !(0.0..=1.0).contains(&v)) {
        return Err(invalid(
            op,
            format!("p-values must lie in [0, 1], got {bad}"),
        ));
    }
    check_unit_open(op, "alpha", alpha)
}

/// Indices that sort `p` ascending (stable).
fn ascending_order(p: &[f64]) -> Vec<usize> {
    let mut idx: Vec<usize> = (0..p.len()).collect();
    idx.sort_by(|&a, &b| p[a].total_cmp(&p[b]));
    idx
}

/// Scatter sorted results back into the input order.
fn unsort(order: &[usize], sorted_p: Vec<f64>, sorted_reject: Vec<bool>) -> Adjusted {
    let mut p_adjusted = vec![0.0; order.len()];
    let mut reject = vec![false; order.len()];
    for (rank, &i) in order.iter().enumerate() {
        p_adjusted[i] = sorted_p[rank].min(1.0);
        reject[i] = sorted_reject[rank];
    }
    Adjusted { p_adjusted, reject }
}

/// Bonferroni: `p̃ᵢ = min(1, m·pᵢ)`, reject where `p̃ᵢ ≤ α`.
/// `statsmodels.stats.multitest.multipletests(p, alpha, method='bonferroni')`.
///
/// ```
/// use symplex::stats::hypothesis::bonferroni;
///
/// // statsmodels: multipletests([0.01, 0.04, 0.03, 0.2], method='bonferroni')
/// //   → pvals_corrected [0.04, 0.16, 0.12, 0.8], reject [True, False, False, False]
/// let a = bonferroni(&[0.01, 0.04, 0.03, 0.2], 0.05)?;
/// assert!((a.p_adjusted[1] - 0.16).abs() < 1e-15);
/// assert_eq!(a.reject, vec![true, false, false, false]);
/// # Ok::<(), symplex::prelude::SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for no p-values, a p-value outside
/// `[0, 1]`, or `alpha ∉ (0, 1)`.
pub fn bonferroni(p: &[f64], alpha: f64) -> Result<Adjusted, SymplexError> {
    check_pvalues("bonferroni", p, alpha)?;
    let m = p.len() as f64;
    let p_adjusted: Vec<f64> = p.iter().map(|&v| (v * m).min(1.0)).collect();
    let reject = p_adjusted.iter().map(|&v| v <= alpha).collect();
    Ok(Adjusted { p_adjusted, reject })
}

/// Holm's step-down procedure: with `p₍₁₎ ≤ … ≤ p₍ₘ₎`, `p̃₍ᵢ₎ = min(1,
/// max_{j ≤ i} (m − j + 1) p₍ⱼ₎)`, reject where `p̃ ≤ α`.
/// `multipletests(p, alpha, method='holm')`.
///
/// ```
/// use symplex::stats::hypothesis::holm;
///
/// // statsmodels: multipletests([0.01, 0.04, 0.03, 0.2], method='holm')
/// //   → pvals_corrected [0.04, 0.09, 0.09, 0.2], reject [True, False, False, False]
/// let a = holm(&[0.01, 0.04, 0.03, 0.2], 0.05)?;
/// assert!((a.p_adjusted[2] - 0.09).abs() < 1e-15);
/// assert_eq!(a.reject, vec![true, false, false, false]);
/// # Ok::<(), symplex::prelude::SymplexError>(())
/// ```
///
/// # Errors
///
/// As [`bonferroni`].
pub fn holm(p: &[f64], alpha: f64) -> Result<Adjusted, SymplexError> {
    check_pvalues("holm", p, alpha)?;
    let m = p.len();
    let order = ascending_order(p);
    let mut sorted_p = Vec::with_capacity(m);
    let mut running = 0.0f64;
    for (rank, &i) in order.iter().enumerate() {
        running = running.max(p[i] * (m - rank) as f64);
        sorted_p.push(running);
    }
    let sorted_reject = sorted_p.iter().map(|&v| v.min(1.0) <= alpha).collect();
    Ok(unsort(&order, sorted_p, sorted_reject))
}

/// Step-up false-discovery-rate control with the critical constants
/// `cᵢ = i/(m·scale)`: `p̃₍ᵢ₎ = min(1, min_{j ≥ i} p₍ⱼ₎/cⱼ)`; reject `p₍ᵢ₎` up to
/// the largest `i` with `p₍ᵢ₎ ≤ α cᵢ` (statsmodels' `fdr_bh` / `fdr_by`).
fn fdr_step_up(
    op: &'static str,
    p: &[f64],
    alpha: f64,
    scale: f64,
) -> Result<Adjusted, SymplexError> {
    check_pvalues(op, p, alpha)?;
    let m = p.len();
    let order = ascending_order(p);
    let factor: Vec<f64> = (1..=m).map(|i| i as f64 / m as f64 / scale).collect();
    let mut sorted_p = vec![0.0; m];
    let mut running = f64::INFINITY;
    for rank in (0..m).rev() {
        running = running.min(p[order[rank]] / factor[rank]);
        sorted_p[rank] = running;
    }
    let last = (0..m)
        .rev()
        .find(|&rank| p[order[rank]] <= alpha * factor[rank]);
    let sorted_reject = (0..m).map(|rank| last.is_some_and(|l| rank <= l)).collect();
    Ok(unsort(&order, sorted_p, sorted_reject))
}

/// Benjamini–Hochberg: `p̃₍ᵢ₎ = min(1, min_{j ≥ i} m·p₍ⱼ₎/j)`; reject the
/// hypotheses up to the largest `i` with `p₍ᵢ₎ ≤ α i/m`.
/// `multipletests(p, alpha, method='fdr_bh')`.
///
/// ```
/// use symplex::stats::hypothesis::benjamini_hochberg;
///
/// // statsmodels: multipletests([0.01, 0.04, 0.03, 0.2], method='fdr_bh')
/// //   → pvals_corrected [0.04, 0.05333333333333334, 0.05333333333333334, 0.2], reject [True, False, False, False]
/// let a = benjamini_hochberg(&[0.01, 0.04, 0.03, 0.2], 0.05)?;
/// assert!((a.p_adjusted[1] - 0.053_333_333_333_333_34).abs() < 1e-15);
/// assert_eq!(a.reject, vec![true, false, false, false]);
/// # Ok::<(), symplex::prelude::SymplexError>(())
/// ```
///
/// # Errors
///
/// As [`bonferroni`].
pub fn benjamini_hochberg(p: &[f64], alpha: f64) -> Result<Adjusted, SymplexError> {
    fdr_step_up("benjamini_hochberg", p, alpha, 1.0)
}

/// Benjamini–Yekutieli (FDR under arbitrary dependence): Benjamini–Hochberg
/// with the constants divided by `Σ_{j=1}^{m} 1/j`.
/// `multipletests(p, alpha, method='fdr_by')`.
///
/// ```
/// use symplex::stats::hypothesis::benjamini_yekutieli;
///
/// // statsmodels: multipletests([0.01, 0.04, 0.03, 0.2], method='fdr_by')
/// //   → pvals_corrected [0.08333333333333331, 0.1111111111111111, 0.1111111111111111, 0.41666666666666663]
/// let a = benjamini_yekutieli(&[0.01, 0.04, 0.03, 0.2], 0.05)?;
/// assert!((a.p_adjusted[0] - 0.083_333_333_333_333_31).abs() < 1e-15);
/// assert_eq!(a.reject, vec![false; 4]);
/// # Ok::<(), symplex::prelude::SymplexError>(())
/// ```
///
/// # Errors
///
/// As [`bonferroni`].
pub fn benjamini_yekutieli(p: &[f64], alpha: f64) -> Result<Adjusted, SymplexError> {
    let harmonic: f64 = (1..=p.len()).map(|j| 1.0 / j as f64).sum();
    fdr_step_up("benjamini_yekutieli", p, alpha, harmonic)
}

// ═══════════════════════════════════════════════════════════════════════════
// 7. Resampling
// ═══════════════════════════════════════════════════════════════════════════

/// The `p`-quantile of an ascending sample by linear interpolation
/// (`numpy.quantile`'s default).
fn quantile_sorted(sorted: &[f64], p: f64) -> f64 {
    let n = sorted.len();
    if n == 1 {
        return sorted[0];
    }
    let h = (n - 1) as f64 * p;
    let lo = (h.floor() as usize).min(n - 1);
    let hi = (lo + 1).min(n - 1);
    let frac = h - lo as f64;
    sorted[lo] + frac * (sorted[hi] - sorted[lo])
}

fn check_f64_data(
    op: &'static str,
    name: &str,
    data: &[f64],
    min: usize,
) -> Result<(), SymplexError> {
    if data.len() < min {
        return Err(invalid(
            op,
            format!("{name} needs at least {min} observations"),
        ));
    }
    if let Some(bad) = data.iter().find(|v| !v.is_finite()) {
        return Err(invalid(
            op,
            format!("{name} contains a non-finite value {bad}"),
        ));
    }
    Ok(())
}

/// A bootstrap confidence interval for `statistic(data)`: `n_resamples`
/// resamples with replacement drawn with the deterministic `rng`, then the
/// `Percentile` interval `(q_{α/2}, q_{1−α/2})` of the resampled statistics
/// or the `Basic` interval `(2θ̂ − q_{1−α/2}, 2θ̂ − q_{α/2})`.  Quantiles are
/// linearly interpolated.  (`scipy.stats.bootstrap(method='percentile' |
/// 'basic')` up to the random stream.)
///
/// ```
/// use symplex::stats::Rng;
/// use symplex::stats::hypothesis::{bootstrap_ci, BootstrapMethod};
///
/// let data: Vec<f64> = (1..=20).map(f64::from).collect();
/// let mean = |x: &[f64]| x.iter().sum::<f64>() / x.len() as f64;
/// let (lo, hi) = bootstrap_ci(&data, mean, 2000, 0.95, &mut Rng::new(7), BootstrapMethod::Percentile)?;
/// assert!(lo < 10.5 && 10.5 < hi);
/// # Ok::<(), symplex::prelude::SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for empty or non-finite data,
/// `n_resamples = 0`, or `confidence ∉ (0, 1)`;
/// [`SymplexError::ComputationFailed`] if the statistic is non-finite.
pub fn bootstrap_ci(
    data: &[f64],
    statistic: impl Fn(&[f64]) -> f64,
    n_resamples: usize,
    confidence: f64,
    rng: &mut Rng,
    method: BootstrapMethod,
) -> Result<(f64, f64), SymplexError> {
    const OP: &str = "bootstrap_ci";
    check_f64_data(OP, "the data", data, 1)?;
    check_unit_open(OP, "confidence", confidence)?;
    if n_resamples == 0 {
        return Err(invalid(OP, "at least one resample is needed"));
    }
    let n = data.len();
    let observed = statistic(data);
    let mut resample = vec![0.0; n];
    let mut stats = Vec::with_capacity(n_resamples);
    for _ in 0..n_resamples {
        for slot in &mut resample {
            *slot = data[rng.below(n)];
        }
        let s = statistic(&resample);
        if !s.is_finite() {
            return Err(SymplexError::computation_failed(
                OP,
                "the statistic of a resample is not finite",
            ));
        }
        stats.push(s);
    }
    stats.sort_by(f64::total_cmp);
    let alpha = 1.0 - confidence;
    let (lo, hi) = (
        quantile_sorted(&stats, alpha / 2.0),
        quantile_sorted(&stats, 1.0 - alpha / 2.0),
    );
    Ok(match method {
        BootstrapMethod::Percentile => (lo, hi),
        BootstrapMethod::Basic => (2.0 * observed - hi, 2.0 * observed - lo),
    })
}

/// A randomised two-sample permutation test of `statistic(x, y)`: the
/// pooled sample is shuffled `n_permutations` times with the deterministic
/// `rng` and split into the original sizes; the p-value is `(#{permuted
/// statistics at least as extreme} + 1) / (n_permutations + 1)` — `T* ≥ T`
/// (`Greater`), `T* ≤ T` (`Less`), `2 min(·, ·)` clipped to `1` two-sided
/// (`scipy.stats.permutation_test(permutation_type='independent')`'s
/// definitions, up to the random stream).
///
/// ```
/// use symplex::stats::Rng;
/// use symplex::stats::hypothesis::{permutation_test, Alternative};
///
/// let x = [1.1, 2.3, 1.9, 2.8, 2.2, 1.7];
/// let y = [4.9, 5.2, 4.4, 5.8, 5.1, 4.7];
/// let diff = |a: &[f64], b: &[f64]| {
///     a.iter().sum::<f64>() / a.len() as f64 - b.iter().sum::<f64>() / b.len() as f64
/// };
/// let r = permutation_test(&x, &y, diff, 4000, &mut Rng::new(1), Alternative::TwoSided)?;
/// assert!(r.p_value < 0.01);
/// # Ok::<(), symplex::prelude::SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for an empty or non-finite sample or
/// `n_permutations = 0`; [`SymplexError::ComputationFailed`] if the
/// statistic is non-finite.
pub fn permutation_test(
    x: &[f64],
    y: &[f64],
    statistic: impl Fn(&[f64], &[f64]) -> f64,
    n_permutations: usize,
    rng: &mut Rng,
    alt: Alternative,
) -> Result<PermutationResult, SymplexError> {
    const OP: &str = "permutation_test";
    check_f64_data(OP, "the first sample", x, 1)?;
    check_f64_data(OP, "the second sample", y, 1)?;
    if n_permutations == 0 {
        return Err(invalid(OP, "at least one permutation is needed"));
    }
    let observed = statistic(x, y);
    if !observed.is_finite() {
        return Err(SymplexError::computation_failed(
            OP,
            "the observed statistic is not finite",
        ));
    }
    let n1 = x.len();
    let mut pooled: Vec<f64> = x.iter().chain(y).copied().collect();
    let slack = 1e-14 * observed.abs();
    let (mut count_ge, mut count_le) = (0usize, 0usize);
    for _ in 0..n_permutations {
        // Fisher–Yates shuffle.
        for i in (1..pooled.len()).rev() {
            let j = rng.below(i + 1);
            pooled.swap(i, j);
        }
        let s = statistic(&pooled[..n1], &pooled[n1..]);
        if !s.is_finite() {
            return Err(SymplexError::computation_failed(
                OP,
                "the statistic of a permutation is not finite",
            ));
        }
        if s >= observed - slack {
            count_ge += 1;
        }
        if s <= observed + slack {
            count_le += 1;
        }
    }
    let denom = (n_permutations + 1) as f64;
    let greater = (count_ge + 1) as f64 / denom;
    let less = (count_le + 1) as f64 / denom;
    let p_value = match alt {
        Alternative::Greater => greater,
        Alternative::Less => less,
        Alternative::TwoSided => (2.0 * greater.min(less)).min(1.0),
    };
    Ok(PermutationResult {
        statistic: observed,
        p_value,
    })
}

// ═══════════════════════════════════════════════════════════════════════════
// 8. Power and sample size
// ═══════════════════════════════════════════════════════════════════════════

/// The sample size for estimating a proportion `p` to within `±margin` at
/// the given confidence: `n = ⌈z²_{1−α/2} p(1−p) / margin²⌉`.
/// `statsmodels.stats.proportion.samplesize_confint_proportion(p, margin, alpha)`
/// (which returns the un-rounded `n`).
///
/// ```
/// use symplex::stats::hypothesis::sample_size_for_proportion;
///
/// // statsmodels: samplesize_confint_proportion(0.5, 0.03) = 1067.0718946372576  → 1068
/// assert_eq!(sample_size_for_proportion(0.03, 0.95, 0.5)?, 1068);
/// # Ok::<(), symplex::prelude::SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for `margin ≤ 0`, `confidence ∉ (0,
/// 1)` or `p ∉ (0, 1)`.
pub fn sample_size_for_proportion(
    margin: f64,
    confidence: f64,
    p: f64,
) -> Result<usize, SymplexError> {
    const OP: &str = "sample_size_for_proportion";
    check_finite(OP, "margin", margin)?;
    if margin <= 0.0 {
        return Err(invalid(OP, "the margin must be positive"));
    }
    check_unit_open(OP, "confidence", confidence)?;
    check_unit_open(OP, "p", p)?;
    let z = norm_isf((1.0 - confidence) / 2.0);
    let n = z * z * p * (1.0 - p) / (margin * margin);
    Ok(n.ceil() as usize)
}

/// statsmodels' two-sided normal power with `nobs = n/2` (equal groups):
/// `Φ̄(z_{α/2} − h√(n/2)) + Φ(−z_{α/2} − h√(n/2))`.
fn normal_power_two_sided(effect: f64, n_per_group: f64, alpha: f64) -> f64 {
    let crit = norm_isf(alpha / 2.0);
    let shift = effect * (n_per_group / 2.0).sqrt();
    norm_sf(crit - shift) + norm_cdf(-crit - shift)
}

fn check_proportions_and_alpha(
    op: &'static str,
    p1: f64,
    p2: f64,
    alpha: f64,
) -> Result<(), SymplexError> {
    check_unit_open(op, "p1", p1)?;
    check_unit_open(op, "p2", p2)?;
    check_unit_open(op, "alpha", alpha)
}

/// The power of the two-sided two-proportion z-test with `n_per_group` per
/// group at level `alpha`, by the normal approximation on Cohen's `h`:
/// `Φ̄(z_{α/2} − h√(n/2)) + Φ(−z_{α/2} − h√(n/2))`.
/// `statsmodels.stats.power.NormalIndPower().power(proportion_effectsize(p1, p2), n, alpha, ratio=1)`.
///
/// ```
/// use symplex::stats::hypothesis::power_two_proportions;
///
/// // statsmodels: NormalIndPower().power(proportion_effectsize(0.5, 0.4), nobs1=200, alpha=0.05, ratio=1) = 0.5214145419211713
/// let p = power_two_proportions(0.5, 0.4, 200, 0.05)?;
/// assert!((p - 0.521_414_541_921_171_3).abs() < 1e-9);
/// # Ok::<(), symplex::prelude::SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for proportions or `alpha` outside
/// `(0, 1)` or `n_per_group = 0`.
pub fn power_two_proportions(
    p1: f64,
    p2: f64,
    n_per_group: usize,
    alpha: f64,
) -> Result<f64, SymplexError> {
    const OP: &str = "power_two_proportions";
    check_proportions_and_alpha(OP, p1, p2, alpha)?;
    if n_per_group == 0 {
        return Err(invalid(OP, "the group size must be positive"));
    }
    let h = 2.0 * p1.sqrt().asin() - 2.0 * p2.sqrt().asin();
    Ok(normal_power_two_sided(h, n_per_group as f64, alpha))
}

/// The smallest integer `n` with `power(n) ≥ target` for an increasing
/// `power`, searched by doubling then bisection from `start`.
fn smallest_n_with_power(
    op: &'static str,
    start: usize,
    target: f64,
    power: impl Fn(usize) -> Result<f64, SymplexError>,
) -> Result<usize, SymplexError> {
    const CAP: usize = 1 << 40;
    let (mut lo, mut hi) = (start, start);
    while power(hi)? < target {
        lo = hi;
        hi = hi.saturating_mul(2);
        if hi > CAP {
            return Err(SymplexError::computation_failed(
                op,
                "the required sample size exceeds 2^40",
            ));
        }
    }
    // Invariant: power(hi) ≥ target; lo = start or power(lo) < target.
    if lo == hi {
        return Ok(hi);
    }
    while hi - lo > 1 {
        let mid = lo + (hi - lo) / 2;
        if power(mid)? >= target {
            hi = mid;
        } else {
            lo = mid;
        }
    }
    Ok(hi)
}

/// The per-group sample size for the two-sided two-proportion z-test to
/// reach `power` at level `alpha`: the smallest `n` with
/// [`power_two_proportions`]`(p1, p2, n, alpha) ≥ power`, i.e. the ceiling of
/// `statsmodels.stats.power.NormalIndPower().solve_power(proportion_effectsize(p1, p2), alpha=alpha, power=power, ratio=1)`.
///
/// ```
/// use symplex::stats::hypothesis::sample_size_two_proportions;
///
/// // statsmodels: NormalIndPower().solve_power(proportion_effectsize(0.5, 0.4), alpha=0.05, power=0.8, ratio=1) = 387.1677468578098
/// assert_eq!(sample_size_two_proportions(0.5, 0.4, 0.05, 0.8)?, 388);
/// # Ok::<(), symplex::prelude::SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for proportions, `alpha` or `power`
/// outside `(0, 1)`, `power ≤ alpha`, or `p1 = p2`.
pub fn sample_size_two_proportions(
    p1: f64,
    p2: f64,
    alpha: f64,
    power: f64,
) -> Result<usize, SymplexError> {
    const OP: &str = "sample_size_two_proportions";
    check_proportions_and_alpha(OP, p1, p2, alpha)?;
    check_unit_open(OP, "power", power)?;
    if power <= alpha {
        return Err(invalid(OP, "the target power must exceed alpha"));
    }
    if p1 == p2 {
        return Err(invalid(OP, "equal proportions have no finite sample size"));
    }
    smallest_n_with_power(OP, 1, power, |n| power_two_proportions(p1, p2, n, alpha))
}

/// The power of the two-sided two-sample Student t-test (equal group
/// sizes) for a standardised effect `d` at level `alpha`, through the
/// noncentral t distribution: with `ν = 2n − 2`, `δ = d√(n/2)` and the
/// critical `t_c = t_{1−α/2, ν}`,
/// `power = 1 − ∫₀^∞ [Φ(t_c√(v/ν) − δ) − Φ(−t_c√(v/ν) − δ)] f_{χ²_ν}(v) dv`,
/// integrated by adaptive Gauss–Kronrod quadrature (about `1e-10`).
/// `statsmodels.stats.power.TTestIndPower().power(d, nobs1=n, alpha=alpha, ratio=1)`.
///
/// ```
/// use symplex::stats::hypothesis::power_t_test_two_sample;
///
/// // statsmodels: TTestIndPower().power(0.5, nobs1=30, alpha=0.05, ratio=1) = 0.47789652076016464
/// let p = power_t_test_two_sample(0.5, 30, 0.05)?;
/// assert!((p - 0.477_896_520_760_164_64).abs() < 1e-8);
/// # Ok::<(), symplex::prelude::SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for a non-finite effect, `n < 2` or
/// `alpha ∉ (0, 1)`; the quantile's or quadrature's error if they fail.
pub fn power_t_test_two_sample(
    effect_size: f64,
    n_per_group: usize,
    alpha: f64,
) -> Result<f64, SymplexError> {
    const OP: &str = "power_t_test_two_sample";
    check_finite(OP, "effect_size", effect_size)?;
    check_unit_open(OP, "alpha", alpha)?;
    if n_per_group < 2 {
        return Err(invalid(OP, "each group needs at least two observations"));
    }
    let df = (2 * n_per_group - 2) as f64;
    let delta = effect_size * (n_per_group as f64 / 2.0).sqrt();
    let ctx = Context::new();
    let t_crit = student_t_quantile_f64(OP, &ctx, df, 1.0 - alpha / 2.0)?;
    // χ²_ν density in log form.
    let half = df / 2.0;
    let log_norm = half * std::f64::consts::LN_2 + lgamma(half);
    let density = move |v: f64| -> f64 {
        if v <= 0.0 {
            return 0.0;
        }
        ((half - 1.0) * v.ln() - v / 2.0 - log_norm).exp()
    };
    let integrand = move |v: f64| -> f64 {
        let scale = (v / df).sqrt();
        (norm_cdf(t_crit * scale - delta) - norm_cdf(-t_crit * scale - delta)) * density(v)
    };
    // The χ²_ν mass outside [ν − 40σ, ν + 40σ + 50] (σ = √(2ν)) is far below
    // double precision; a finite window keeps the adaptive rule on the peak.
    let sd = (2.0 * df).sqrt();
    let lo = (df - 40.0 * sd).max(0.0);
    let hi = df + 40.0 * sd + 50.0;
    let opts = QuadOpts::default();
    let (accept, _) = quadrature(&integrand, lo, hi, &opts)?;
    Ok((1.0 - accept).clamp(0.0, 1.0))
}

/// The per-group sample size for the two-sided two-sample t-test to reach
/// `power` at level `alpha`: the smallest `n ≥ 2` with
/// [`power_t_test_two_sample`]`(d, n, alpha) ≥ power`, i.e. the ceiling of
/// `statsmodels.stats.power.TTestIndPower().solve_power(d, alpha=alpha, power=power, ratio=1)`.
///
/// ```
/// use symplex::stats::hypothesis::sample_size_t_test_two_sample;
///
/// // statsmodels: TTestIndPower().solve_power(0.5, alpha=0.05, power=0.8, ratio=1) = 63.765610588911635  → 64
/// assert_eq!(sample_size_t_test_two_sample(0.5, 0.05, 0.8)?, 64);
/// # Ok::<(), symplex::prelude::SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] for `d = 0`, or `alpha` / `power`
/// outside `(0, 1)` with `power ≤ alpha`.
pub fn sample_size_t_test_two_sample(
    effect_size: f64,
    alpha: f64,
    power: f64,
) -> Result<usize, SymplexError> {
    const OP: &str = "sample_size_t_test_two_sample";
    check_finite(OP, "effect_size", effect_size)?;
    check_unit_open(OP, "alpha", alpha)?;
    check_unit_open(OP, "power", power)?;
    if power <= alpha {
        return Err(invalid(OP, "the target power must exceed alpha"));
    }
    if effect_size == 0.0 {
        return Err(invalid(OP, "a zero effect has no finite sample size"));
    }
    smallest_n_with_power(OP, 2, power, |n| {
        power_t_test_two_sample(effect_size, n, alpha)
    })
}