u_analytics/

distribution.rs

//! Distribution analysis.
//!
//! Empirical distribution functions, histogram binning, QQ-plot data,
//! distribution fitting/testing, and kernel density estimation.
//!
//! # Examples
//!
//! ```
//! use u_analytics::distribution::{ecdf, histogram_bins, BinMethod};
//!
//! let data = [1.0, 2.0, 3.0, 4.0, 5.0];
//! let (values, probs) = ecdf(&data).unwrap();
//! assert_eq!(values.len(), 5);
//! assert!((probs[4] - 1.0).abs() < 1e-10);
//!
//! let bins = histogram_bins(&data, BinMethod::Sturges).unwrap();
//! assert!(bins.n_bins >= 2);
//! ```

use u_numflow::distributions::{
    BetaDistribution, Exponential, GammaDistribution, LogNormal, Normal, Sample, Weibull,
};
use u_numflow::special;
use u_numflow::stats;

/// A point on the empirical CDF.
///
/// Computes the ECDF: F_n(x) = (number of observations ≤ x) / n.
///
/// # Returns
///
/// Tuple of (sorted unique values, cumulative probabilities). Returns
/// `None` if the data is empty or contains non-finite values.
///
/// # Examples
///
/// ```
/// use u_analytics::distribution::ecdf;
///
/// let data = [3.0, 1.0, 2.0, 1.0, 4.0];
/// let (vals, probs) = ecdf(&data).unwrap();
/// assert_eq!(vals, vec![1.0, 2.0, 3.0, 4.0]);
/// assert!((probs[0] - 0.4).abs() < 1e-10); // 2 values ≤ 1.0
/// assert!((probs[3] - 1.0).abs() < 1e-10);
/// ```
pub fn ecdf(data: &[f64]) -> Option<(Vec<f64>, Vec<f64>)> {
    if data.is_empty() || data.iter().any(|v| !v.is_finite()) {
        return None;
    }

    let n = data.len() as f64;
    let mut sorted: Vec<f64> = data.to_vec();
    sorted.sort_by(|a, b| a.partial_cmp(b).expect("finite values"));

    let mut values = Vec::new();
    let mut probs = Vec::new();

    let mut i = 0;
    while i < sorted.len() {
        let val = sorted[i];
        // Count all occurrences of this value
        let mut j = i;
        while j < sorted.len() && (sorted[j] - val).abs() < 1e-300 {
            j += 1;
        }
        values.push(val);
        probs.push(j as f64 / n);
        i = j;
    }

    Some((values, probs))
}

// ---------------------------------------------------------------------------
// Histogram binning
// ---------------------------------------------------------------------------

/// Method for computing optimal number of histogram bins.
#[derive(Debug, Clone, Copy)]
pub enum BinMethod {
    /// Sturges' rule: k = ⌈log₂(n)⌉ + 1. Best for near-normal data.
    Sturges,
    /// Scott's rule: h = 3.49·σ·n^(-1/3). Width-based.
    Scott,
    /// Freedman-Diaconis rule: h = 2·IQR·n^(-1/3). Robust to outliers.
    FreedmanDiaconis,
    /// Explicit bin count chosen by the caller. Must be ≥ 1.
    Fixed(usize),
}

/// Result of histogram bin computation.
#[derive(Debug, Clone)]
pub struct HistogramBins {
    /// Number of bins.
    pub n_bins: usize,
    /// Bin width.
    pub bin_width: f64,
    /// Bin edges (length = n_bins + 1).
    pub edges: Vec<f64>,
    /// Bin counts.
    pub counts: Vec<usize>,
}

/// Computes optimal histogram bins using the specified method.
///
/// # Returns
///
/// `None` if fewer than 2 data points, non-finite values, or zero range.
///
/// # Examples
///
/// ```
/// use u_analytics::distribution::{histogram_bins, BinMethod};
///
/// let data = [1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0];
/// let result = histogram_bins(&data, BinMethod::Sturges).unwrap();
/// assert!(result.n_bins >= 3);
/// assert_eq!(result.edges.len(), result.n_bins + 1);
/// assert_eq!(result.counts.len(), result.n_bins);
/// ```
pub fn histogram_bins(data: &[f64], method: BinMethod) -> Option<HistogramBins> {
    let n = data.len();
    if n < 2 || data.iter().any(|v| !v.is_finite()) {
        return None;
    }

    let min_val = data.iter().cloned().reduce(f64::min)?;
    let max_val = data.iter().cloned().reduce(f64::max)?;
    let range = max_val - min_val;

    if range < 1e-300 {
        return None; // all same value
    }

    let nf = n as f64;

    let n_bins = match method {
        BinMethod::Sturges => {
            let k = (nf.log2()).ceil() as usize + 1;
            k.max(2)
        }
        BinMethod::Scott => {
            let sd = stats::std_dev(data)?;
            if sd < 1e-300 {
                return None;
            }
            let h = 3.49 * sd * nf.powf(-1.0 / 3.0);
            ((range / h).ceil() as usize).max(2)
        }
        BinMethod::FreedmanDiaconis => {
            let q1 = stats::quantile(data, 0.25)?;
            let q3 = stats::quantile(data, 0.75)?;
            let iqr = q3 - q1;
            if iqr < 1e-300 {
                // Fall back to Sturges if IQR is zero
                let k = (nf.log2()).ceil() as usize + 1;
                k.max(2)
            } else {
                let h = 2.0 * iqr * nf.powf(-1.0 / 3.0);
                ((range / h).ceil() as usize).max(2)
            }
        }
        // Caller-specified count is honored exactly (no automatic floor);
        // zero bins is invalid input.
        BinMethod::Fixed(0) => return None,
        BinMethod::Fixed(k) => k,
    };

    let bin_width = range / n_bins as f64;

    // Compute edges
    let mut edges = Vec::with_capacity(n_bins + 1);
    for i in 0..=n_bins {
        edges.push(min_val + i as f64 * bin_width);
    }

    // Count data in each bin
    let mut counts = vec![0_usize; n_bins];
    for &x in data {
        let bin = ((x - min_val) / bin_width).floor() as usize;
        let bin = bin.min(n_bins - 1); // last point goes in last bin
        counts[bin] += 1;
    }

    Some(HistogramBins {
        n_bins,
        bin_width,
        edges,
        counts,
    })
}

// ---------------------------------------------------------------------------
// QQ-plot
// ---------------------------------------------------------------------------

/// Generates QQ-plot data (theoretical quantiles vs sample quantiles).
///
/// Uses the standard normal distribution as the theoretical distribution.
/// Sample quantiles are the sorted data. Theoretical quantiles are computed
/// using the plotting position (i - 0.5) / n.
///
/// # Returns
///
/// Tuple of (theoretical_quantiles, sample_quantiles). Returns `None` if
/// fewer than 3 data points.
///
/// # Examples
///
/// ```
/// use u_analytics::distribution::qq_plot_normal;
///
/// let data = [-1.5, -0.5, 0.0, 0.5, 1.5];
/// let (theoretical, sample) = qq_plot_normal(&data).unwrap();
/// assert_eq!(theoretical.len(), 5);
/// assert_eq!(sample.len(), 5);
/// // Sample and theoretical should be roughly aligned for normal data
/// ```
pub fn qq_plot_normal(data: &[f64]) -> Option<(Vec<f64>, Vec<f64>)> {
    let n = data.len();
    if n < 3 || data.iter().any(|v| !v.is_finite()) {
        return None;
    }

    let mut sample: Vec<f64> = data.to_vec();
    sample.sort_by(|a, b| a.partial_cmp(b).expect("finite values"));

    let nf = n as f64;
    let theoretical: Vec<f64> = (0..n)
        .map(|i| {
            let p = (i as f64 + 0.5) / nf;
            special::inverse_normal_cdf(p)
        })
        .collect();

    Some((theoretical, sample))
}

// ---------------------------------------------------------------------------
// Kolmogorov-Smirnov test
// ---------------------------------------------------------------------------

/// Kolmogorov-Smirnov one-sample test against the standard normal distribution.
///
/// # Algorithm
///
/// D = max|F_n(x) - Φ(x)| where F_n is the ECDF and Φ is the CDF of N(0,1).
/// P-value approximation: Kolmogorov distribution (Marsaglia et al., 2003).
///
/// # Returns
///
/// `None` if fewer than 5 observations or non-finite values.
///
/// # Examples
///
/// ```
/// use u_analytics::distribution::ks_test_normal;
///
/// // Data from approximately N(0,1)
/// let data = [-1.2, -0.8, -0.3, 0.1, 0.5, 0.7, 1.1, 1.4];
/// let (d_stat, p_value) = ks_test_normal(&data).unwrap();
/// assert!(p_value > 0.05); // cannot reject normality
/// ```
pub fn ks_test_normal(data: &[f64]) -> Option<(f64, f64)> {
    let n = data.len();
    if n < 5 || data.iter().any(|v| !v.is_finite()) {
        return None;
    }

    let mean = stats::mean(data)?;
    let sd = stats::std_dev(data)?;
    if sd < 1e-300 {
        return None;
    }

    let mut sorted: Vec<f64> = data.to_vec();
    sorted.sort_by(|a, b| a.partial_cmp(b).expect("finite values"));

    let nf = n as f64;
    let mut d_stat = 0.0_f64;

    for (i, &x) in sorted.iter().enumerate() {
        let z = (x - mean) / sd;
        let cdf = special::standard_normal_cdf(z);
        let ecdf_above = (i + 1) as f64 / nf;
        let ecdf_below = i as f64 / nf;
        d_stat = d_stat.max((ecdf_above - cdf).abs());
        d_stat = d_stat.max((ecdf_below - cdf).abs());
    }

    // P-value approximation (Kolmogorov distribution, simplified)
    // For large n: P(D > x) ≈ 2·Σ (-1)^(k-1) exp(-2k²n·x²)
    let lambda = (nf.sqrt() + 0.12 + 0.11 / nf.sqrt()) * d_stat;
    let mut p_value = 0.0;
    for k in 1..=100 {
        let kf = k as f64;
        let sign = if k % 2 == 1 { 1.0 } else { -1.0 };
        let term = sign * (-2.0 * kf * kf * lambda * lambda).exp();
        p_value += term;
        if term.abs() < 1e-15 {
            break;
        }
    }
    p_value = (2.0 * p_value).clamp(0.0, 1.0);

    Some((d_stat, p_value))
}

// ---------------------------------------------------------------------------
// Kernel Density Estimation
// ---------------------------------------------------------------------------

/// Bandwidth selection method for kernel density estimation.
#[derive(Debug, Clone, Copy)]
pub enum BandwidthMethod {
    /// Silverman's rule of thumb: h = 0.9 * min(σ, IQR/1.34) * n^(-1/5).
    /// Robust to outliers and multimodal distributions.
    ///
    /// Reference: Silverman (1986), "Density Estimation for Statistics and
    /// Data Analysis"
    Silverman,
    /// Scott's rule: h = 1.06 * σ * n^(-1/5).
    /// Slightly smoother, assumes approximately normal data.
    ///
    /// Reference: Scott (1992), "Multivariate Density Estimation"
    Scott,
    /// Manual bandwidth specification.
    Manual(f64),
}

/// Result of kernel density estimation.
#[derive(Debug, Clone)]
pub struct KdeResult {
    /// Evaluation points (x-axis).
    pub x: Vec<f64>,
    /// Density estimates at each evaluation point (y-axis).
    pub density: Vec<f64>,
    /// Bandwidth used.
    pub bandwidth: f64,
}

/// Gaussian kernel density estimation.
///
/// # Algorithm
///
/// f̂(x) = (1/nh) Σᵢ K((x - xᵢ)/h)
///
/// where K(u) = (1/√(2π)) exp(-u²/2) is the Gaussian kernel and
/// h is the bandwidth (smoothing parameter).
///
/// The evaluation grid extends 3h beyond the data range to capture tail
/// contributions from boundary points (3σ covers 99.7% of a Gaussian).
///
/// Reference: Silverman (1986), "Density Estimation for Statistics and
/// Data Analysis"
///
/// # Parameters
///
/// - `data`: Sample observations (must have ≥ 2 finite values)
/// - `method`: Bandwidth selection method
/// - `n_points`: Number of evaluation grid points (typical: 256–1024)
///
/// # Returns
///
/// `None` if fewer than 2 data points, fewer than 2 grid points,
/// non-finite values, or zero variance (for automatic bandwidth methods).
///
/// # Examples
///
/// ```
/// use u_analytics::distribution::{kde, BandwidthMethod};
///
/// let data = [1.0, 1.1, 1.2, 2.0, 2.1, 2.2, 5.0];
/// let result = kde(&data, BandwidthMethod::Silverman, 512).unwrap();
/// assert_eq!(result.x.len(), 512);
/// assert_eq!(result.density.len(), 512);
/// assert!(result.bandwidth > 0.0);
/// // Density should integrate approximately to 1
/// let dx = result.x[1] - result.x[0];
/// let integral: f64 = result.density.iter().sum::<f64>() * dx;
/// assert!((integral - 1.0).abs() < 0.05);
/// ```
pub fn kde(data: &[f64], method: BandwidthMethod, n_points: usize) -> Option<KdeResult> {
    let n = data.len();
    if n < 2 || n_points < 2 || data.iter().any(|v| !v.is_finite()) {
        return None;
    }

    let bandwidth = kde_bandwidth(data, method)?;

    // Evaluation grid: extend 3*bandwidth beyond data range
    let min_val = data.iter().cloned().reduce(f64::min)?;
    let max_val = data.iter().cloned().reduce(f64::max)?;
    let x_min = min_val - 3.0 * bandwidth;
    let x_max = max_val + 3.0 * bandwidth;
    let step = (x_max - x_min) / (n_points - 1) as f64;

    let x: Vec<f64> = (0..n_points).map(|i| x_min + i as f64 * step).collect();

    // Evaluate density at each grid point
    let inv_h = 1.0 / bandwidth;
    let inv_nh = inv_h / n as f64;
    let inv_sqrt_2pi = 1.0 / (2.0 * std::f64::consts::PI).sqrt();

    let density: Vec<f64> = x
        .iter()
        .map(|&xi| {
            let sum: f64 = data
                .iter()
                .map(|&xj| {
                    let u = (xi - xj) * inv_h;
                    inv_sqrt_2pi * (-0.5 * u * u).exp()
                })
                .sum();
            sum * inv_nh
        })
        .collect();

    Some(KdeResult {
        x,
        density,
        bandwidth,
    })
}

/// Evaluates the kernel density estimate at a single point.
///
/// Useful for point-wise evaluation without generating a full grid.
///
/// # Returns
///
/// `None` if data is empty, contains non-finite values, or bandwidth
/// is invalid.
pub fn kde_evaluate(data: &[f64], bandwidth: f64, x: f64) -> Option<f64> {
    let n = data.len();
    if n == 0 || bandwidth <= 0.0 || !bandwidth.is_finite() || !x.is_finite() {
        return None;
    }
    if data.iter().any(|v| !v.is_finite()) {
        return None;
    }

    let inv_h = 1.0 / bandwidth;
    let inv_nh = inv_h / n as f64;
    let inv_sqrt_2pi = 1.0 / (2.0 * std::f64::consts::PI).sqrt();

    let sum: f64 = data
        .iter()
        .map(|&xj| {
            let u = (x - xj) * inv_h;
            inv_sqrt_2pi * (-0.5 * u * u).exp()
        })
        .sum();

    Some(sum * inv_nh)
}

/// Computes the bandwidth for KDE using the specified method.
///
/// Useful when you want to inspect the bandwidth before running full KDE.
///
/// # Returns
///
/// `None` if fewer than 2 data points, non-finite values, or zero
/// variance (for automatic methods).
pub fn kde_bandwidth(data: &[f64], method: BandwidthMethod) -> Option<f64> {
    let n = data.len();
    if n < 2 || data.iter().any(|v| !v.is_finite()) {
        return None;
    }

    match method {
        BandwidthMethod::Silverman => {
            let sd = stats::std_dev(data)?;
            if sd < 1e-300 {
                return None;
            }
            let q1 = stats::quantile(data, 0.25)?;
            let q3 = stats::quantile(data, 0.75)?;
            let iqr = q3 - q1;
            let spread = if iqr > 1e-300 { sd.min(iqr / 1.34) } else { sd };
            Some(0.9 * spread * (n as f64).powf(-0.2))
        }
        BandwidthMethod::Scott => {
            let sd = stats::std_dev(data)?;
            if sd < 1e-300 {
                return None;
            }
            Some(1.06 * sd * (n as f64).powf(-0.2))
        }
        BandwidthMethod::Manual(h) => {
            if h <= 0.0 || !h.is_finite() {
                None
            } else {
                Some(h)
            }
        }
    }
}

// ---------------------------------------------------------------------------
// MLE Distribution Fitting
// ---------------------------------------------------------------------------

/// Result of a distribution fit via maximum likelihood estimation.
#[derive(Debug, Clone)]
pub struct FitResult {
    /// Name of the fitted distribution (e.g., "Normal", "Exponential").
    pub distribution: String,
    /// Estimated parameters (name, value) pairs.
    pub parameters: Vec<(String, f64)>,
    /// Log-likelihood at MLE parameters.
    pub log_likelihood: f64,
    /// Akaike Information Criterion: -2·ℓ + 2k.
    pub aic: f64,
    /// Bayesian Information Criterion: -2·ℓ + k·ln(n).
    pub bic: f64,
    /// Number of estimated parameters (k).
    pub n_params: usize,
}

/// Fits a Normal distribution N(μ, σ²) via MLE.
///
/// # Estimators
///
/// - μ̂ = x̄ (sample mean)
/// - σ̂ = √((1/n) Σ(xᵢ - x̄)²) (biased MLE, not sample std dev)
///
/// # Returns
///
/// `None` if fewer than 2 data points, non-finite values, or zero variance.
///
/// # Examples
///
/// ```
/// use u_analytics::distribution::fit_normal;
///
/// let data = [2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0];
/// let fit = fit_normal(&data).unwrap();
/// assert_eq!(fit.distribution, "Normal");
/// assert!((fit.parameters[0].1 - 5.0).abs() < 1e-10); // μ̂ = 5.0
/// assert!(fit.parameters[1].1 > 0.0); // σ̂ > 0
/// ```
pub fn fit_normal(data: &[f64]) -> Option<FitResult> {
    let n = data.len();
    if n < 2 || data.iter().any(|v| !v.is_finite()) {
        return None;
    }

    let nf = n as f64;
    let mu = stats::mean(data)?;

    // Biased MLE variance (denominator n, not n-1)
    let sum_sq: f64 = data.iter().map(|&x| (x - mu).powi(2)).sum();
    let sigma_sq = sum_sq / nf;
    if sigma_sq < 1e-300 {
        return None;
    }
    let sigma = sigma_sq.sqrt();

    // Log-likelihood: -n/2·ln(2π) - n·ln(σ) - n/2
    let log_lik = -nf / 2.0 * (2.0 * std::f64::consts::PI).ln() - nf * sigma.ln() - nf / 2.0;

    let k = 2;
    let aic = -2.0 * log_lik + 2.0 * k as f64;
    let bic = -2.0 * log_lik + k as f64 * nf.ln();

    Some(FitResult {
        distribution: "Normal".to_string(),
        parameters: vec![("mu".to_string(), mu), ("sigma".to_string(), sigma)],
        log_likelihood: log_lik,
        aic,
        bic,
        n_params: k,
    })
}

/// Fits an Exponential distribution Exp(λ) via MLE.
///
/// # Estimators
///
/// - λ̂ = 1/x̄ (rate parameter)
///
/// # Returns
///
/// `None` if fewer than 2 data points, non-positive values, or non-finite values.
///
/// # Examples
///
/// ```
/// use u_analytics::distribution::fit_exponential;
///
/// let data = [0.5, 1.2, 0.8, 2.1, 0.3, 1.5, 0.9, 1.8];
/// let fit = fit_exponential(&data).unwrap();
/// assert_eq!(fit.distribution, "Exponential");
/// assert!(fit.parameters[0].1 > 0.0); // λ̂ > 0
/// ```
pub fn fit_exponential(data: &[f64]) -> Option<FitResult> {
    let n = data.len();
    if n < 2 || data.iter().any(|v| !v.is_finite()) {
        return None;
    }
    // Exponential requires x > 0
    if data.iter().any(|&v| v <= 0.0) {
        return None;
    }

    let nf = n as f64;
    let mean = stats::mean(data)?;
    if mean < 1e-300 {
        return None;
    }

    let lambda = 1.0 / mean;

    // Log-likelihood: n·ln(λ) - λ·n·x̄ = n·ln(λ) - n
    let log_lik = nf * lambda.ln() - nf;

    let k = 1;
    let aic = -2.0 * log_lik + 2.0 * k as f64;
    let bic = -2.0 * log_lik + k as f64 * nf.ln();

    Some(FitResult {
        distribution: "Exponential".to_string(),
        parameters: vec![("lambda".to_string(), lambda)],
        log_likelihood: log_lik,
        aic,
        bic,
        n_params: k,
    })
}

/// Digamma function ψ(x) = d/dx ln(Γ(x)).
///
/// Uses the asymptotic expansion for x ≥ 8, with recurrence for smaller x.
/// Reference: Abramowitz & Stegun (1972), formula 6.3.18.
fn digamma(x: f64) -> f64 {
    if x <= 0.0 {
        return f64::NAN;
    }

    // Recurrence: ψ(x) = ψ(x+1) - 1/x
    let mut val = 0.0;
    let mut x = x;
    while x < 8.0 {
        val -= 1.0 / x;
        x += 1.0;
    }

    // Asymptotic expansion for x ≥ 8
    // ψ(x) ≈ ln(x) - 1/(2x) - 1/(12x²) + 1/(120x⁴) - 1/(252x⁶) + ...
    let inv_x = 1.0 / x;
    let inv_x2 = inv_x * inv_x;
    val += x.ln()
        - 0.5 * inv_x
        - inv_x2 * (1.0 / 12.0 - inv_x2 * (1.0 / 120.0 - inv_x2 * (1.0 / 252.0)));

    val
}

/// Trigamma function ψ'(x) = d²/dx² ln(Γ(x)).
///
/// Asymptotic expansion for x ≥ 8, with recurrence for smaller x.
fn trigamma(x: f64) -> f64 {
    if x <= 0.0 {
        return f64::NAN;
    }

    // Recurrence: ψ'(x) = ψ'(x+1) + 1/x²
    let mut val = 0.0;
    let mut x = x;
    while x < 8.0 {
        val += 1.0 / (x * x);
        x += 1.0;
    }

    // Asymptotic expansion for x ≥ 8
    // ψ'(x) ≈ 1/x + 1/(2x²) + 1/(6x³) - 1/(30x⁵) + 1/(42x⁷) - ...
    let inv_x = 1.0 / x;
    let inv_x2 = inv_x * inv_x;
    val += inv_x
        + 0.5 * inv_x2
        + inv_x2 * inv_x * (1.0 / 6.0 - inv_x2 * (1.0 / 30.0 - inv_x2 * (1.0 / 42.0)));

    val
}

/// Fits a Gamma distribution Gamma(α, β) via MLE.
///
/// # Algorithm
///
/// Rate parameter β̂ = α̂/x̄ (closed-form once α is known).
/// Shape parameter α̂ solved via Newton-Raphson on:
/// ln(α) - ψ(α) = ln(x̄) - (1/n)Σln(xᵢ)
///
/// Initial guess from method of moments: α₀ = x̄²/s².
///
/// Reference: Minka (2002), "Estimating a Gamma distribution"
///
/// # Returns
///
/// `None` if fewer than 2 data points, non-positive values, non-finite values,
/// or Newton-Raphson fails to converge.
///
/// # Examples
///
/// ```
/// use u_analytics::distribution::fit_gamma;
///
/// let data = [2.1, 3.5, 1.8, 4.2, 2.9, 3.1, 2.5, 3.8, 1.5, 2.7];
/// let fit = fit_gamma(&data).unwrap();
/// assert_eq!(fit.distribution, "Gamma");
/// assert!(fit.parameters[0].1 > 0.0); // α̂ > 0
/// assert!(fit.parameters[1].1 > 0.0); // β̂ > 0
/// ```
pub fn fit_gamma(data: &[f64]) -> Option<FitResult> {
    let n = data.len();
    if n < 2 || data.iter().any(|v| !v.is_finite()) {
        return None;
    }
    // Gamma requires x > 0
    if data.iter().any(|&v| v <= 0.0) {
        return None;
    }

    let nf = n as f64;
    let mean = stats::mean(data)?;
    if mean < 1e-300 {
        return None;
    }

    // Biased variance (denominator n)
    let var: f64 = data.iter().map(|&x| (x - mean).powi(2)).sum::<f64>() / nf;
    if var < 1e-300 {
        return None;
    }

    let sum_log: f64 = data.iter().map(|&x| x.ln()).sum();
    let mean_log = sum_log / nf;
    let target = mean.ln() - mean_log; // ln(x̄) - (1/n)Σln(xᵢ)

    // Method of moments initial guess
    let mut alpha = mean * mean / var;
    if alpha < 1e-10 {
        alpha = 0.5;
    }

    // Newton-Raphson: solve ln(α) - ψ(α) = target
    for _ in 0..200 {
        let f = alpha.ln() - digamma(alpha) - target;
        let f_prime = 1.0 / alpha - trigamma(alpha);
        if f_prime.abs() < 1e-30 {
            break;
        }

        let delta = f / f_prime;
        alpha -= delta;

        // Keep alpha positive
        if alpha <= 0.0 {
            alpha = 1e-10;
        }

        if delta.abs() < 1e-10 {
            break;
        }
    }

    if !alpha.is_finite() || alpha <= 0.0 {
        return None;
    }

    let beta = alpha / mean; // rate parameter

    // Log-likelihood: n·α·ln(β) - n·ln(Γ(α)) + (α-1)·Σln(xᵢ) - β·Σxᵢ
    let sum_x: f64 = data.iter().sum();
    let log_lik = nf * alpha * beta.ln() - nf * special::ln_gamma(alpha) + (alpha - 1.0) * sum_log
        - beta * sum_x;

    let k = 2;
    let aic = -2.0 * log_lik + 2.0 * k as f64;
    let bic = -2.0 * log_lik + k as f64 * nf.ln();

    Some(FitResult {
        distribution: "Gamma".to_string(),
        parameters: vec![("alpha".to_string(), alpha), ("beta".to_string(), beta)],
        log_likelihood: log_lik,
        aic,
        bic,
        n_params: k,
    })
}

/// Fits a LogNormal distribution via MLE.
///
/// # Estimators
///
/// If Y = ln(X), then Y ~ N(μ, σ²):
/// - μ̂ = (1/n) Σ ln(xᵢ)
/// - σ̂ = √((1/n) Σ (ln(xᵢ) - μ̂)²)
///
/// # Returns
///
/// `None` if fewer than 2 data points, non-positive values, or non-finite values.
///
/// # Examples
///
/// ```
/// use u_analytics::distribution::fit_lognormal;
///
/// let data = [1.2, 2.5, 1.8, 3.1, 0.9, 2.0, 1.5, 4.0];
/// let fit = fit_lognormal(&data).unwrap();
/// assert_eq!(fit.distribution, "LogNormal");
/// assert!(fit.parameters[0].1.is_finite()); // μ̂
/// assert!(fit.parameters[1].1 > 0.0); // σ̂ > 0
/// ```
pub fn fit_lognormal(data: &[f64]) -> Option<FitResult> {
    let n = data.len();
    if n < 2 || data.iter().any(|v| !v.is_finite()) {
        return None;
    }
    if data.iter().any(|&v| v <= 0.0) {
        return None;
    }

    let nf = n as f64;
    let log_data: Vec<f64> = data.iter().map(|&x| x.ln()).collect();

    let mu = log_data.iter().sum::<f64>() / nf;
    let sigma_sq = log_data.iter().map(|&y| (y - mu).powi(2)).sum::<f64>() / nf;
    if sigma_sq < 1e-300 {
        return None;
    }
    let sigma = sigma_sq.sqrt();

    // Log-likelihood: -n/2·ln(2π) - n·ln(σ) - (1/2σ²)Σ(ln(xᵢ)-μ)² - Σln(xᵢ)
    let sum_log: f64 = log_data.iter().sum();
    let log_lik =
        -nf / 2.0 * (2.0 * std::f64::consts::PI).ln() - nf * sigma.ln() - nf / 2.0 - sum_log;

    let k = 2;
    let aic = -2.0 * log_lik + 2.0 * k as f64;
    let bic = -2.0 * log_lik + k as f64 * nf.ln();

    Some(FitResult {
        distribution: "LogNormal".to_string(),
        parameters: vec![("mu".to_string(), mu), ("sigma".to_string(), sigma)],
        log_likelihood: log_lik,
        aic,
        bic,
        n_params: k,
    })
}

/// Fits a Poisson distribution via MLE.
///
/// # Estimators
///
/// λ̂ = x̄ (sample mean). Data are counts: non-negative integers
/// (represented as f64 for API consistency).
///
/// # Returns
///
/// `None` if fewer than 2 data points, a value that is not a non-negative
/// integer, all zeros, or non-finite values. (A non-integer has no Poisson
/// probability; `ln Γ(x + 1)` would extend the formula past where it means one.)
///
/// # Examples
///
/// ```
/// use u_analytics::distribution::fit_poisson;
///
/// let data = [2.0, 3.0, 1.0, 4.0, 2.0, 3.0, 1.0, 2.0];
/// let fit = fit_poisson(&data).unwrap();
/// assert_eq!(fit.distribution, "Poisson");
/// assert!((fit.parameters[0].1 - 2.25).abs() < 1e-10); // λ̂ = mean
/// ```
pub fn fit_poisson(data: &[f64]) -> Option<FitResult> {
    let n = data.len();
    if n < 2 || data.iter().any(|v| !v.is_finite()) {
        return None;
    }
    if data.iter().any(|&v| v < 0.0 || v.fract() != 0.0) {
        return None;
    }

    let nf = n as f64;
    let lambda = stats::mean(data)?;
    if lambda < 1e-300 {
        return None; // All zeros → degenerate
    }

    // Log-likelihood: Σ(xᵢ·ln(λ) - λ - ln(xᵢ!))
    // = n·x̄·ln(λ) - n·λ - Σln(xᵢ!)
    let sum_log_fact: f64 = data.iter().map(|&x| special::ln_gamma(x + 1.0)).sum();
    let log_lik = nf * lambda * lambda.ln() - nf * lambda - sum_log_fact;

    let k = 1;
    let aic = -2.0 * log_lik + 2.0 * k as f64;
    let bic = -2.0 * log_lik + k as f64 * nf.ln();

    Some(FitResult {
        distribution: "Poisson".to_string(),
        parameters: vec![("lambda".to_string(), lambda)],
        log_likelihood: log_lik,
        aic,
        bic,
        n_params: k,
    })
}

/// Fits a Beta distribution Beta(α, β) via MLE.
///
/// # Algorithm
///
/// Newton-Raphson on the system:
/// - ψ(α) - ψ(α + β) = (1/n) Σ ln(xᵢ)
/// - ψ(β) - ψ(α + β) = (1/n) Σ ln(1 - xᵢ)
///
/// Initial guess from method of moments:
/// - α₀ = x̄ · ((x̄(1-x̄)/s²) - 1)
/// - β₀ = (1-x̄) · ((x̄(1-x̄)/s²) - 1)
///
/// Data must be in (0, 1).
///
/// # Returns
///
/// `None` if fewer than 2 data points, values outside (0,1), or convergence failure.
///
/// # Examples
///
/// ```
/// use u_analytics::distribution::fit_beta;
///
/// let data = [0.2, 0.35, 0.5, 0.15, 0.45, 0.3, 0.6, 0.25, 0.4, 0.55];
/// let fit = fit_beta(&data).unwrap();
/// assert_eq!(fit.distribution, "Beta");
/// assert!(fit.parameters[0].1 > 0.0); // α̂ > 0
/// assert!(fit.parameters[1].1 > 0.0); // β̂ > 0
/// ```
pub fn fit_beta(data: &[f64]) -> Option<FitResult> {
    let n = data.len();
    if n < 2 || data.iter().any(|v| !v.is_finite()) {
        return None;
    }
    // Beta requires 0 < x < 1
    if data.iter().any(|&v| v <= 0.0 || v >= 1.0) {
        return None;
    }

    let nf = n as f64;
    let mean = stats::mean(data)?;
    let var = data.iter().map(|&x| (x - mean).powi(2)).sum::<f64>() / nf;
    if var < 1e-300 {
        return None;
    }

    let sum_log_x: f64 = data.iter().map(|&x| x.ln()).sum();
    let sum_log_1mx: f64 = data.iter().map(|&x| (1.0 - x).ln()).sum();
    let mean_log_x = sum_log_x / nf;
    let mean_log_1mx = sum_log_1mx / nf;

    // Method of moments initial guess
    let factor = mean * (1.0 - mean) / var - 1.0;
    if factor <= 0.0 {
        return None;
    }
    let mut alpha = mean * factor;
    let mut beta = (1.0 - mean) * factor;

    if alpha < 0.01 {
        alpha = 0.5;
    }
    if beta < 0.01 {
        beta = 0.5;
    }

    // Newton-Raphson (joint update)
    for _ in 0..200 {
        let psi_ab = digamma(alpha + beta);
        let tri_ab = trigamma(alpha + beta);

        let f1 = digamma(alpha) - psi_ab - mean_log_x;
        let f2 = digamma(beta) - psi_ab - mean_log_1mx;

        let j11 = trigamma(alpha) - tri_ab;
        let j12 = -tri_ab;
        let j21 = -tri_ab;
        let j22 = trigamma(beta) - tri_ab;

        let det = j11 * j22 - j12 * j21;
        if det.abs() < 1e-30 {
            break;
        }

        let da = (j22 * f1 - j12 * f2) / det;
        let db = (j11 * f2 - j21 * f1) / det;

        alpha -= da;
        beta -= db;

        if alpha <= 0.0 {
            alpha = 1e-10;
        }
        if beta <= 0.0 {
            beta = 1e-10;
        }

        if da.abs() < 1e-10 && db.abs() < 1e-10 {
            break;
        }
    }

    if !alpha.is_finite() || !beta.is_finite() || alpha <= 0.0 || beta <= 0.0 {
        return None;
    }

    // Log-likelihood: n·[ln(Γ(α+β)) - ln(Γ(α)) - ln(Γ(β))]
    //                + (α-1)·Σln(xᵢ) + (β-1)·Σln(1-xᵢ)
    let log_lik = nf
        * (special::ln_gamma(alpha + beta) - special::ln_gamma(alpha) - special::ln_gamma(beta))
        + (alpha - 1.0) * sum_log_x
        + (beta - 1.0) * sum_log_1mx;

    let k = 2;
    let aic = -2.0 * log_lik + 2.0 * k as f64;
    let bic = -2.0 * log_lik + k as f64 * nf.ln();

    Some(FitResult {
        distribution: "Beta".to_string(),
        parameters: vec![("alpha".to_string(), alpha), ("beta".to_string(), beta)],
        log_likelihood: log_lik,
        aic,
        bic,
        n_params: k,
    })
}

/// Fits a two-parameter Weibull distribution by maximum likelihood
/// ([`crate::weibull::weibull_mle`]), as a [`FitResult`] for comparison with
/// the other families.
///
/// # Returns
///
/// `None` if fewer than 2 data points, a value that is not strictly positive
/// and finite, or the fit does not converge.
///
/// # Examples
///
/// ```
/// use u_analytics::distribution::fit_weibull;
///
/// let data = [12.0, 25.0, 31.0, 44.0, 52.0, 60.0, 71.0, 85.0];
/// let fit = fit_weibull(&data).unwrap();
/// assert_eq!(fit.distribution, "Weibull");
/// assert_eq!(fit.parameters[0].0, "shape");
/// ```
pub fn fit_weibull(data: &[f64]) -> Option<FitResult> {
    if data.len() < 2 || data.iter().any(|&v| !(v.is_finite() && v > 0.0)) {
        return None;
    }
    let r = crate::weibull::weibull_mle(data)?;
    let nf = data.len() as f64;
    let k = 2;
    Some(FitResult {
        distribution: "Weibull".to_string(),
        parameters: vec![
            ("shape".to_string(), r.shape),
            ("scale".to_string(), r.scale),
        ],
        log_likelihood: r.log_likelihood,
        aic: -2.0 * r.log_likelihood + 2.0 * k as f64,
        bic: -2.0 * r.log_likelihood + k as f64 * nf.ln(),
        n_params: k,
    })
}

/// Fits the continuous families that apply to the data and returns them sorted
/// by AIC (best first).
///
/// Tries Normal, Exponential, Gamma, LogNormal, Weibull and Beta; each is kept
/// only if it fits (the positive families need data > 0, Beta needs data in
/// (0, 1)).
///
/// Poisson is not a candidate: its likelihood is a probability mass, the
/// others' a density, and AIC compares likelihoods only on the same footing --
/// a density's value changes with the unit of measurement, a mass does not.
/// Fit count data with [`fit_poisson`] directly.
///
/// # Examples
///
/// ```
/// use u_analytics::distribution::fit_best;
///
/// let data = [2.1, 3.5, 1.8, 4.2, 2.9, 3.1, 2.5, 3.8, 1.5, 2.7];
/// let fits = fit_best(&data);
/// assert!(!fits.is_empty());
/// // Results sorted by AIC (lowest = best fit)
/// if fits.len() >= 2 {
///     assert!(fits[0].aic <= fits[1].aic);
/// }
/// ```
pub fn fit_best(data: &[f64]) -> Vec<FitResult> {
    let mut results = Vec::new();

    if let Some(r) = fit_normal(data) {
        results.push(r);
    }
    if let Some(r) = fit_exponential(data) {
        results.push(r);
    }
    if let Some(r) = fit_gamma(data) {
        results.push(r);
    }
    if let Some(r) = fit_lognormal(data) {
        results.push(r);
    }
    if let Some(r) = fit_weibull(data) {
        results.push(r);
    }
    if let Some(r) = fit_beta(data) {
        results.push(r);
    }

    results.sort_by(|a, b| {
        a.aic
            .partial_cmp(&b.aic)
            .unwrap_or(std::cmp::Ordering::Equal)
    });
    results
}

// ---------------------------------------------------------------------------
// Goodness of fit of a fitted family
// ---------------------------------------------------------------------------

/// A family [`fit_best`] reports, with its fitted parameters.
#[derive(Debug, Clone)]
enum Fitted {
    Normal(Normal),
    LogNormal(LogNormal),
    Weibull(Weibull),
    Gamma(GammaDistribution),
    Exponential(Exponential),
    Beta(BetaDistribution),
}

impl Fitted {
    /// The distribution a [`FitResult`] names, or `None` for a family outside
    /// [`fit_best`]'s (Poisson) or parameters it cannot be built from.
    fn of(fit: &FitResult) -> Option<Self> {
        let param = |name: &str| {
            fit.parameters
                .iter()
                .find(|(n, _)| n == name)
                .map(|&(_, v)| v)
        };
        Some(match fit.distribution.as_str() {
            "Normal" => Self::Normal(Normal::new(param("mu")?, param("sigma")?).ok()?),
            "LogNormal" => Self::LogNormal(LogNormal::new(param("mu")?, param("sigma")?).ok()?),
            "Weibull" => Self::Weibull(Weibull::new(param("shape")?, param("scale")?).ok()?),
            // `fit_gamma` reports the rate as `beta`.
            "Gamma" => Self::Gamma(GammaDistribution::new(param("alpha")?, param("beta")?).ok()?),
            "Exponential" => Self::Exponential(Exponential::new(param("lambda")?).ok()?),
            "Beta" => Self::Beta(BetaDistribution::new(param("alpha")?, param("beta")?).ok()?),
            _ => return None,
        })
    }

    /// Refits the same family to `data`.
    fn refit(&self, data: &[f64]) -> Option<FitResult> {
        match self {
            Self::Normal(_) => fit_normal(data),
            Self::LogNormal(_) => fit_lognormal(data),
            Self::Weibull(_) => fit_weibull(data),
            Self::Gamma(_) => fit_gamma(data),
            Self::Exponential(_) => fit_exponential(data),
            Self::Beta(_) => fit_beta(data),
        }
    }

    fn quantile(&self, p: f64) -> Option<f64> {
        match self {
            Self::Normal(d) => d.quantile(p),
            Self::LogNormal(d) => d.quantile(p),
            Self::Weibull(d) => d.quantile(p),
            Self::Gamma(d) => d.quantile(p),
            Self::Exponential(d) => d.quantile(p),
            Self::Beta(d) => d.quantile(p),
        }
    }

    /// `(ln F(x), ln (1 − F(x)))`, each tail evaluated directly so that a value
    /// far out contributes its true log-probability rather than `ln 0`.
    fn ln_tails(&self, x: f64) -> (f64, f64) {
        let (cdf, sf) = match self {
            Self::Normal(d) => {
                let z = (x - d.mu()) / d.sigma();
                (
                    special::standard_normal_cdf(z),
                    special::standard_normal_sf(z),
                )
            }
            Self::LogNormal(d) => {
                if x <= 0.0 {
                    (0.0, 1.0)
                } else {
                    let z = (x.ln() - d.mu()) / d.sigma();
                    (
                        special::standard_normal_cdf(z),
                        special::standard_normal_sf(z),
                    )
                }
            }
            Self::Weibull(d) => {
                if x <= 0.0 {
                    (0.0, 1.0)
                } else {
                    let h = (x / d.scale()).powf(d.shape());
                    (-(-h).exp_m1(), (-h).exp())
                }
            }
            Self::Gamma(d) => {
                if x <= 0.0 {
                    (0.0, 1.0)
                } else {
                    let t = d.rate() * x;
                    (
                        special::regularized_lower_gamma(d.shape(), t),
                        special::regularized_upper_gamma(d.shape(), t),
                    )
                }
            }
            Self::Exponential(d) => (d.cdf(x), d.survival(x)),
            Self::Beta(d) => {
                if x <= 0.0 {
                    (0.0, 1.0)
                } else if x >= 1.0 {
                    (1.0, 0.0)
                } else {
                    (
                        special::regularized_incomplete_beta(x, d.alpha, d.beta),
                        special::regularized_incomplete_beta(1.0 - x, d.beta, d.alpha),
                    )
                }
            }
        };
        (
            cdf.max(f64::MIN_POSITIVE).ln(),
            sf.max(f64::MIN_POSITIVE).ln(),
        )
    }

    fn sample_n(&self, n: usize, seed: u64) -> Vec<f64> {
        match self {
            Self::Normal(d) => d.sample_n(n, seed),
            Self::LogNormal(d) => d.sample_n(n, seed),
            Self::Weibull(d) => d.sample_n(n, seed),
            Self::Gamma(d) => d.sample_n(n, seed),
            Self::Exponential(d) => d.sample_n(n, seed),
            Self::Beta(d) => d.sample_n(n, seed),
        }
    }

    /// A² of `data` against this distribution.
    fn anderson_darling(&self, data: &[f64]) -> f64 {
        let mut x = data.to_vec();
        x.sort_by(|a, b| a.partial_cmp(b).expect("finite values"));
        let n = x.len();
        let tails: Vec<(f64, f64)> = x.iter().map(|&v| self.ln_tails(v)).collect();
        let s: f64 = (0..n)
            .map(|i| (2 * i + 1) as f64 * (tails[i].0 + tails[n - 1 - i].1))
            .sum();
        -(n as f64) - s / n as f64
    }
}

/// The 0.135 %, 50 % and 99.865 % quantiles (or any others) of the family a
/// [`FitResult`] names — for the capability methods that read a fitted model.
pub(crate) fn fitted_quantile(fit: &FitResult, p: f64) -> Option<f64> {
    Fitted::of(fit)?.quantile(p)
}

/// Anderson–Darling statistic A² of `data` against the distribution `fit`
/// describes (its fitted parameters, not re-estimated).
///
/// A² = −n − (1/n) Σᵢ (2i − 1)·[ln F(x₍ᵢ₎) + ln(1 − F(x₍ₙ₊₁₋ᵢ₎))]. It weights
/// the tails, where a wrong family misplaces quantiles most. Smaller is closer;
/// among families fitted to the same data it ranks them by fit rather than by
/// likelihood. To judge whether the family fits at all, use
/// [`goodness_of_fit`] — the null distribution of A² depends on the family and
/// on the parameters having been estimated from the same data.
///
/// # Returns
///
/// `None` if `data` is empty or has a non-finite value, or `fit` is not one of
/// the families [`fit_best`] reports.
///
/// # References
///
/// Anderson, T.W. & Darling, D.A. (1954). "A test of goodness of fit."
/// *JASA*, 49(268), 765–769.
///
/// # Examples
///
/// ```
/// use u_analytics::distribution::{anderson_darling_fit, fit_normal};
///
/// let data = [4.1, 5.3, 4.8, 5.9, 5.0, 4.4, 5.6, 4.9, 5.2, 4.7];
/// let fit = fit_normal(&data).unwrap();
/// let a2 = anderson_darling_fit(&data, &fit).unwrap();
/// assert!(a2 > 0.0 && a2 < 1.0);
/// ```
pub fn anderson_darling_fit(data: &[f64], fit: &FitResult) -> Option<f64> {
    if data.is_empty() || data.iter().any(|v| !v.is_finite()) {
        return None;
    }
    Some(Fitted::of(fit)?.anderson_darling(data))
}

/// Whether a fitted family describes the data: the Anderson–Darling statistic
/// against the fitted distribution and its p-value.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct GoodnessOfFit {
    /// A² of the data against the fitted distribution ([`anderson_darling_fit`]).
    pub statistic: f64,
    /// Probability, if the data came from the family, of an A² at least this
    /// large after fitting the family to a sample of the same size — the
    /// parametric-bootstrap estimate `(1 + #{A²_b ≥ A²}) / (B + 1)`. Small
    /// values say the family does not describe the data.
    pub p_value: f64,
    /// Bootstrap replicates `B` behind `p_value`; the smallest p it can report
    /// is `1 / (B + 1)`.
    pub replicates: usize,
}

/// Replicates [`goodness_of_fit`] uses when a caller has no reason to choose:
/// with 199 the p-value is resolved to 0.005, and a 0.05 or 0.01 cut falls on
/// a whole number of replicates ((B + 1)·α integer — Davison & Hinkley §4.2).
pub const DEFAULT_GOF_REPLICATES: usize = 199;

/// Tests whether the family `fit` describes `data`, with the family's
/// parameters estimated from that same data.
///
/// The statistic is A² against the fitted distribution
/// ([`anderson_darling_fit`]). Its null distribution is not the textbook one
/// when parameters are estimated — fitting pulls the model toward the sample —
/// and it depends on the family (and, for Gamma and Beta, on the shape). The
/// p-value is therefore the **parametric bootstrap**: draw `replicates` samples
/// of size n from the fitted distribution, refit the same family to each, and
/// count how often their A² reaches the observed one. This applies the same
/// way to every family, where tables exist only for some.
///
/// Deterministic: the same `(data, fit, replicates, seed)` gives the same
/// p-value.
///
/// # Complexity
///
/// O(B · (n log n + fit)), where `fit` is the family's maximum-likelihood
/// iteration (closed form for Normal, LogNormal and Exponential).
///
/// # Returns
///
/// `None` if `replicates` is 0, `data` has fewer than 2 values or a non-finite
/// one, or `fit` is not one of the families [`fit_best`] reports.
///
/// # References
///
/// - Stute, W., González Manteiga, W. & Presedo Quindimil, M. (1993).
///   "Bootstrap based goodness-of-fit-tests." *Metrika*, 40, 243–256.
/// - Davison, A.C. & Hinkley, D.V. (1997). *Bootstrap Methods and their
///   Application*, §4.2 (the `(1 + count) / (B + 1)` estimate).
/// - Stephens, M.A. (1974). "EDF statistics for goodness of fit and some
///   comparisons." *JASA*, 69(347), 730–737 (A² with estimated parameters).
///
/// # Examples
///
/// ```
/// use u_analytics::distribution::{fit_best, goodness_of_fit, DEFAULT_GOF_REPLICATES};
///
/// // Uniform data: no listed family describes it.
/// let data: Vec<f64> = (1..=500).map(f64::from).collect();
/// let best = &fit_best(&data)[0];
/// let gof = goodness_of_fit(&data, best, DEFAULT_GOF_REPLICATES, 0).unwrap();
/// assert!(gof.p_value < 0.01);
/// ```
pub fn goodness_of_fit(
    data: &[f64],
    fit: &FitResult,
    replicates: usize,
    seed: u64,
) -> Option<GoodnessOfFit> {
    if replicates == 0 || data.len() < 2 || data.iter().any(|v| !v.is_finite()) {
        return None;
    }
    let model = Fitted::of(fit)?;
    let statistic = model.anderson_darling(data);
    let n = data.len();
    let at_least = (0..replicates as u64)
        .filter(|&b| {
            let sample = model.sample_n(n, seed.wrapping_add(b));
            // A replicate the family cannot refit (a draw that degenerates)
            // counts as at least as extreme: the estimate stays conservative.
            match model.refit(&sample).and_then(|f| Fitted::of(&f)) {
                Some(refitted) => refitted.anderson_darling(&sample) >= statistic,
                None => true,
            }
        })
        .count();
    Some(GoodnessOfFit {
        statistic,
        p_value: (1 + at_least) as f64 / (replicates + 1) as f64,
        replicates,
    })
}

#[cfg(test)]
mod tests {
    use super::*;

    // -----------------------------------------------------------------------
    // ECDF
    // -----------------------------------------------------------------------

    #[test]
    fn ecdf_basic() {
        let (vals, probs) = ecdf(&[3.0, 1.0, 2.0]).expect("should compute");
        assert_eq!(vals, vec![1.0, 2.0, 3.0]);
        assert!((probs[0] - 1.0 / 3.0).abs() < 1e-10);
        assert!((probs[1] - 2.0 / 3.0).abs() < 1e-10);
        assert!((probs[2] - 1.0).abs() < 1e-10);
    }

    #[test]
    fn ecdf_with_ties() {
        let (vals, probs) = ecdf(&[1.0, 2.0, 2.0, 3.0]).expect("should compute");
        assert_eq!(vals, vec![1.0, 2.0, 3.0]);
        assert!((probs[0] - 0.25).abs() < 1e-10);
        assert!((probs[1] - 0.75).abs() < 1e-10); // 3 out of 4 ≤ 2.0
        assert!((probs[2] - 1.0).abs() < 1e-10);
    }

    #[test]
    fn ecdf_single() {
        let (vals, probs) = ecdf(&[5.0]).expect("should compute");
        assert_eq!(vals, vec![5.0]);
        assert!((probs[0] - 1.0).abs() < 1e-10);
    }

    #[test]
    fn ecdf_empty() {
        assert!(ecdf(&[]).is_none());
    }

    #[test]
    fn ecdf_nan() {
        assert!(ecdf(&[1.0, f64::NAN, 3.0]).is_none());
    }

    // -----------------------------------------------------------------------
    // Histogram bins
    // -----------------------------------------------------------------------

    #[test]
    fn hist_sturges_basic() {
        let data: Vec<f64> = (0..100).map(|i| i as f64).collect();
        let r = histogram_bins(&data, BinMethod::Sturges).expect("should compute");
        assert!(r.n_bins >= 5);
        assert_eq!(r.edges.len(), r.n_bins + 1);
        assert_eq!(r.counts.len(), r.n_bins);
        let total: usize = r.counts.iter().sum();
        assert_eq!(total, 100);
    }

    #[test]
    fn hist_scott() {
        let data: Vec<f64> = (0..100).map(|i| i as f64).collect();
        let r = histogram_bins(&data, BinMethod::Scott).expect("should compute");
        assert!(r.n_bins >= 3);
        let total: usize = r.counts.iter().sum();
        assert_eq!(total, 100);
    }

    #[test]
    fn hist_fd() {
        let data: Vec<f64> = (0..100).map(|i| i as f64).collect();
        let r = histogram_bins(&data, BinMethod::FreedmanDiaconis).expect("should compute");
        assert!(r.n_bins >= 3);
        let total: usize = r.counts.iter().sum();
        assert_eq!(total, 100);
    }

    #[test]
    fn hist_edge_cases() {
        assert!(histogram_bins(&[1.0], BinMethod::Sturges).is_none()); // < 2
        assert!(histogram_bins(&[5.0, 5.0, 5.0], BinMethod::Sturges).is_none());
        // zero range
    }

    #[test]
    fn hist_fixed_exact_count() {
        let data: Vec<f64> = (0..100).map(|i| i as f64).collect();
        let r = histogram_bins(&data, BinMethod::Fixed(5)).expect("should compute");
        assert_eq!(r.n_bins, 5);
        assert_eq!(r.edges.len(), 6);
        assert_eq!(r.counts.len(), 5);
        let total: usize = r.counts.iter().sum();
        assert_eq!(total, 100);
    }

    #[test]
    fn hist_fixed_single_bin() {
        let data = [1.0, 2.0, 3.0];
        let r = histogram_bins(&data, BinMethod::Fixed(1)).expect("should compute");
        assert_eq!(r.n_bins, 1);
        assert_eq!(r.counts, vec![3]);
    }

    #[test]
    fn hist_fixed_zero_is_invalid() {
        let data = [1.0, 2.0, 3.0];
        assert!(histogram_bins(&data, BinMethod::Fixed(0)).is_none());
    }

    #[test]
    fn hist_fixed_shares_input_guards() {
        assert!(histogram_bins(&[1.0], BinMethod::Fixed(5)).is_none()); // < 2 points
        assert!(histogram_bins(&[5.0, 5.0, 5.0], BinMethod::Fixed(5)).is_none());
        // zero range
    }

    // -----------------------------------------------------------------------
    // QQ-plot
    // -----------------------------------------------------------------------

    #[test]
    fn qq_basic() {
        let data = [-1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5];
        let (theo, samp) = qq_plot_normal(&data).expect("should compute");
        assert_eq!(theo.len(), 7);
        assert_eq!(samp.len(), 7);
        // Sorted sample
        assert!((samp[0] + 1.5).abs() < 1e-10);
        assert!((samp[6] - 1.5).abs() < 1e-10);
        // Theoretical: median should be near 0
        assert!(theo[3].abs() < 0.1);
    }

    #[test]
    fn qq_edge_cases() {
        assert!(qq_plot_normal(&[1.0, 2.0]).is_none()); // < 3
    }

    // -----------------------------------------------------------------------
    // KS test
    // -----------------------------------------------------------------------

    #[test]
    fn ks_normal_data() {
        // Approximately normal data
        let data = [-2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5];
        let (d, p) = ks_test_normal(&data).expect("should compute");
        assert!(d > 0.0 && d < 1.0, "D = {d}");
        assert!(p > 0.05, "p = {p}");
    }

    #[test]
    fn ks_non_normal_data() {
        // Clearly non-normal (uniform 0-100)
        let data: Vec<f64> = (0..50).map(|i| i as f64 * 2.0).collect();
        let (d, _p) = ks_test_normal(&data).expect("should compute");
        assert!(d > 0.0);
    }

    #[test]
    fn ks_edge_cases() {
        assert!(ks_test_normal(&[1.0, 2.0, 3.0, 4.0]).is_none()); // < 5
        assert!(ks_test_normal(&[5.0, 5.0, 5.0, 5.0, 5.0]).is_none()); // zero var
    }

    // -----------------------------------------------------------------------
    // KDE
    // -----------------------------------------------------------------------

    #[test]
    fn kde_silverman_basic() {
        let data: Vec<f64> = (0..100).map(|i| i as f64).collect();
        let r = kde(&data, BandwidthMethod::Silverman, 256).expect("should compute");
        assert_eq!(r.x.len(), 256);
        assert_eq!(r.density.len(), 256);
        assert!(r.bandwidth > 0.0);
        // All densities non-negative
        assert!(r.density.iter().all(|&d| d >= 0.0));
        // Integral ≈ 1
        let dx = r.x[1] - r.x[0];
        let integral: f64 = r.density.iter().sum::<f64>() * dx;
        assert!(
            (integral - 1.0).abs() < 0.05,
            "integral = {integral}, expected ≈ 1.0"
        );
    }

    #[test]
    fn kde_scott_basic() {
        let data: Vec<f64> = (0..100).map(|i| i as f64).collect();
        let r = kde(&data, BandwidthMethod::Scott, 256).expect("should compute");
        assert!(r.bandwidth > 0.0);
        let dx = r.x[1] - r.x[0];
        let integral: f64 = r.density.iter().sum::<f64>() * dx;
        assert!(
            (integral - 1.0).abs() < 0.05,
            "integral = {integral}, expected ≈ 1.0"
        );
    }

    #[test]
    fn kde_manual_bandwidth() {
        let data = [1.0, 2.0, 3.0, 4.0, 5.0];
        let r = kde(&data, BandwidthMethod::Manual(0.5), 128).expect("should compute");
        assert!((r.bandwidth - 0.5).abs() < 1e-10);
        assert_eq!(r.x.len(), 128);
    }

    #[test]
    fn kde_edge_cases() {
        // Too few points
        assert!(kde(&[1.0], BandwidthMethod::Silverman, 256).is_none());
        // Empty
        assert!(kde(&[], BandwidthMethod::Silverman, 256).is_none());
        // NaN
        assert!(kde(&[1.0, f64::NAN], BandwidthMethod::Silverman, 256).is_none());
        // Constant data (zero variance)
        assert!(kde(&[5.0, 5.0, 5.0], BandwidthMethod::Silverman, 256).is_none());
        // Too few grid points
        assert!(kde(&[1.0, 2.0, 3.0], BandwidthMethod::Manual(1.0), 1).is_none());
        // Invalid manual bandwidth
        assert!(kde(&[1.0, 2.0, 3.0], BandwidthMethod::Manual(0.0), 256).is_none());
        assert!(kde(&[1.0, 2.0, 3.0], BandwidthMethod::Manual(-1.0), 256).is_none());
    }

    #[test]
    fn kde_bimodal() {
        // Two clusters: around 0 and around 10
        let mut data = Vec::new();
        for i in 0..50 {
            data.push(i as f64 * 0.1); // 0.0 to 4.9
        }
        for i in 0..50 {
            data.push(10.0 + i as f64 * 0.1); // 10.0 to 14.9
        }
        let r = kde(&data, BandwidthMethod::Silverman, 512).expect("should compute");
        // Should have higher density near the two clusters
        // Find density at x ≈ 2.5 (in cluster 1) and x ≈ 7.5 (between clusters)
        let idx_cluster = r.x.iter().position(|&x| x >= 2.5).expect("grid point");
        let idx_valley = r.x.iter().position(|&x| x >= 7.5).expect("grid point");
        assert!(
            r.density[idx_cluster] > r.density[idx_valley],
            "cluster density ({}) should exceed valley density ({})",
            r.density[idx_cluster],
            r.density[idx_valley]
        );
    }

    #[test]
    fn kde_evaluate_single() {
        let data = [0.0, 1.0, 2.0, 3.0, 4.0];
        let h = 1.0;
        let d = kde_evaluate(&data, h, 2.0).expect("should compute");
        assert!(d > 0.0);
        // At the center of the data, density should be highest
        let d_edge = kde_evaluate(&data, h, -5.0).expect("should compute");
        assert!(d > d_edge, "center density ({d}) > edge density ({d_edge})");
    }

    #[test]
    fn kde_evaluate_edge_cases() {
        assert!(kde_evaluate(&[], 1.0, 0.0).is_none());
        assert!(kde_evaluate(&[1.0], 0.0, 0.0).is_none());
        assert!(kde_evaluate(&[1.0], -1.0, 0.0).is_none());
        assert!(kde_evaluate(&[1.0], 1.0, f64::NAN).is_none());
        assert!(kde_evaluate(&[f64::NAN], 1.0, 0.0).is_none());
    }

    #[test]
    fn kde_bandwidth_methods() {
        let data: Vec<f64> = (0..50).map(|i| i as f64).collect();
        let h_silv = kde_bandwidth(&data, BandwidthMethod::Silverman).expect("silverman");
        let h_scott = kde_bandwidth(&data, BandwidthMethod::Scott).expect("scott");
        let h_manual = kde_bandwidth(&data, BandwidthMethod::Manual(2.5)).expect("manual");
        assert!(h_silv > 0.0);
        assert!(h_scott > 0.0);
        assert!((h_manual - 2.5).abs() < 1e-10);
        // Scott produces slightly larger bandwidth than Silverman for uniform data
        // (1.06σ vs 0.9*min(σ, IQR/1.34))
        assert!(h_scott > h_silv);
    }

    #[test]
    fn kde_bandwidth_edge_cases() {
        assert!(kde_bandwidth(&[1.0], BandwidthMethod::Silverman).is_none());
        assert!(kde_bandwidth(&[5.0, 5.0], BandwidthMethod::Silverman).is_none());
        assert!(kde_bandwidth(&[1.0, 2.0], BandwidthMethod::Manual(0.0)).is_none());
    }

    #[test]
    fn kde_grid_extends_beyond_data() {
        let data = [10.0, 20.0, 30.0];
        let r = kde(&data, BandwidthMethod::Manual(2.0), 64).expect("should compute");
        // Grid should extend 3*h = 6 beyond data range [10, 30]
        assert!(r.x[0] < 10.0 - 2.0, "grid starts at {}", r.x[0]);
        assert!(
            *r.x.last().expect("non-empty") > 30.0 + 2.0,
            "grid ends at {}",
            r.x.last().expect("non-empty")
        );
    }

    // -----------------------------------------------------------------------
    // Distribution Fitting
    // -----------------------------------------------------------------------

    #[test]
    fn fit_normal_basic() {
        let data = [2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0];
        let r = fit_normal(&data).expect("should compute");
        assert_eq!(r.distribution, "Normal");
        assert_eq!(r.n_params, 2);
        assert!((r.parameters[0].1 - 5.0).abs() < 1e-10, "μ̂ = 5.0");
        // Biased MLE σ̂ = √(Σ(xᵢ-x̄)²/n) = √(32/8) = 2.0
        // (xᵢ-5)² = 9,1,1,1,0,0,4,16 → sum = 32
        let expected_sigma = (32.0_f64 / 8.0).sqrt();
        assert!(
            (r.parameters[1].1 - expected_sigma).abs() < 1e-10,
            "σ̂ = {}, expected {}",
            r.parameters[1].1,
            expected_sigma
        );
        assert!(r.log_likelihood.is_finite());
        assert!(r.aic.is_finite());
        assert!(r.bic.is_finite());
    }

    #[test]
    fn fit_normal_edge_cases() {
        assert!(fit_normal(&[1.0]).is_none()); // < 2
        assert!(fit_normal(&[]).is_none());
        assert!(fit_normal(&[1.0, f64::NAN]).is_none());
        assert!(fit_normal(&[5.0, 5.0, 5.0]).is_none()); // zero variance
    }

    #[test]
    fn fit_exponential_basic() {
        let data = [0.5, 1.2, 0.8, 2.1, 0.3, 1.5, 0.9, 1.8];
        let r = fit_exponential(&data).expect("should compute");
        assert_eq!(r.distribution, "Exponential");
        assert_eq!(r.n_params, 1);
        let mean: f64 = data.iter().sum::<f64>() / data.len() as f64;
        let expected_lambda = 1.0 / mean;
        assert!(
            (r.parameters[0].1 - expected_lambda).abs() < 1e-10,
            "λ̂ = {}, expected {}",
            r.parameters[0].1,
            expected_lambda
        );
    }

    #[test]
    fn fit_exponential_edge_cases() {
        assert!(fit_exponential(&[1.0]).is_none()); // < 2
        assert!(fit_exponential(&[1.0, -1.0]).is_none()); // negative
        assert!(fit_exponential(&[0.0, 1.0]).is_none()); // zero value
        assert!(fit_exponential(&[1.0, f64::NAN]).is_none());
    }

    #[test]
    fn fit_gamma_basic() {
        // Generate Gamma-like data (known α ≈ 4, β ≈ 2, mean = α/β = 2)
        let data = [
            1.5, 2.1, 1.8, 2.5, 1.9, 2.3, 2.0, 1.7, 2.4, 2.2, 1.6, 2.6, 1.4, 2.8, 2.1, 1.9, 2.3,
            2.0, 1.8, 2.5,
        ];
        let r = fit_gamma(&data).expect("should compute");
        assert_eq!(r.distribution, "Gamma");
        assert_eq!(r.n_params, 2);
        assert!(r.parameters[0].1 > 0.0, "α̂ > 0");
        assert!(r.parameters[1].1 > 0.0, "β̂ > 0");
        assert!(r.log_likelihood.is_finite());
    }

    #[test]
    fn fit_gamma_edge_cases() {
        assert!(fit_gamma(&[1.0]).is_none()); // < 2
        assert!(fit_gamma(&[1.0, -1.0]).is_none()); // negative
        assert!(fit_gamma(&[0.0, 1.0]).is_none()); // zero
        assert!(fit_gamma(&[1.0, f64::NAN]).is_none());
        assert!(fit_gamma(&[5.0, 5.0, 5.0]).is_none()); // zero variance
    }

    #[test]
    fn fit_gamma_mean_rate_relationship() {
        let data = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0];
        let r = fit_gamma(&data).expect("should compute");
        let alpha = r.parameters[0].1;
        let beta = r.parameters[1].1;
        let mean: f64 = data.iter().sum::<f64>() / data.len() as f64;
        // β̂ = α̂/x̄, so α̂/β̂ ≈ x̄
        assert!(
            (alpha / beta - mean).abs() < 1e-6,
            "α/β = {}, mean = {}",
            alpha / beta,
            mean
        );
    }

    #[test]
    fn fit_best_sorted_by_aic() {
        let data = [2.1, 3.5, 1.8, 4.2, 2.9, 3.1, 2.5, 3.8, 1.5, 2.7];
        let fits = fit_best(&data);
        assert!(!fits.is_empty());
        // Should include Normal and Gamma (and Exponential since all > 0)
        assert!(fits.len() >= 2, "got {} fits", fits.len());
        // Verify AIC ordering
        for w in fits.windows(2) {
            assert!(
                w[0].aic <= w[1].aic,
                "AIC not sorted: {} > {}",
                w[0].aic,
                w[1].aic
            );
        }
    }

    #[test]
    fn fit_best_mixed_sign_data() {
        // Negative values → only Normal should succeed
        let data = [-2.0, -1.0, 0.5, 1.0, 2.0, 3.0, 1.5, -0.5];
        let fits = fit_best(&data);
        assert_eq!(fits.len(), 1, "only Normal should fit");
        assert_eq!(fits[0].distribution, "Normal");
    }

    #[test]
    fn digamma_known_values() {
        // ψ(1) = -γ ≈ -0.5772156649
        assert!((digamma(1.0) + 0.5772156649).abs() < 1e-8);
        // ψ(2) = 1 - γ ≈ 0.4227843351
        assert!((digamma(2.0) - 0.4227843351).abs() < 1e-8);
        // ψ(0.5) = -γ - 2·ln(2) ≈ -1.9635100260
        assert!((digamma(0.5) + 1.9635100260).abs() < 1e-7);
    }

    #[test]
    fn trigamma_known_values() {
        // ψ'(1) = π²/6 ≈ 1.6449340668
        assert!((trigamma(1.0) - std::f64::consts::PI.powi(2) / 6.0).abs() < 1e-7);
        // ψ'(2) = π²/6 - 1 ≈ 0.6449340668
        assert!((trigamma(2.0) - (std::f64::consts::PI.powi(2) / 6.0 - 1.0)).abs() < 1e-7);
    }

    // -----------------------------------------------------------------------
    // LogNormal fitting
    // -----------------------------------------------------------------------

    #[test]
    fn fit_lognormal_basic() {
        let data = [1.2, 2.5, 1.8, 3.1, 0.9, 2.0, 1.5, 4.0];
        let r = fit_lognormal(&data).expect("should compute");
        assert_eq!(r.distribution, "LogNormal");
        assert_eq!(r.n_params, 2);
        assert!(r.parameters[1].1 > 0.0); // σ̂ > 0
        assert!(r.log_likelihood.is_finite());
    }

    #[test]
    fn fit_lognormal_params_match_log_transform() {
        let data = [1.0, 2.0, 3.0, 4.0, 5.0];
        let r = fit_lognormal(&data).expect("should compute");
        let mu = r.parameters[0].1;
        let sigma = r.parameters[1].1;
        // μ̂ should be mean of ln(data)
        let log_data: Vec<f64> = data.iter().map(|&x| x.ln()).collect();
        let expected_mu: f64 = log_data.iter().sum::<f64>() / 5.0;
        assert!(
            (mu - expected_mu).abs() < 1e-10,
            "μ = {}, expected {}",
            mu,
            expected_mu
        );
        assert!(sigma > 0.0);
    }

    #[test]
    fn fit_lognormal_edge_cases() {
        assert!(fit_lognormal(&[1.0]).is_none()); // < 2
        assert!(fit_lognormal(&[1.0, -1.0]).is_none()); // negative
        assert!(fit_lognormal(&[0.0, 1.0]).is_none()); // zero
    }

    // -----------------------------------------------------------------------
    // Poisson fitting
    // -----------------------------------------------------------------------

    #[test]
    fn fit_poisson_basic() {
        let data = [2.0, 3.0, 1.0, 4.0, 2.0, 3.0, 1.0, 2.0];
        let r = fit_poisson(&data).expect("should compute");
        assert_eq!(r.distribution, "Poisson");
        assert_eq!(r.n_params, 1);
        let mean: f64 = data.iter().sum::<f64>() / data.len() as f64;
        assert!(
            (r.parameters[0].1 - mean).abs() < 1e-10,
            "λ̂ = {}, expected {}",
            r.parameters[0].1,
            mean
        );
    }

    #[test]
    fn fit_poisson_edge_cases() {
        assert!(fit_poisson(&[1.0]).is_none()); // < 2
        assert!(fit_poisson(&[1.0, -1.0]).is_none()); // negative
        assert!(fit_poisson(&[0.0, 0.0]).is_none()); // all zeros
        assert!(fit_poisson(&[1.0, 2.5]).is_none()); // not a count
    }

    // -----------------------------------------------------------------------
    // Beta fitting
    // -----------------------------------------------------------------------

    #[test]
    fn fit_beta_basic() {
        let data = [0.2, 0.35, 0.5, 0.15, 0.45, 0.3, 0.6, 0.25, 0.4, 0.55];
        let r = fit_beta(&data).expect("should compute");
        assert_eq!(r.distribution, "Beta");
        assert_eq!(r.n_params, 2);
        assert!(r.parameters[0].1 > 0.0, "α̂ > 0");
        assert!(r.parameters[1].1 > 0.0, "β̂ > 0");
        assert!(r.log_likelihood.is_finite());
    }

    #[test]
    fn fit_beta_symmetric() {
        // Symmetric data around 0.5 → α ≈ β
        let data = [0.3, 0.4, 0.5, 0.6, 0.7, 0.35, 0.45, 0.55, 0.65, 0.5];
        let r = fit_beta(&data).expect("should compute");
        let alpha = r.parameters[0].1;
        let beta = r.parameters[1].1;
        // α and β should be roughly equal for symmetric data
        assert!(
            (alpha - beta).abs() / alpha < 0.3,
            "α = {alpha}, β = {beta}, should be roughly equal"
        );
    }

    #[test]
    fn fit_beta_edge_cases() {
        assert!(fit_beta(&[0.5]).is_none()); // < 2
        assert!(fit_beta(&[0.5, 1.0]).is_none()); // value = 1.0 (boundary)
        assert!(fit_beta(&[0.0, 0.5]).is_none()); // value = 0.0 (boundary)
        assert!(fit_beta(&[0.5, 1.5]).is_none()); // value > 1
        assert!(fit_beta(&[0.5, -0.5]).is_none()); // value < 0
    }

    // -----------------------------------------------------------------------
    // fit_best (updated)
    // -----------------------------------------------------------------------

    #[test]
    fn fit_best_compares_the_continuous_families_that_apply() {
        let data = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0];
        let fits = fit_best(&data);
        let names: Vec<&str> = fits.iter().map(|f| f.distribution.as_str()).collect();
        // All > 0 and above 1: every positive family, not Beta.
        for family in ["Normal", "Exponential", "Gamma", "LogNormal", "Weibull"] {
            assert!(
                names.contains(&family),
                "should include {family}: {names:?}"
            );
        }
        assert!(!names.contains(&"Beta"), "{names:?}");
        // Counts are still data a density is fitted to; a mass is not ranked
        // against densities.
        assert!(!names.contains(&"Poisson"), "{names:?}");

        let proportions = [0.2, 0.35, 0.5, 0.15, 0.45, 0.3, 0.6, 0.25, 0.4, 0.55];
        let names: Vec<String> = fit_best(&proportions)
            .into_iter()
            .map(|f| f.distribution)
            .collect();
        assert!(names.iter().any(|n| n == "Beta"), "{names:?}");
    }

    /// The ranking a mass and a density cannot share: rescaling continuous
    /// data moves every density's log-likelihood by n·ln(scale) and leaves the
    /// order among densities alone -- so the winner must not depend on units.
    #[test]
    fn fit_best_ranking_does_not_depend_on_the_unit() {
        let data = [2.0, 3.0, 1.0, 4.0, 2.0, 3.0, 5.0, 2.0, 3.0, 4.0];
        let in_thousandths: Vec<f64> = data.iter().map(|v| v * 1000.0).collect();
        let order = |d: &[f64]| -> Vec<String> {
            fit_best(d).into_iter().map(|f| f.distribution).collect()
        };
        assert_eq!(order(&data), order(&in_thousandths));
    }

    #[test]
    fn fit_weibull_matches_the_weibull_mle() {
        let data = [12.0, 25.0, 31.0, 44.0, 52.0, 60.0, 71.0, 85.0];
        let fit = fit_weibull(&data).expect("positive data");
        let mle = crate::weibull::weibull_mle(&data).expect("same data");
        assert_eq!(fit.parameters[0], ("shape".to_string(), mle.shape));
        assert_eq!(fit.parameters[1], ("scale".to_string(), mle.scale));
        assert!((fit.aic - (-2.0 * mle.log_likelihood + 4.0)).abs() < 1e-12);
        assert!(fit_weibull(&[1.0, 0.0, 2.0]).is_none());
    }
}

#[cfg(test)]
mod proptests {
    use super::*;
    use proptest::prelude::*;

    proptest! {
        #[test]
        fn ecdf_last_is_one(
            data in proptest::collection::vec(-1e3_f64..1e3, 1..=50)
        ) {
            if let Some((_, probs)) = ecdf(&data) {
                let last = *probs.last().expect("non-empty");
                prop_assert!((last - 1.0).abs() < 1e-10, "last prob = {last}");
            }
        }

        #[test]
        fn ecdf_monotonic(
            data in proptest::collection::vec(-1e3_f64..1e3, 2..=50)
        ) {
            if let Some((_, probs)) = ecdf(&data) {
                for w in probs.windows(2) {
                    prop_assert!(w[1] >= w[0] - 1e-15, "{} < {}", w[1], w[0]);
                }
            }
        }

        #[test]
        fn hist_counts_sum_to_n(
            data in proptest::collection::vec(-1e3_f64..1e3, 5..=100)
        ) {
            if let Some(r) = histogram_bins(&data, BinMethod::Sturges) {
                let total: usize = r.counts.iter().sum();
                prop_assert_eq!(total, data.len());
            }
        }

        #[test]
        fn kde_density_non_negative(
            data in proptest::collection::vec(-100.0_f64..100.0, 5..=50)
        ) {
            if let Some(r) = kde(&data, BandwidthMethod::Silverman, 128) {
                for &d in &r.density {
                    prop_assert!(d >= 0.0, "negative density: {d}");
                }
            }
        }

        #[test]
        fn kde_integral_approx_one(
            data in proptest::collection::vec(-100.0_f64..100.0, 10..=80)
        ) {
            if let Some(r) = kde(&data, BandwidthMethod::Silverman, 512) {
                let dx = r.x[1] - r.x[0];
                let integral: f64 = r.density.iter().sum::<f64>() * dx;
                prop_assert!(
                    (integral - 1.0).abs() < 0.1,
                    "integral = {integral}, expected ≈ 1.0"
                );
            }
        }

        #[test]
        fn fit_normal_aic_finite(
            data in proptest::collection::vec(-100.0_f64..100.0, 5..=50)
        ) {
            if let Some(r) = fit_normal(&data) {
                prop_assert!(r.aic.is_finite(), "AIC = {}", r.aic);
                prop_assert!(r.bic.is_finite(), "BIC = {}", r.bic);
                prop_assert!(r.log_likelihood.is_finite(), "LL = {}", r.log_likelihood);
            }
        }

        #[test]
        fn fit_gamma_alpha_positive(
            data in proptest::collection::vec(0.1_f64..100.0, 5..=50)
        ) {
            if let Some(r) = fit_gamma(&data) {
                let alpha = r.parameters[0].1;
                let beta = r.parameters[1].1;
                prop_assert!(alpha > 0.0, "α = {alpha}");
                prop_assert!(beta > 0.0, "β = {beta}");
                prop_assert!(r.aic.is_finite(), "AIC = {}", r.aic);
            }
        }
    }
}

/// `fit_best` against the maximum-likelihood fits of an independent
/// computation: closed forms for Normal, Exponential and LogNormal, and the
/// score equations solved to 1e-15 (SciPy `brentq` / `root`, digamma from
/// SciPy) for Gamma, Weibull and Beta, with the log-likelihood summed from
/// SciPy's own `logpdf`. `scipy.stats.<family>.fit` is not the reference — it
/// stops short of the optimum (Weibull, Gamma).
#[cfg(test)]
mod scipy_reference {
    use super::*;

    type Want = (&'static str, &'static [f64], f64, f64, f64);

    fn close(got: f64, want: f64) -> bool {
        (got - want).abs() <= 1e-9 * want.abs().max(1.0)
    }

    fn check(data: &[f64], want: &[Want], what: &str) {
        let got = fit_best(data);
        let names: Vec<&str> = got.iter().map(|f| f.distribution.as_str()).collect();
        let wanted: Vec<&str> = want.iter().map(|w| w.0).collect();
        assert_eq!(names, wanted, "{what}: families in AIC order");
        for (f, w) in got.iter().zip(want) {
            let name = &f.distribution;
            assert_eq!(f.parameters.len(), w.1.len(), "{what} {name}: parameters");
            for ((pname, p), q) in f.parameters.iter().zip(w.1) {
                assert!(
                    close(*p, *q),
                    "{what} {name} {pname}: got {p}, reference {q}"
                );
            }
            assert!(
                close(f.log_likelihood, w.2),
                "{what} {name} log-likelihood: got {}, reference {}",
                f.log_likelihood,
                w.2
            );
            assert!(
                close(f.aic, w.3),
                "{what} {name} AIC: got {}, reference {}",
                f.aic,
                w.3
            );
            assert!(
                close(f.bic, w.4),
                "{what} {name} BIC: got {}, reference {}",
                f.bic,
                w.4
            );
        }
    }

    #[test]
    fn data_inside_the_unit_interval_is_fitted_by_all_six_families() {
        check(
            &[
                0.12, 0.23, 0.31, 0.38, 0.44, 0.52, 0.57, 0.63, 0.71, 0.78, 0.86, 0.93,
            ],
            &[
                (
                    "Beta",
                    &[1.8438565181620945, 1.5583703583443078],
                    1.035200520768992,
                    1.929598958462016,
                    2.8994122580380166,
                ),
                (
                    "Weibull",
                    &[2.391258961247938, 0.6091890535880672],
                    0.10397825424951779,
                    3.7920434915009644,
                    4.761856791076966,
                ),
                (
                    "Normal",
                    &[0.54, 0.24348168993444524],
                    -0.0746999829213777,
                    4.149399965842756,
                    5.119213265418756,
                ),
                (
                    "Gamma",
                    &[3.7996792156263166, 7.036442991900586],
                    -0.4996467398561213,
                    4.999293479712243,
                    5.969106779288243,
                ),
                (
                    "LogNormal",
                    &[-0.7535094301720809, 0.5764955804748335],
                    -1.3756979711767767,
                    6.751395942353554,
                    7.721209241929554,
                ),
                (
                    "Exponential",
                    &[1.8518518518518516],
                    -4.605766326914196,
                    11.211532653828392,
                    11.696439303616392,
                ),
            ],
            "unit interval",
        );
    }

    #[test]
    fn positive_data_is_fitted_by_the_five_positive_families() {
        check(
            &[2.1, 3.5, 1.8, 4.2, 2.9, 3.1, 2.5, 3.8, 1.5, 2.7, 6.4, 9.8],
            &[
                (
                    "LogNormal",
                    &[1.1668644506389008, 0.5022892551220499],
                    -22.766686368395682,
                    49.533372736791364,
                    50.50318603636737,
                ),
                (
                    "Gamma",
                    &[3.7501260541767762, 1.0158355000027386],
                    -23.629758495985087,
                    51.25951699197017,
                    52.22933029154618,
                ),
                (
                    "Weibull",
                    &[1.8208938249614839, 4.190960652455631],
                    -24.62416845020696,
                    53.24833690041392,
                    54.21815019998992,
                ),
                (
                    "Normal",
                    &[3.6916666666666664, 2.2220329780531065],
                    -26.60833279104041,
                    57.21666558208082,
                    58.186478881656825,
                ),
                (
                    "Exponential",
                    &[0.2708803611738149],
                    -27.672936327157075,
                    57.34587265431415,
                    57.83077930410215,
                ),
            ],
            "positive",
        );
    }

    #[test]
    fn data_of_mixed_sign_is_fitted_by_the_normal_alone() {
        check(
            &[-1.2, 0.4, 2.3, -0.7, 1.9, 0.8, -2.1, 3.1],
            &[(
                "Normal",
                &[0.5625, 1.6999540434964704],
                -15.596318005312547,
                35.192636010625094,
                35.35151909398476,
            )],
            "mixed sign",
        );
    }
}

/// Goodness of fit against SciPy 1.17.1 (`ad_ref.py` / `p_ref.py` in the
/// cycle log): A² from `logcdf` / `logsf` of the same distribution, and the
/// parametric bootstrap of `scipy.stats.goodness_of_fit(statistic="ad")`.
#[cfg(test)]
mod goodness_of_fit_reference {
    use super::*;

    const POS: &[f64] = &[
        3.2, 4.8, 5.1, 6.0, 6.7, 7.3, 8.1, 9.4, 10.2, 12.5, 15.8, 41.0,
    ];
    const UNIT: &[f64] = &[
        0.04, 0.11, 0.19, 0.23, 0.31, 0.38, 0.42, 0.55, 0.61, 0.74, 0.88, 0.97,
    ];

    fn fixed(distribution: &str, parameters: &[(&str, f64)]) -> FitResult {
        FitResult {
            distribution: distribution.to_string(),
            parameters: parameters
                .iter()
                .map(|&(n, v)| (n.to_string(), v))
                .collect(),
            log_likelihood: 0.0,
            aic: 0.0,
            bic: 0.0,
            n_params: parameters.len(),
        }
    }

    /// Every family's A² against fixed parameters — including a value (41) in
    /// the far upper tail of the Normal and Exponential, where `1 − F` taken by
    /// subtraction would lose its digits.
    #[test]
    fn anderson_darling_statistic_reproduces_scipy_for_every_family() {
        type Case = (
            &'static str,
            &'static [(&'static str, f64)],
            &'static [f64],
            f64,
        );
        let cases: [Case; 6] = [
            (
                "Normal",
                &[("mu", 9.0), ("sigma", 4.0)],
                POS,
                3.0398709332150187,
            ),
            (
                "LogNormal",
                &[("mu", 2.0), ("sigma", 0.6)],
                POS,
                0.4154008172250965,
            ),
            (
                "Weibull",
                &[("shape", 1.7), ("scale", 11.0)],
                POS,
                0.8259972063661056,
            ),
            (
                "Gamma",
                &[("alpha", 2.5), ("beta", 0.25)],
                POS,
                0.5575965468003456,
            ),
            ("Exponential", &[("lambda", 0.1)], POS, 1.235635893786549),
            (
                "Beta",
                &[("alpha", 1.3), ("beta", 1.6)],
                UNIT,
                0.34159812117638744,
            ),
        ];
        for (name, params, data, want) in cases {
            let got = anderson_darling_fit(data, &fixed(name, params)).unwrap();
            assert!(
                ((got - want) / want).abs() < 1e-9,
                "{name}: A² {got}, SciPy {want}"
            );
        }
    }

    /// The bootstrap p-value agrees with SciPy's Monte Carlo of the same test
    /// (99 999 samples there, 9 999 here: standard errors 0.001 and 0.004).
    /// Exponential only: SciPy's Normal fit takes the n − 1 standard deviation
    /// where `fit_normal` takes the maximum-likelihood one, so its A² (0.4988)
    /// is a different test's statistic (A² at the MLE fit: 0.5086, as here).
    #[test]
    fn bootstrap_p_value_agrees_with_scipy_goodness_of_fit() {
        let data = [
            3.2, 4.8, 5.1, 6.0, 6.7, 7.3, 8.1, 9.4, 10.2, 12.5, 15.8, 21.0,
        ];
        let fit = fit_exponential(&data).unwrap();
        let (statistic, p) = (1.4502868637902981, 0.03158);
        let g = goodness_of_fit(&data, &fit, 9999, 3).unwrap();
        assert!(
            ((g.statistic - statistic) / statistic).abs() < 1e-9,
            "A² {} vs {statistic}",
            g.statistic
        );
        assert!(
            (g.p_value - p).abs() < 0.015,
            "p {} vs SciPy {p}",
            g.p_value
        );
    }

    /// The case that motivated the test: uniform data, whose best-AIC family
    /// places the median at 4400 of 1..10 000 — no listed family fits it, and
    /// the test says so; a sample drawn from a family is not rejected.
    #[test]
    fn rejects_the_best_family_for_uniform_data_and_keeps_a_true_one() {
        use u_numflow::distributions::Sample;
        let uniform: Vec<f64> = (1..=2000).map(f64::from).collect();
        let best = &fit_best(&uniform)[0];
        let g = goodness_of_fit(&uniform, best, DEFAULT_GOF_REPLICATES, 0).unwrap();
        assert_eq!(g.p_value, 1.0 / 200.0, "{}: {g:?}", best.distribution);

        let weibull = Weibull::new(2.0, 10.0).unwrap().sample_n(300, 11);
        let fit = fit_weibull(&weibull).unwrap();
        let g = goodness_of_fit(&weibull, &fit, DEFAULT_GOF_REPLICATES, 0).unwrap();
        assert!(g.p_value > 0.05, "{g:?}");
    }

    #[test]
    fn same_seed_same_p_value_and_refusals() {
        let data = [2.1, 3.5, 1.8, 4.2, 2.9, 3.1, 2.5, 3.8, 1.5, 2.7];
        let fit = fit_gamma(&data).unwrap();
        let a = goodness_of_fit(&data, &fit, 99, 7).unwrap();
        assert_eq!(a, goodness_of_fit(&data, &fit, 99, 7).unwrap());
        assert_eq!(a.replicates, 99);
        assert!(goodness_of_fit(&data, &fit, 0, 7).is_none());
        assert!(goodness_of_fit(&[1.0], &fit, 9, 7).is_none());
        assert!(goodness_of_fit(&[1.0, f64::NAN, 2.0], &fit, 9, 7).is_none());
        let poisson = fit_poisson(&[1.0, 2.0, 2.0, 3.0]).unwrap();
        assert!(anderson_darling_fit(&data, &poisson).is_none());
    }

    /// The statistic and the p-value do not depend on the unit of the data.
    #[test]
    fn unit_invariant() {
        let data = [2.1, 3.5, 1.8, 4.2, 2.9, 3.1, 2.5, 3.8, 1.5, 2.7, 6.3, 2.2];
        for scale in [1e-9, 1.0, 1e9] {
            let x: Vec<f64> = data.iter().map(|v| v * scale).collect();
            for fit in fit_best(&x) {
                let base: Vec<f64> = data.to_vec();
                let reference = fit_best(&base)
                    .into_iter()
                    .find(|f| f.distribution == fit.distribution);
                let Some(reference) = reference else { continue };
                let a = goodness_of_fit(&x, &fit, 49, 1).unwrap();
                let b = goodness_of_fit(&base, &reference, 49, 1).unwrap();
                assert!(
                    (a.statistic - b.statistic).abs() < 1e-6 * b.statistic.max(1.0),
                    "{} at ×{scale}: {} vs {}",
                    fit.distribution,
                    a.statistic,
                    b.statistic
                );
            }
        }
    }
}