symplex 0.41.0

//! Order statistics: the distribution of the `k`-th smallest of `n`
//! independent draws from a distribution, as a [`Family`] of its own.
//!
//! For a continuous parent with distribution function `F` and density `f`
//! (David & Nagaraja, *Order Statistics*, §2.1),
//!
//! * density: `f_(k)(x) = k·C(n, k)·F(x)^{k−1}·(1 − F(x))^{n−k}·f(x)`,
//! * distribution function: `F_(k)(x) = I_{F(x)}(k, n − k + 1)`, the
//!   regularised incomplete beta function (`Σ_{j≥k} C(n, j) F^j (1−F)^{n−j}`),
//!
//! and the support is the parent's.  Moments go through the generic
//! machinery of [`Distribution`] (exact integration of the density), so
//! for `Uniform(0, 1)` the `k`-th of `n` has mean `k/(n + 1)` exactly.
//!
//! For a discrete parent on a finite set of values the order statistic is
//! enumerated exactly into a [`Finite`](super::Finite) table:
//! `P(X_(k) = v) = P(X_(k) ≤ v) − P(X_(k) ≤ v⁻)` with
//! `P(X_(k) ≤ v) = Σ_{j≥k} C(n, j) F(v)^j (1 − F(v))^{n−j}`.  On an
//! infinite integer lattice the same difference `I_{F(x)} − I_{F(x−1)}` is
//! the mass function of the [`OrderStatistic`] family.
//!
//! ```
//! use symplex::prelude::*;
//! use symplex::stats::Distribution;
//! use symplex::stats::order::{maximum_of, minimum_of};
//!
//! let ctx = Context::new();
//! let u = Distribution::uniform(ctx.int(0), ctx.int(1));
//! let m = maximum_of(&u, 3)?;                       // density 3x² on [0, 1]
//! assert_eq!(m.mean(), ctx.rational(3, 4));          // n/(n+1)
//! let d = minimum_of(&Distribution::die(ctx.int(6)), 2)?;
//! assert_eq!(d.density(&ctx.int(1)).eval(), ctx.rational(11, 36));
//! # Ok::<(), SymplexError>(())
//! ```

use std::fmt;

use crate::api::context::Context;
use crate::api::expr::Ex;
use crate::base::errors::SymplexError;
use crate::base::interval::{Interval, IntervalKind};

use super::data::binomial_q;
use super::family::{Distribution, Family, Sampler, TailOrders, family_boilerplate};
use super::sample::Rng;
use super::support::{Kind, Support, is_neg_inf, is_pos_inf};

const OP: &str = "stats::order";

fn invalid(reason: impl Into<String>) -> SymplexError {
    SymplexError::invalid_argument(OP, reason)
}

/// The largest finite discrete support enumerated into an exact table;
/// beyond it the lattice formula is used instead.
const MAX_ENUMERATED: i64 = 10_000;

/// The `k`-th smallest of `n` independent copies of `inner` (`1 ≤ k ≤ n`;
/// `k = 1` is the minimum, `k = n` the maximum).
#[derive(Clone, Debug, PartialEq)]
pub struct OrderStatistic {
    /// The parent distribution.
    pub inner: Distribution,
    /// Sample size `n ≥ 1`.
    pub n: usize,
    /// Rank `k`, `1 ≤ k ≤ n`.
    pub k: usize,
}

impl OrderStatistic {
    /// The family after validating `1 ≤ k ≤ n`.
    ///
    /// # Errors
    ///
    /// [`SymplexError::InvalidArgument`] if `n = 0` or `k ∉ [1, n]`.
    pub fn try_new(inner: Distribution, n: usize, k: usize) -> Result<Self, SymplexError> {
        if n == 0 {
            return Err(invalid("an order statistic needs a sample of size n ≥ 1"));
        }
        if k == 0 || k > n {
            return Err(invalid(format!(
                "the rank k must satisfy 1 ≤ k ≤ n, got k = {k}, n = {n}"
            )));
        }
        Ok(OrderStatistic { inner, n, k })
    }

    fn ctx(&self) -> Context {
        self.inner.context()
    }

    /// `n` and `k` as expressions.
    fn nk(&self) -> (Ex, Ex) {
        let ctx = self.ctx();
        (ctx.int(usize_to_i64(self.n)), ctx.int(usize_to_i64(self.k)))
    }

    /// The parent's `F(x)` as a single expression: `0` / `1` when `x` is
    /// decidably below / above the support, else the parent's
    /// [`cdf_on_support`](Distribution::cdf_on_support) — the closed form
    /// when the family has one, else the integral / sum of the density
    /// from the support's lower end.  (The whole-line [`Distribution::cdf`]
    /// would wrap a symbolic `x` in a `Piecewise`.)
    fn parent_cdf(&self, x: &Ex) -> Ex {
        let ctx = self.ctx();
        if let Some(iv) = self.inner.support().as_interval() {
            let (lo, hi) = (&iv.lower, &iv.upper);
            if !is_neg_inf(lo) && (x - lo).is_negative() == Some(true) {
                return ctx.zero();
            }
            if !is_pos_inf(hi) && (x - hi).is_nonnegative() == Some(true) {
                return ctx.one();
            }
        }
        self.inner.cdf_on_support(x)
    }

    /// The parent's `P(X > x)` as a single expression, as
    /// [`parent_cdf`](Self::parent_cdf): `1` / `0` when `x` is decidably
    /// below / above the support, else the parent's `sf_on_support` (its
    /// non-cancelling survival function in the far upper tail).
    fn parent_sf(&self, x: &Ex) -> Ex {
        let ctx = self.ctx();
        if let Some(iv) = self.inner.support().as_interval() {
            let (lo, hi) = (&iv.lower, &iv.upper);
            if !is_neg_inf(lo) && (x - lo).is_negative() == Some(true) {
                return ctx.one();
            }
            if !is_pos_inf(hi) && (x - hi).is_nonnegative() == Some(true) {
                return ctx.zero();
            }
        }
        self.inner.sf_on_support(x)
    }

    /// `I_u(k, n − k + 1)`: `P(X_(k) ≤ x)` as a function of `u = F(x)`.
    fn beta_cdf(&self, u: &Ex) -> Ex {
        let ctx = self.ctx();
        let (n, k) = self.nk();
        u.betainc_regularized(&k, &(n - &k + 1), &ctx.zero())
    }
}

fn usize_to_i64(n: usize) -> i64 {
    i64::try_from(n).unwrap_or(i64::MAX)
}

impl Family for OrderStatistic {
    family_boilerplate!(OrderStatistic, "OrderStatistic");

    fn context(&self) -> Context {
        self.ctx()
    }

    fn parameters(&self) -> Vec<(&'static str, Ex)> {
        let (n, k) = self.nk();
        let mut p = vec![("n", n), ("k", k)];
        p.extend(self.inner.parameters());
        p
    }

    fn support(&self) -> Support {
        self.inner.support()
    }

    // Continuous: k·C(n,k)·F^{k−1}·(1−F)^{n−k}·f.
    // Lattice: I_{F(x)}(k, n−k+1) − I_{F(x−1)}(k, n−k+1).
    fn density(&self, x: &Ex) -> Ex {
        let ctx = self.ctx();
        let f_x = self.parent_cdf(x);
        match self.inner.kind() {
            Kind::Continuous => {
                let (_, k) = self.nk();
                let choose = ctx.from_ratio(binomial_q(self.n, self.k));
                let mut acc = k * choose * self.inner.density(x);
                if self.k > 1 {
                    acc *= f_x.powi(usize_to_i64(self.k - 1));
                }
                if self.n > self.k {
                    acc *= (ctx.one() - &f_x).powi(usize_to_i64(self.n - self.k));
                }
                acc
            }
            Kind::Discrete => {
                let below = self.parent_cdf(&(x - 1));
                self.beta_cdf(&f_x) - self.beta_cdf(&below)
            }
        }
    }

    fn cdf(&self, x: &Ex) -> Option<Ex> {
        Some(self.beta_cdf(&self.parent_cdf(x)))
    }

    // P(X_(k) > x) = 1 − I_F(k, n − k + 1) = I_{1−F}(n − k + 1, k), with
    // 1 − F the parent's survival function: the maximum of n normals
    // beyond 20 is n·½ erfc(10√2) to first order, not 1 − (1 − …)ⁿ.
    fn sf(&self, x: &Ex) -> Option<Ex> {
        let ctx = self.ctx();
        let (n, k) = self.nk();
        Some(
            self.parent_sf(x)
                .betainc_regularized(&(n - &k + 1), &k, &ctx.zero()),
        )
    }

    // The tails of X_(k): P(X_(k) > x) ~ C(n, k−1) S(x)^{n−k+1} and
    // P(X_(k) ≤ x) ~ C(n, k) F(x)^k, so a parent's exponential rate scales
    // by n − k + 1 on the right and by k on the left (a heavy tail, end 0,
    // stays heavy).  A scaled finite end is taken as closed: whether the
    // order statistic's mgf is finite exactly there depends on the
    // parent's polynomial factor (the minimum of two Gamma(1/10) has a
    // finite mgf at t = 2/θ), so that one point keeps the integral.
    fn mgf_domain(&self) -> Option<Interval<Ex>> {
        let d = self.inner.family().mgf_domain()?;
        let (n, k) = self.nk();
        let m = &n - &k + 1;
        let scale = |v: Ex, factor: &Ex, open: bool| -> (Ex, bool) {
            if is_pos_inf(&v) || is_neg_inf(&v) || v.is_zero() == Some(true) {
                (v, open)
            } else {
                ((factor * v).simplify(), false)
            }
        };
        let (lower, lo_open) = scale(d.lower, &k, d.kind.lower_open());
        let (upper, hi_open) = scale(d.upper, &m, d.kind.upper_open());
        Some(Interval {
            lower,
            upper,
            kind: IntervalKind::from_open_ends(lo_open, hi_open),
        })
    }

    // The same tails scale the parent's orders: S(x)^{n−k+1} with
    // S(x) ~ x^{−r} has order (n−k+1)·r on the right, F(x)^k order k·r on
    // the left (the maximum of two Cauchy variables has a left tail of
    // order 2 and no mean only from the right).
    fn tail_orders(&self) -> Option<TailOrders> {
        let parent = self.inner.family().tail_orders()?;
        let (n, k) = self.nk();
        let m = &n - &k + 1;
        Some(TailOrders {
            left: parent.left.map(|r| (&k * r).simplify()),
            right: parent.right.map(|r| (m * r).simplify()),
        })
    }

    // n draws from the parent, sorted; the k-th smallest.
    fn sampler(&self) -> Option<Result<Sampler, SymplexError>> {
        let mut inner = match self.inner.sampler() {
            Ok(s) => s,
            Err(e) => return Some(Err(e)),
        };
        let (n, k) = (self.n, self.k);
        let mut buf = vec![0.0_f64; n];
        Some(Ok(Box::new(move |rng: &mut Rng| {
            for slot in buf.iter_mut() {
                *slot = inner(rng);
            }
            buf.sort_by(f64::total_cmp);
            buf.get(k.saturating_sub(1)).copied().unwrap_or(f64::NAN)
        })))
    }

    fn fmt_display(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(
            f,
            "OrderStatistic(k={}, n={}, {})",
            self.k, self.n, self.inner
        )
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// Constructors
// ═══════════════════════════════════════════════════════════════════════════

/// The distribution of the `k`-th smallest of `n` independent draws from
/// `dist`: an [`OrderStatistic`] family for a continuous parent (or a
/// parent on an infinite integer lattice), and an exactly enumerated
/// [`Finite`](super::Finite) table for a parent with finitely many values.
///
/// ```
/// use symplex::prelude::*;
/// use symplex::stats::Distribution;
/// use symplex::stats::order::order_statistic;
///
/// let ctx = Context::new();
/// let u = Distribution::uniform(ctx.int(0), ctx.int(1));
/// // The second smallest of five: mean 2/6, P(X_(2) ≤ 1/2) = I_{1/2}(2, 4) = 13/16
/// let x2 = order_statistic(&u, 5, 2)?;
/// assert_eq!(x2.mean(), ctx.rational(1, 3));
/// assert_eq!(x2.cdf(&ctx.rational(1, 2)), ctx.rational(13, 16));
/// # Ok::<(), SymplexError>(())
/// ```
///
/// # Errors
///
/// [`SymplexError::InvalidArgument`] if `n = 0` or `k ∉ [1, n]`, or a
/// finite parent has values that cannot be ordered numerically.
pub fn order_statistic(
    dist: &Distribution,
    n: usize,
    k: usize,
) -> Result<Distribution, SymplexError> {
    let family = OrderStatistic::try_new(dist.clone(), n, k)?;
    if dist.kind() == Kind::Discrete
        && let Some(table) = finite_table(dist)?
    {
        return enumerate(dist, &table, n, k);
    }
    Ok(Distribution::from_family(family))
}

/// The minimum of `n` independent draws (`k = 1`).
///
/// # Errors
///
/// As [`order_statistic`].
pub fn minimum_of(dist: &Distribution, n: usize) -> Result<Distribution, SymplexError> {
    order_statistic(dist, n, 1)
}

/// The maximum of `n` independent draws (`k = n`).
///
/// # Errors
///
/// As [`order_statistic`].
pub fn maximum_of(dist: &Distribution, n: usize) -> Result<Distribution, SymplexError> {
    order_statistic(dist, n, n)
}

/// The `(value, mass)` table of a discrete distribution with finitely many
/// values, sorted by value; `None` for an infinite (or huge) lattice.
fn finite_table(dist: &Distribution) -> Result<Option<Vec<(Ex, Ex)>>, SymplexError> {
    let support = dist.support();
    let values: Vec<Ex> = if let Some(points) = support.as_points() {
        points
    } else if let Some(iv) = support.as_interval() {
        let (Some(lo), Some(hi)) = (iv.lower.eval().as_i64(), iv.upper.eval().as_i64()) else {
            return Ok(None);
        };
        if hi < lo || hi - lo >= MAX_ENUMERATED {
            return Ok(None);
        }
        let ctx = dist.context();
        (lo..=hi).map(|v| ctx.int(v)).collect()
    } else {
        return Ok(None);
    };
    let mut keyed: Vec<(f64, Ex)> = Vec::with_capacity(values.len());
    for v in values {
        let key = v.eval_f64().map_err(|_| {
            invalid(format!(
                "the values of a finite parent must be numeric to order them, got `{v}`"
            ))
        })?;
        keyed.push((key, v));
    }
    keyed.sort_by(|a, b| a.0.total_cmp(&b.0));
    let table = keyed
        .into_iter()
        .map(|(_, v)| {
            let p = dist.density(&v).eval();
            (v, p)
        })
        .collect();
    Ok(Some(table))
}

/// The exact table of the `k`-th of `n` draws from a finite parent.
fn enumerate(
    dist: &Distribution,
    table: &[(Ex, Ex)],
    n: usize,
    k: usize,
) -> Result<Distribution, SymplexError> {
    let ctx = dist.context();
    // G(u) = Σ_{j≥k} C(n, j) uʲ (1−u)ⁿ⁻ʲ
    let g = |u: &Ex| -> Ex {
        let mut acc = ctx.zero();
        for j in k..=n {
            let mut term = ctx.from_ratio(binomial_q(n, j));
            if j > 0 {
                term *= u.powi(usize_to_i64(j));
            }
            if n > j {
                term *= (ctx.one() - u).powi(usize_to_i64(n - j));
            }
            acc += term;
        }
        acc.eval()
    };
    let mut cumulative = ctx.zero();
    let mut previous = ctx.zero();
    let mut out = Vec::with_capacity(table.len());
    for (v, p) in table {
        cumulative = (&cumulative + p).eval();
        let current = g(&cumulative);
        let mass = (&current - &previous).simplify();
        previous = current;
        out.push((v.clone(), mass));
    }
    Ok(Distribution::finite(&ctx, out))
}