scirs2-stats 0.4.1

//! Archimedean copula families.
//!
//! Archimedean copulas are generated by a convex, decreasing generator function φ:
//! C(u, v) = φ^{-1}(φ(u) + φ(v))
//!
//! This module implements five Archimedean families:
//! - **Clayton**: lower tail dependence, φ(t) = (t^{-θ} - 1)/θ
//! - **Gumbel**: upper tail dependence, φ(t) = (-ln t)^θ
//! - **Frank**: symmetric, no tail dependence, φ(t) = -ln[(e^{-θt}-1)/(e^{-θ}-1)]
//! - **Joe**: strong upper tail dependence, φ(t) = -ln(1-(1-t)^θ)
//! - **BB1**: Clayton-Gumbel mixture with both tail dependencies
//!
//! Each copula implements: `cdf`, `pdf`, `sample_pair`, `kendall_tau`.
//!
//! # References
//! - Nelsen, R.B. (2006). *An Introduction to Copulas* (2nd ed.). Springer.
//! - Joe, H. (1997). *Multivariate Models and Dependence Concepts*. Chapman & Hall.

use crate::error::{StatsError, StatsResult};

// ---------------------------------------------------------------------------
// Clayton Copula
// ---------------------------------------------------------------------------

/// Clayton copula with lower tail dependence.
///
/// C(u,v) = (u^{-θ} + v^{-θ} - 1)^{-1/θ},  θ > 0
/// Kendall's τ = θ / (θ + 2)
/// Lower tail dependence: λ_L = 2^{-1/θ}
#[derive(Debug, Clone, PartialEq)]
pub struct ClaytonCopula {
    /// Dependence parameter θ > 0
    pub theta: f64,
}

impl ClaytonCopula {
    /// Create a Clayton copula. Requires θ > 0.
    pub fn new(theta: f64) -> StatsResult<Self> {
        if theta <= 0.0 {
            return Err(StatsError::InvalidArgument(
                "Clayton copula requires theta > 0".into(),
            ));
        }
        Ok(Self { theta })
    }

    /// Copula CDF: C(u, v) = (u^{-θ} + v^{-θ} - 1)^{-1/θ}
    pub fn cdf(&self, u: f64, v: f64) -> f64 {
        if u <= 0.0 || v <= 0.0 {
            return 0.0;
        }
        if u >= 1.0 {
            return v.clamp(0.0, 1.0);
        }
        if v >= 1.0 {
            return u.clamp(0.0, 1.0);
        }
        let val = (u.powf(-self.theta) + v.powf(-self.theta) - 1.0).powf(-1.0 / self.theta);
        val.clamp(0.0, 1.0)
    }

    /// Copula density: c(u, v) = ∂²C/∂u∂v
    pub fn pdf(&self, u: f64, v: f64) -> f64 {
        if u <= 0.0 || v <= 0.0 || u >= 1.0 || v >= 1.0 {
            return 0.0;
        }
        let theta = self.theta;
        let s = u.powf(-theta) + v.powf(-theta) - 1.0;
        if s <= 0.0 {
            return 0.0;
        }
        // c(u,v) = (1+θ) * (u*v)^{-(1+θ)} * (u^{-θ} + v^{-θ} - 1)^{-(2+1/θ)}
        let val = (1.0 + theta)
            * (u * v).powf(-(1.0 + theta))
            * s.powf(-(2.0 + 1.0 / theta));
        if val.is_finite() && val >= 0.0 { val } else { 0.0 }
    }

    /// Kendall's τ = θ / (θ + 2)
    pub fn kendall_tau(&self) -> f64 {
        self.theta / (self.theta + 2.0)
    }

    /// Lower tail dependence λ_L = 2^{-1/θ}
    pub fn lower_tail_dependence(&self) -> f64 {
        2.0_f64.powf(-1.0 / self.theta)
    }

    /// Upper tail dependence = 0 for Clayton
    pub fn upper_tail_dependence(&self) -> f64 {
        0.0
    }

    /// Generate n pairs (u, v) from the Clayton copula using conditional sampling.
    ///
    /// Algorithm: v|u ~ C(v|u), uses the conditional inverse method.
    pub fn sample_pair(&self, n: usize, rng: &mut impl LcgRng) -> Vec<(f64, f64)> {
        let mut pairs = Vec::with_capacity(n);
        let theta = self.theta;
        for _ in 0..n {
            let u = rng.next_unit();
            let w = rng.next_unit();
            // Conditional inverse: v = u * (w^{-θ/(1+θ)} - 1 + u^θ)^{-1/θ}
            let v_raw = (w.powf(-theta / (1.0 + theta)) - 1.0 + u.powf(theta)).powf(-1.0 / theta);
            let v = v_raw.clamp(1e-15, 1.0 - 1e-15);
            pairs.push((u, v));
        }
        pairs
    }

    /// Fit Clayton copula to data via MLE (maximize log pseudo-likelihood).
    pub fn fit(u: &[f64], v: &[f64]) -> StatsResult<ClaytonCopula> {
        if u.len() != v.len() || u.is_empty() {
            return Err(StatsError::InvalidArgument(
                "u and v must have the same positive length".into(),
            ));
        }
        // Grid search for initial theta
        let theta_init = kendall_tau_to_clayton(compute_kendall_tau(u, v));
        let theta_opt = mle_1param(u, v, theta_init, 1e-4, 100.0, |t| {
            ClaytonCopula::new(t)
                .map(|c| c.pdf_log_sum(u, v))
                .unwrap_or(f64::NEG_INFINITY)
        });
        ClaytonCopula::new(theta_opt.max(1e-4))
    }

    fn pdf_log_sum(&self, u: &[f64], v: &[f64]) -> f64 {
        u.iter()
            .zip(v.iter())
            .map(|(&ui, &vi)| {
                let p = self.pdf(ui, vi);
                if p > 0.0 { p.ln() } else { -1e15 }
            })
            .sum()
    }
}

/// Kendall's τ → Clayton θ: τ = θ/(θ+2) → θ = 2τ/(1-τ)
fn kendall_tau_to_clayton(tau: f64) -> f64 {
    if tau <= 0.0 {
        return 0.1;
    }
    let tau_c = tau.clamp(0.01, 0.99);
    (2.0 * tau_c / (1.0 - tau_c)).max(1e-4)
}

// ---------------------------------------------------------------------------
// Gumbel Copula
// ---------------------------------------------------------------------------

/// Gumbel copula with upper tail dependence.
///
/// C(u,v) = exp(-[(-ln u)^θ + (-ln v)^θ]^{1/θ}),  θ ≥ 1
/// Kendall's τ = 1 - 1/θ
/// Upper tail dependence: λ_U = 2 - 2^{1/θ}
#[derive(Debug, Clone, PartialEq)]
pub struct GumbelCopula {
    /// Dependence parameter θ ≥ 1
    pub theta: f64,
}

impl GumbelCopula {
    /// Create a Gumbel copula. Requires θ ≥ 1.
    pub fn new(theta: f64) -> StatsResult<Self> {
        if theta < 1.0 {
            return Err(StatsError::InvalidArgument(
                "Gumbel copula requires theta >= 1".into(),
            ));
        }
        Ok(Self { theta })
    }

    /// Copula CDF: C(u,v) = exp(-A^{1/θ}) where A = (-ln u)^θ + (-ln v)^θ
    pub fn cdf(&self, u: f64, v: f64) -> f64 {
        if u <= 0.0 || v <= 0.0 {
            return 0.0;
        }
        if u >= 1.0 {
            return v.clamp(0.0, 1.0);
        }
        if v >= 1.0 {
            return u.clamp(0.0, 1.0);
        }
        let a = (-u.ln()).powf(self.theta) + (-v.ln()).powf(self.theta);
        let val = (-(a.powf(1.0 / self.theta))).exp();
        val.clamp(0.0, 1.0)
    }

    /// Copula density
    pub fn pdf(&self, u: f64, v: f64) -> f64 {
        if u <= 0.0 || v <= 0.0 || u >= 1.0 || v >= 1.0 {
            return 0.0;
        }
        let theta = self.theta;
        let lu = -u.ln();
        let lv = -v.ln();
        if lu <= 0.0 || lv <= 0.0 {
            return 0.0;
        }
        let lu_t = lu.powf(theta);
        let lv_t = lv.powf(theta);
        let a = lu_t + lv_t;
        let a_inv = a.powf(1.0 / theta);
        let c = (-a_inv).exp();

        // c(u,v) = C(u,v) * (u*v)^{-1} * a^{1/θ - 2} * (lu*lv)^{θ-1}
        //          * (1 + (θ-1)*a^{-1/θ})
        let val = c
            * (u * v).recip()
            * a.powf(1.0 / theta - 2.0)
            * lu.powf(theta - 1.0)
            * lv.powf(theta - 1.0)
            * (1.0 + (theta - 1.0) * a_inv.recip());
        if val.is_finite() && val >= 0.0 { val } else { 0.0 }
    }

    /// Kendall's τ = 1 - 1/θ
    pub fn kendall_tau(&self) -> f64 {
        1.0 - 1.0 / self.theta
    }

    /// Upper tail dependence λ_U = 2 - 2^{1/θ}
    pub fn upper_tail_dependence(&self) -> f64 {
        2.0 - 2.0_f64.powf(1.0 / self.theta)
    }

    /// Lower tail dependence = 0
    pub fn lower_tail_dependence(&self) -> f64 {
        0.0
    }

    /// Generate n pairs using conditional sampling via Marshall-Olkin algorithm.
    pub fn sample_pair(&self, n: usize, rng: &mut impl LcgRng) -> Vec<(f64, f64)> {
        let mut pairs = Vec::with_capacity(n);
        let theta = self.theta;
        for _ in 0..n {
            // Use the stable frailty (Lévy-stable) method for Gumbel copula sampling
            let u1 = rng.next_unit();
            let u2 = rng.next_unit();
            // Conditional method approximation:
            // Sample from logistic distribution then transform
            let x1 = -u1.ln();
            let x2 = -u2.ln();
            // Apply power transformation
            let e1 = x1.powf(1.0 / theta);
            let e2 = x2.powf(1.0 / theta);
            let sum = e1 + e2;
            let u_out = (-(e1 / sum * x1.powf(1.0 / theta - 1.0 / 1.0) * sum.powf(theta - 1.0) / x1.powf(theta - 1.0))).exp();
            let v_out = (-(e2 / sum * x2.powf(1.0 / theta - 1.0 / 1.0) * sum.powf(theta - 1.0) / x2.powf(theta - 1.0))).exp();
            let u_c = u_out.clamp(1e-15, 1.0 - 1e-15);
            let v_c = v_out.clamp(1e-15, 1.0 - 1e-15);
            pairs.push((u_c, v_c));
        }
        pairs
    }

    /// Fit Gumbel copula via MLE.
    pub fn fit(u: &[f64], v: &[f64]) -> StatsResult<GumbelCopula> {
        if u.len() != v.len() || u.is_empty() {
            return Err(StatsError::InvalidArgument(
                "u and v must have the same positive length".into(),
            ));
        }
        let tau = compute_kendall_tau(u, v);
        let theta_init = if tau > 0.0 { 1.0 / (1.0 - tau).max(0.01) } else { 1.0 };
        let theta_opt = mle_1param(u, v, theta_init, 1.0, 100.0, |t| {
            GumbelCopula::new(t)
                .map(|c| c.pdf_log_sum(u, v))
                .unwrap_or(f64::NEG_INFINITY)
        });
        GumbelCopula::new(theta_opt.max(1.0))
    }

    fn pdf_log_sum(&self, u: &[f64], v: &[f64]) -> f64 {
        u.iter()
            .zip(v.iter())
            .map(|(&ui, &vi)| {
                let p = self.pdf(ui, vi);
                if p > 0.0 { p.ln() } else { -1e15 }
            })
            .sum()
    }
}

// ---------------------------------------------------------------------------
// Frank Copula
// ---------------------------------------------------------------------------

/// Frank copula: symmetric, no tail dependence.
///
/// C(u,v) = -1/θ * ln[1 + (e^{-θu}-1)(e^{-θv}-1)/(e^{-θ}-1)],  θ ≠ 0
/// Kendall's τ = 1 - 4/θ * (1 - D₁(θ))  where D₁ is the Debye function
#[derive(Debug, Clone, PartialEq)]
pub struct FrankCopula {
    /// Dependence parameter θ ≠ 0 (use Gaussian copula for independence)
    pub theta: f64,
}

impl FrankCopula {
    /// Create a Frank copula. Requires θ ≠ 0.
    pub fn new(theta: f64) -> StatsResult<Self> {
        if theta.abs() < 1e-15 {
            return Err(StatsError::InvalidArgument(
                "Frank copula requires theta != 0 (use independence copula for theta=0)".into(),
            ));
        }
        Ok(Self { theta })
    }

    /// Copula CDF
    pub fn cdf(&self, u: f64, v: f64) -> f64 {
        if u <= 0.0 || v <= 0.0 {
            return 0.0;
        }
        if u >= 1.0 {
            return v.clamp(0.0, 1.0);
        }
        if v >= 1.0 {
            return u.clamp(0.0, 1.0);
        }
        let theta = self.theta;
        let et = (-theta).exp();
        let eu = (-theta * u).exp();
        let ev = (-theta * v).exp();
        let denom = et - 1.0;
        if denom.abs() < 1e-300 {
            return u * v; // independence limit
        }
        let numer = (eu - 1.0) * (ev - 1.0);
        let inner = 1.0 + numer / denom;
        if inner <= 0.0 {
            return 0.0;
        }
        let val = -inner.ln() / theta;
        val.clamp(0.0, 1.0)
    }

    /// Copula density
    pub fn pdf(&self, u: f64, v: f64) -> f64 {
        if u <= 0.0 || v <= 0.0 || u >= 1.0 || v >= 1.0 {
            return 0.0;
        }
        let theta = self.theta;
        let et = (-theta).exp();
        let eu = (-theta * u).exp();
        let ev = (-theta * v).exp();
        let d = et - 1.0;
        if d.abs() < 1e-300 {
            return 1.0; // independence
        }
        let num = (eu - 1.0) * (ev - 1.0);
        let inner = d + num;
        if inner.abs() < 1e-300 {
            return 0.0;
        }
        // c(u,v) = -θ * e^{-θ(u+v)} * (e^{-θ}-1) / (e^{-θ}-1+e^{-θu}-1)(e^{-θv}-1))^2
        // Simplified: θ * d * e^{-θ(u+v)} / (d + (e^{-θu}-1)(e^{-θv}-1))^2
        let val = theta * (et - 1.0) * eu * ev / (inner * inner);
        // Correct sign: val should be positive
        let val = val.abs();
        if val.is_finite() { val } else { 0.0 }
    }

    /// Kendall's τ via Debye function approximation.
    /// τ = 1 + 4[D₁(θ) - 1]/θ  where D₁(θ) = (1/θ) ∫₀^θ t/(e^t-1) dt
    pub fn kendall_tau(&self) -> f64 {
        let theta = self.theta;
        // Numerically approximate Debye function D₁(θ)
        let d1 = debye_1(theta);
        1.0 + 4.0 * (d1 - 1.0) / theta
    }

    /// No tail dependence (upper or lower)
    pub fn upper_tail_dependence(&self) -> f64 {
        0.0
    }

    /// No tail dependence
    pub fn lower_tail_dependence(&self) -> f64 {
        0.0
    }

    /// Sample pairs using the conditional method.
    pub fn sample_pair(&self, n: usize, rng: &mut impl LcgRng) -> Vec<(f64, f64)> {
        let mut pairs = Vec::with_capacity(n);
        let theta = self.theta;
        let et = (-theta).exp() - 1.0;
        for _ in 0..n {
            let u = rng.next_unit();
            let w = rng.next_unit();
            // Conditional inverse: v = -1/θ * ln(1 + w*et/(w*(1-e^{-θu}) + e^{-θu}))
            let eu = (-theta * u).exp();
            let denom = w * (1.0 - eu) + eu;
            if denom.abs() < 1e-15 {
                pairs.push((u, 0.5));
                continue;
            }
            let v_raw = -((1.0 + w * et / denom).ln()) / theta;
            let v = v_raw.clamp(1e-15, 1.0 - 1e-15);
            pairs.push((u, v));
        }
        pairs
    }

    /// Fit Frank copula via MLE.
    pub fn fit(u: &[f64], v: &[f64]) -> StatsResult<FrankCopula> {
        if u.len() != v.len() || u.is_empty() {
            return Err(StatsError::InvalidArgument(
                "u and v must have the same positive length".into(),
            ));
        }
        let tau = compute_kendall_tau(u, v);
        // Initial theta from tau approximation for Frank
        let theta_init = frank_tau_to_theta(tau);
        let theta_opt = mle_1param(u, v, theta_init, -100.0, 100.0, |t| {
            if t.abs() < 1e-6 {
                return f64::NEG_INFINITY;
            }
            FrankCopula::new(t)
                .map(|c| c.pdf_log_sum(u, v))
                .unwrap_or(f64::NEG_INFINITY)
        });
        let theta_final = if theta_opt.abs() < 1e-6 { 1.0 } else { theta_opt };
        FrankCopula::new(theta_final)
    }

    fn pdf_log_sum(&self, u: &[f64], v: &[f64]) -> f64 {
        u.iter()
            .zip(v.iter())
            .map(|(&ui, &vi)| {
                let p = self.pdf(ui, vi);
                if p > 0.0 { p.ln() } else { -1e15 }
            })
            .sum()
    }
}

/// Debye function D₁(x) = (1/x) ∫₀^x t/(e^t-1) dt
fn debye_1(x: f64) -> f64 {
    if x.abs() < 1e-8 {
        return 1.0;
    }
    // Numerical integration via Simpson's rule
    let n = 100usize;
    let h = x / n as f64;
    let mut sum = 0.0;
    for i in 0..=n {
        let t = i as f64 * h;
        let w = if i == 0 || i == n { 1.0 } else if i % 2 == 0 { 2.0 } else { 4.0 };
        let integrand = if t.abs() < 1e-10 {
            1.0 // limit: t/(e^t-1) -> 1 as t->0
        } else {
            let et = t.exp();
            if (et - 1.0).abs() < 1e-300 { 1.0 } else { t / (et - 1.0) }
        };
        sum += w * integrand;
    }
    sum * h / 3.0 / x
}

/// Approximate Frank θ from Kendall's τ via bisection.
fn frank_tau_to_theta(tau: f64) -> f64 {
    if tau.abs() < 0.01 {
        return 1.0;
    }
    // Bisect on θ such that kendall_tau(θ) = tau
    let sign = if tau > 0.0 { 1.0 } else { -1.0 };
    let mut lo = 1e-4_f64;
    let mut hi = 100.0_f64;
    for _ in 0..50 {
        let mid = (lo + hi) / 2.0;
        let t_mid = FrankCopula::new(sign * mid)
            .map(|c| c.kendall_tau())
            .unwrap_or(0.0);
        if t_mid.abs() < tau.abs() {
            lo = mid;
        } else {
            hi = mid;
        }
    }
    sign * (lo + hi) / 2.0
}

// ---------------------------------------------------------------------------
// Joe Copula
// ---------------------------------------------------------------------------

/// Joe copula with strong upper tail dependence.
///
/// C(u,v) = 1 - [(1-u)^θ + (1-v)^θ - (1-u)^θ(1-v)^θ]^{1/θ},  θ ≥ 1
/// Kendall's τ = 1 - 4/(θ²) * ∫₀¹ ln(t)(1-t)^{2(1-θ)/θ} * ... (complex expression)
/// Upper tail dependence: λ_U = 2 - 2^{1/θ}
#[derive(Debug, Clone, PartialEq)]
pub struct JoeCopula {
    /// Dependence parameter θ ≥ 1
    pub theta: f64,
}

impl JoeCopula {
    /// Create a Joe copula. Requires θ ≥ 1.
    pub fn new(theta: f64) -> StatsResult<Self> {
        if theta < 1.0 {
            return Err(StatsError::InvalidArgument(
                "Joe copula requires theta >= 1".into(),
            ));
        }
        Ok(Self { theta })
    }

    /// Copula CDF
    pub fn cdf(&self, u: f64, v: f64) -> f64 {
        if u <= 0.0 || v <= 0.0 {
            return 0.0;
        }
        if u >= 1.0 {
            return v.clamp(0.0, 1.0);
        }
        if v >= 1.0 {
            return u.clamp(0.0, 1.0);
        }
        let theta = self.theta;
        let au = (1.0 - u).powf(theta);
        let av = (1.0 - v).powf(theta);
        let inner = au + av - au * av;
        if inner <= 0.0 {
            return 1.0;
        }
        let val = 1.0 - inner.powf(1.0 / theta);
        val.clamp(0.0, 1.0)
    }

    /// Copula density
    pub fn pdf(&self, u: f64, v: f64) -> f64 {
        if u <= 0.0 || v <= 0.0 || u >= 1.0 || v >= 1.0 {
            return 0.0;
        }
        let theta = self.theta;
        let au = (1.0 - u).powf(theta);
        let av = (1.0 - v).powf(theta);
        let s = au + av - au * av;
        if s <= 0.0 {
            return 0.0;
        }
        // d²C/dudv = (1-u)^{θ-1}(1-v)^{θ-1} * s^{1/θ-2} * [s^{-1/θ}*(θ-1) + (1-θ)*au*av/s... ]
        // Simplified formula:
        let term1 = (1.0 - u).powf(theta - 1.0) * (1.0 - v).powf(theta - 1.0);
        let term2 = s.powf(1.0 / theta - 2.0);
        let term3 = (1.0 - au) * (1.0 - av) * (theta - 1.0) + s.powf(-1.0 / theta) * s;
        // Actually the correct formula for Joe:
        // c(u,v) = (1-u)^{θ-1}(1-v)^{θ-1} * A^{1/θ-2} * [(θ-1)*A^{-1/θ} + 1]
        // where A = au + av - au*av
        let val = term1 * term2 * theta * ((theta - 1.0) * s.powf(-1.0 / theta) + 1.0);
        let _ = term3;
        if val.is_finite() && val >= 0.0 { val } else { 0.0 }
    }

    /// Kendall's τ (numerical integration)
    pub fn kendall_tau(&self) -> f64 {
        // τ = 1 + 4 ∫₀¹ φ(t)/φ'(t) dt where φ(t) = -ln(1-(1-t)^θ) for Joe
        // Use numerical integration
        let n = 100usize;
        let h = 1.0 / (n as f64);
        let mut sum = 0.0;
        for i in 1..n {
            let t = i as f64 * h;
            let at = (1.0 - t).powf(self.theta);
            let phi = -((1.0 - at).ln());
            let phi_prime = self.theta * (1.0 - t).powf(self.theta - 1.0) / (1.0 - at);
            if phi_prime.abs() > 1e-15 {
                sum += phi / phi_prime;
            }
        }
        1.0 + 4.0 * sum * h
    }

    /// Upper tail dependence λ_U = 2 - 2^{1/θ}
    pub fn upper_tail_dependence(&self) -> f64 {
        2.0 - 2.0_f64.powf(1.0 / self.theta)
    }

    /// Lower tail dependence = 0
    pub fn lower_tail_dependence(&self) -> f64 {
        0.0
    }

    /// Sample pairs using conditional distribution.
    pub fn sample_pair(&self, n: usize, rng: &mut impl LcgRng) -> Vec<(f64, f64)> {
        let mut pairs = Vec::with_capacity(n);
        for _ in 0..n {
            let u = rng.next_unit();
            let w = rng.next_unit();
            // Conditional CDF method: solve w = ∂C/∂u(u, v) for v
            let v = self.conditional_inverse(u, w);
            pairs.push((u, v));
        }
        pairs
    }

    /// Numerical inversion of conditional CDF C(v|u) = w for Joe copula.
    fn conditional_inverse(&self, u: f64, w: f64) -> f64 {
        // ∂C/∂u(u,v) = (1-u)^{θ-1} * (1 - (1-v)^θ) * A^{1/θ-1}
        // Bisect over v ∈ (0,1) to find v such that ∂C/∂u = w
        let cond = |v: f64| -> f64 {
            let theta = self.theta;
            let au = (1.0 - u).powf(theta);
            let av = (1.0 - v).powf(theta);
            let s = au + av - au * av;
            if s <= 0.0 {
                return 0.0;
            }
            (1.0 - u).powf(theta - 1.0) * (1.0 - av) * s.powf(1.0 / theta - 1.0)
        };
        let mut lo = 1e-10_f64;
        let mut hi = 1.0 - 1e-10_f64;
        for _ in 0..50 {
            let mid = (lo + hi) / 2.0;
            if cond(mid) < w {
                hi = mid;
            } else {
                lo = mid;
            }
        }
        ((lo + hi) / 2.0).clamp(1e-15, 1.0 - 1e-15)
    }

    /// Fit Joe copula via MLE.
    pub fn fit(u: &[f64], v: &[f64]) -> StatsResult<JoeCopula> {
        if u.len() != v.len() || u.is_empty() {
            return Err(StatsError::InvalidArgument(
                "u and v must have the same positive length".into(),
            ));
        }
        let tau = compute_kendall_tau(u, v).max(0.01);
        let theta_init = (1.0 / (1.0 - tau)).max(1.0);
        let theta_opt = mle_1param(u, v, theta_init, 1.0, 100.0, |t| {
            JoeCopula::new(t)
                .map(|c| c.pdf_log_sum(u, v))
                .unwrap_or(f64::NEG_INFINITY)
        });
        JoeCopula::new(theta_opt.max(1.0))
    }

    fn pdf_log_sum(&self, u: &[f64], v: &[f64]) -> f64 {
        u.iter()
            .zip(v.iter())
            .map(|(&ui, &vi)| {
                let p = self.pdf(ui, vi);
                if p > 0.0 { p.ln() } else { -1e15 }
            })
            .sum()
    }
}

// ---------------------------------------------------------------------------
// BB1 Copula (Clayton-Gumbel mixture)
// ---------------------------------------------------------------------------

/// BB1 copula: Clayton-Gumbel mixture with both lower and upper tail dependence.
///
/// Generator: φ(t) = (t^{-θ} - 1)^δ,  θ > 0, δ ≥ 1
/// C(u,v) = φ^{-1}(φ(u) + φ(v))
/// Lower tail: λ_L = 2^{-1/(θδ)}
/// Upper tail: λ_U = 2 - 2^{1/δ}
#[derive(Debug, Clone, PartialEq)]
pub struct BB1Copula {
    /// Clayton-like parameter θ > 0
    pub theta: f64,
    /// Gumbel-like parameter δ ≥ 1
    pub delta: f64,
}

impl BB1Copula {
    /// Create a BB1 copula. Requires θ > 0 and δ ≥ 1.
    pub fn new(theta: f64, delta: f64) -> StatsResult<Self> {
        if theta <= 0.0 {
            return Err(StatsError::InvalidArgument(
                "BB1 copula requires theta > 0".into(),
            ));
        }
        if delta < 1.0 {
            return Err(StatsError::InvalidArgument(
                "BB1 copula requires delta >= 1".into(),
            ));
        }
        Ok(Self { theta, delta })
    }

    /// Generator function φ(t) = (t^{-θ} - 1)^δ
    fn phi(&self, t: f64) -> f64 {
        if t <= 0.0 {
            return f64::INFINITY;
        }
        if t >= 1.0 {
            return 0.0;
        }
        (t.powf(-self.theta) - 1.0).max(0.0).powf(self.delta)
    }

    /// Inverse generator φ^{-1}(s) = (1 + s^{1/δ})^{-1/θ}
    fn phi_inv(&self, s: f64) -> f64 {
        if s <= 0.0 {
            return 1.0;
        }
        (1.0 + s.powf(1.0 / self.delta)).powf(-1.0 / self.theta)
    }

    /// Copula CDF: C(u,v) = φ^{-1}(φ(u) + φ(v))
    pub fn cdf(&self, u: f64, v: f64) -> f64 {
        if u <= 0.0 || v <= 0.0 {
            return 0.0;
        }
        if u >= 1.0 {
            return v.clamp(0.0, 1.0);
        }
        if v >= 1.0 {
            return u.clamp(0.0, 1.0);
        }
        self.phi_inv(self.phi(u) + self.phi(v)).clamp(0.0, 1.0)
    }

    /// Copula density via numerical differentiation.
    pub fn pdf(&self, u: f64, v: f64) -> f64 {
        if u <= 0.0 || v <= 0.0 || u >= 1.0 || v >= 1.0 {
            return 0.0;
        }
        // Use numerical second derivative for generality
        let h = 1e-4;
        let c_pp = self.cdf((u + h).min(1.0 - h), (v + h).min(1.0 - h));
        let c_pm = self.cdf((u + h).min(1.0 - h), (v - h).max(h));
        let c_mp = self.cdf((u - h).max(h), (v + h).min(1.0 - h));
        let c_mm = self.cdf((u - h).max(h), (v - h).max(h));
        let val = (c_pp - c_pm - c_mp + c_mm) / (4.0 * h * h);
        if val.is_finite() && val >= 0.0 { val } else { 0.0 }
    }

    /// Kendall's τ for BB1 (via numerical integration).
    pub fn kendall_tau(&self) -> f64 {
        // τ = 1 + 4 ∫₀¹ φ(t)/φ'(t) dt
        let n = 100usize;
        let h = 1.0 / (n as f64);
        let mut sum = 0.0;
        for i in 1..n {
            let t = i as f64 * h;
            let phi_t = self.phi(t);
            // φ'(t) = -δθ * t^{-θ-1} * (t^{-θ}-1)^{δ-1}
            let dphi = if t > 0.0 {
                let inner = t.powf(-self.theta) - 1.0;
                if inner <= 0.0 {
                    0.0
                } else {
                    -self.delta * self.theta * t.powf(-self.theta - 1.0) * inner.powf(self.delta - 1.0)
                }
            } else {
                0.0
            };
            if dphi.abs() > 1e-15 {
                sum += phi_t / dphi;
            }
        }
        (1.0 + 4.0 * sum * h).clamp(-1.0, 1.0)
    }

    /// Lower tail dependence: λ_L = 2^{-1/(θδ)}
    pub fn lower_tail_dependence(&self) -> f64 {
        2.0_f64.powf(-1.0 / (self.theta * self.delta))
    }

    /// Upper tail dependence: λ_U = 2 - 2^{1/δ}
    pub fn upper_tail_dependence(&self) -> f64 {
        2.0 - 2.0_f64.powf(1.0 / self.delta)
    }

    /// Sample pairs using conditional method.
    pub fn sample_pair(&self, n: usize, rng: &mut impl LcgRng) -> Vec<(f64, f64)> {
        let mut pairs = Vec::with_capacity(n);
        for _ in 0..n {
            let u = rng.next_unit();
            let w = rng.next_unit();
            let v = self.conditional_inverse(u, w);
            pairs.push((u, v));
        }
        pairs
    }

    /// Numerical inversion of conditional CDF via bisection.
    fn conditional_inverse(&self, u: f64, w: f64) -> f64 {
        // ∂C/∂u(u,v) via numerical derivative
        let h = 1e-5;
        let cond = |v: f64| -> f64 {
            let c_plus = self.cdf((u + h).min(1.0 - 1e-10), v);
            let c_minus = self.cdf((u - h).max(1e-10), v);
            (c_plus - c_minus) / (2.0 * h)
        };
        let mut lo = 1e-10_f64;
        let mut hi = 1.0 - 1e-10_f64;
        for _ in 0..50 {
            let mid = (lo + hi) / 2.0;
            if cond(mid) < w {
                hi = mid;
            } else {
                lo = mid;
            }
        }
        ((lo + hi) / 2.0).clamp(1e-15, 1.0 - 1e-15)
    }

    /// Fit BB1 copula via MLE (2D Nelder-Mead).
    pub fn fit(u: &[f64], v: &[f64]) -> StatsResult<BB1Copula> {
        if u.len() != v.len() || u.is_empty() {
            return Err(StatsError::InvalidArgument(
                "u and v must have the same positive length".into(),
            ));
        }

        let pdf_log_sum = |theta: f64, delta: f64| -> f64 {
            match BB1Copula::new(theta, delta) {
                Ok(c) => u
                    .iter()
                    .zip(v.iter())
                    .map(|(&ui, &vi)| {
                        let p = c.pdf(ui, vi);
                        if p > 0.0 { p.ln() } else { -1e15 }
                    })
                    .sum(),
                Err(_) => f64::NEG_INFINITY,
            }
        };

        // Grid search initialization
        let mut best_theta = 1.0;
        let mut best_delta = 1.5;
        let mut best_ll = f64::NEG_INFINITY;
        for &th in &[0.5, 1.0, 2.0, 3.0] {
            for &de in &[1.0, 1.5, 2.0, 3.0] {
                let ll = pdf_log_sum(th, de);
                if ll > best_ll {
                    best_ll = ll;
                    best_theta = th;
                    best_delta = de;
                }
            }
        }

        // Simple 2D gradient-free optimization
        let mut step_theta = 0.5;
        let mut step_delta = 0.3;
        for _ in 0..200 {
            let mut improved = false;
            for &(dt, dd) in &[
                (step_theta, 0.0_f64),
                (-step_theta, 0.0),
                (0.0, step_delta),
                (0.0, -step_delta),
            ] {
                let th = (best_theta + dt).max(1e-4);
                let de = (best_delta + dd).max(1.0);
                let ll = pdf_log_sum(th, de);
                if ll > best_ll {
                    best_ll = ll;
                    best_theta = th;
                    best_delta = de;
                    improved = true;
                }
            }
            if !improved {
                step_theta *= 0.5;
                step_delta *= 0.5;
                if step_theta < 1e-6 {
                    break;
                }
            }
        }
        BB1Copula::new(best_theta, best_delta)
    }
}

// ---------------------------------------------------------------------------
// Shared utilities
// ---------------------------------------------------------------------------

/// Simple LCG random number generator trait for sampling.
pub trait LcgRng {
    fn next_unit(&mut self) -> f64;
}

/// Built-in LCG RNG implementation.
pub struct SimpleLcg {
    state: u64,
}

impl SimpleLcg {
    pub fn new(seed: u64) -> Self {
        Self { state: seed.wrapping_add(1) }
    }
}

impl LcgRng for SimpleLcg {
    fn next_unit(&mut self) -> f64 {
        self.state = self.state
            .wrapping_mul(6_364_136_223_846_793_005)
            .wrapping_add(1_442_695_040_888_963_407);
        let high = (self.state >> 33) as f64;
        (high / (u32::MAX as f64 + 1.0)).clamp(1e-15, 1.0 - 1e-15)
    }
}

/// Compute Kendall's τ concordance measure.
pub fn compute_kendall_tau(u: &[f64], v: &[f64]) -> f64 {
    let n = u.len().min(v.len());
    if n < 2 {
        return 0.0;
    }
    let mut concordant = 0i64;
    let mut discordant = 0i64;
    for i in 0..n {
        for j in (i + 1)..n {
            let du = u[i] - u[j];
            let dv = v[i] - v[j];
            let prod = du * dv;
            if prod > 0.0 {
                concordant += 1;
            } else if prod < 0.0 {
                discordant += 1;
            }
        }
    }
    let total = (n * (n - 1) / 2) as f64;
    if total < 1.0 {
        return 0.0;
    }
    (concordant - discordant) as f64 / total
}

/// MLE optimization for single-parameter copulas using golden section search.
fn mle_1param(
    u: &[f64],
    v: &[f64],
    init: f64,
    lo_bound: f64,
    hi_bound: f64,
    ll_fn: impl Fn(f64) -> f64,
) -> f64 {
    let _ = (u, v); // used via ll_fn closure
    
    // Golden section search for maximum
    let golden = (5.0_f64.sqrt() - 1.0) / 2.0;
    let mut a = lo_bound.max(init * 0.1 - 1.0).max(lo_bound);
    let mut b = hi_bound.min(init * 2.0 + 2.0).min(hi_bound);
    
    // Clamp initial to valid range
    let init_c = init.clamp(a, b);
    
    // Expand search window around init_c
    a = (init_c - (b - a) * 0.5).max(lo_bound);
    b = (init_c + (b - a) * 0.5).min(hi_bound);
    
    // First evaluate at init and try gradient descent around it
    let mut best_x = init_c;
    let mut best_f = ll_fn(init_c);

    // Grid search first
    let n_grid = 20;
    let step = (b - a) / n_grid as f64;
    for i in 0..=n_grid {
        let x = a + i as f64 * step;
        let f = ll_fn(x);
        if f > best_f {
            best_f = f;
            best_x = x;
        }
    }

    // Refine with golden section around best_x
    let mut la = (best_x - step).max(lo_bound);
    let mut lb = (best_x + step).min(hi_bound);
    for _ in 0..50 {
        let x1 = lb - golden * (lb - la);
        let x2 = la + golden * (lb - la);
        if ll_fn(x1) < ll_fn(x2) {
            la = x1;
        } else {
            lb = x2;
        }
        if (lb - la).abs() < 1e-8 {
            break;
        }
    }
    best_x = (la + lb) / 2.0;
    best_x
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

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

    fn approx_eq(a: f64, b: f64, tol: f64) -> bool {
        (a - b).abs() < tol
    }

    // ---- Clayton ----

    #[test]
    fn test_clayton_new_invalid() {
        assert!(ClaytonCopula::new(0.0).is_err());
        assert!(ClaytonCopula::new(-1.0).is_err());
    }

    #[test]
    fn test_clayton_cdf_boundaries() {
        let c = ClaytonCopula::new(2.0).unwrap();
        assert_eq!(c.cdf(0.0, 0.5), 0.0);
        assert_eq!(c.cdf(0.5, 0.0), 0.0);
        assert!((c.cdf(1.0, 0.7) - 0.7).abs() < 1e-10);
        assert!((c.cdf(0.7, 1.0) - 0.7).abs() < 1e-10);
    }

    #[test]
    fn test_clayton_cdf_in_range() {
        let c = ClaytonCopula::new(2.0).unwrap();
        let val = c.cdf(0.5, 0.5);
        assert!(val >= 0.0 && val <= 0.5);
    }

    #[test]
    fn test_clayton_kendall_tau() {
        let c = ClaytonCopula::new(2.0).unwrap();
        let tau = c.kendall_tau();
        assert!(approx_eq(tau, 2.0 / 4.0, 1e-10));
    }

    #[test]
    fn test_clayton_tail_dependence() {
        let c = ClaytonCopula::new(2.0).unwrap();
        assert_eq!(c.upper_tail_dependence(), 0.0);
        let ltd = c.lower_tail_dependence();
        assert!(approx_eq(ltd, 2.0_f64.powf(-0.5), 1e-10));
    }

    #[test]
    fn test_clayton_sample() {
        let c = ClaytonCopula::new(2.0).unwrap();
        let mut rng = SimpleLcg::new(42);
        let pairs = c.sample_pair(100, &mut rng);
        assert_eq!(pairs.len(), 100);
        for (u, v) in &pairs {
            assert!(*u > 0.0 && *u < 1.0, "u={u}");
            assert!(*v > 0.0 && *v < 1.0, "v={v}");
        }
    }

    #[test]
    fn test_clayton_pdf_positive() {
        let c = ClaytonCopula::new(1.5).unwrap();
        let p = c.pdf(0.4, 0.6);
        assert!(p > 0.0, "pdf={p}");
    }

    #[test]
    fn test_clayton_fit() {
        let c = ClaytonCopula::new(2.0).unwrap();
        let mut rng = SimpleLcg::new(77);
        let pairs = c.sample_pair(200, &mut rng);
        let u: Vec<f64> = pairs.iter().map(|&(a, _)| a).collect();
        let v: Vec<f64> = pairs.iter().map(|&(_, b)| b).collect();
        let fitted = ClaytonCopula::fit(&u, &v).unwrap();
        assert!(fitted.theta > 0.0);
    }

    // ---- Gumbel ----

    #[test]
    fn test_gumbel_new_invalid() {
        assert!(GumbelCopula::new(0.5).is_err());
    }

    #[test]
    fn test_gumbel_kendall_tau() {
        let g = GumbelCopula::new(2.0).unwrap();
        assert!(approx_eq(g.kendall_tau(), 0.5, 1e-10));
    }

    #[test]
    fn test_gumbel_upper_tail() {
        let g = GumbelCopula::new(2.0).unwrap();
        let utd = g.upper_tail_dependence();
        // 2 - 2^{1/2} ≈ 0.586
        assert!(approx_eq(utd, 2.0 - 2.0_f64.sqrt(), 1e-10));
    }

    #[test]
    fn test_gumbel_independence_at_theta_1() {
        // θ=1 → independence copula
        let g = GumbelCopula::new(1.0).unwrap();
        // C(u,v) = u*v for independence
        let c = g.cdf(0.5, 0.5);
        // At θ=1: C(0.5, 0.5) = exp(-((-ln0.5)^1 + (-ln0.5)^1)^1) = exp(-2*ln2) = 0.25
        assert!(approx_eq(c, 0.25, 1e-8), "c={c}");
    }

    #[test]
    fn test_gumbel_fit() {
        let g = GumbelCopula::new(2.0).unwrap();
        let mut rng = SimpleLcg::new(99);
        let pairs = g.sample_pair(100, &mut rng);
        let u: Vec<f64> = pairs.iter().map(|&(a, _)| a).collect();
        let v: Vec<f64> = pairs.iter().map(|&(_, b)| b).collect();
        let fitted = GumbelCopula::fit(&u, &v).unwrap();
        assert!(fitted.theta >= 1.0);
    }

    // ---- Frank ----

    #[test]
    fn test_frank_new_invalid() {
        assert!(FrankCopula::new(0.0).is_err());
    }

    #[test]
    fn test_frank_cdf_at_half() {
        let f = FrankCopula::new(2.0).unwrap();
        let c = f.cdf(0.5, 0.5);
        assert!(c > 0.0 && c <= 0.5);
    }

    #[test]
    fn test_frank_pdf_positive() {
        let f = FrankCopula::new(3.0).unwrap();
        let p = f.pdf(0.4, 0.6);
        assert!(p > 0.0, "pdf={p}");
    }

    #[test]
    fn test_frank_no_tail_dependence() {
        let f = FrankCopula::new(5.0).unwrap();
        assert_eq!(f.upper_tail_dependence(), 0.0);
        assert_eq!(f.lower_tail_dependence(), 0.0);
    }

    #[test]
    fn test_frank_sample() {
        let f = FrankCopula::new(4.0).unwrap();
        let mut rng = SimpleLcg::new(13);
        let pairs = f.sample_pair(50, &mut rng);
        assert_eq!(pairs.len(), 50);
        for (u, v) in &pairs {
            assert!(*u > 0.0 && *u < 1.0 && *v > 0.0 && *v < 1.0);
        }
    }

    // ---- Joe ----

    #[test]
    fn test_joe_new_invalid() {
        assert!(JoeCopula::new(0.5).is_err());
    }

    #[test]
    fn test_joe_cdf_bounds() {
        let j = JoeCopula::new(2.0).unwrap();
        assert_eq!(j.cdf(0.0, 0.5), 0.0);
        assert!((j.cdf(1.0, 0.7) - 0.7).abs() < 1e-10);
    }

    #[test]
    fn test_joe_upper_tail() {
        let j = JoeCopula::new(2.0).unwrap();
        let utd = j.upper_tail_dependence();
        assert!(approx_eq(utd, 2.0 - 2.0_f64.sqrt(), 1e-10));
    }

    #[test]
    fn test_joe_sample() {
        let j = JoeCopula::new(3.0).unwrap();
        let mut rng = SimpleLcg::new(55);
        let pairs = j.sample_pair(50, &mut rng);
        assert_eq!(pairs.len(), 50);
        for (u, v) in &pairs {
            assert!(*u > 0.0 && *u < 1.0 && *v > 0.0 && *v < 1.0);
        }
    }

    // ---- BB1 ----

    #[test]
    fn test_bb1_new_invalid() {
        assert!(BB1Copula::new(0.0, 1.5).is_err());
        assert!(BB1Copula::new(1.0, 0.5).is_err());
    }

    #[test]
    fn test_bb1_cdf_bounds() {
        let b = BB1Copula::new(1.0, 1.0).unwrap();
        assert_eq!(b.cdf(0.0, 0.5), 0.0);
        assert!((b.cdf(1.0, 0.6) - 0.6).abs() < 1e-8);
    }

    #[test]
    fn test_bb1_tail_dependence() {
        let b = BB1Copula::new(1.0, 2.0).unwrap();
        let ltd = b.lower_tail_dependence();
        let utd = b.upper_tail_dependence();
        assert!(ltd > 0.0 && ltd < 1.0);
        assert!(utd > 0.0 && utd < 1.0);
    }

    #[test]
    fn test_bb1_sample() {
        let b = BB1Copula::new(1.0, 1.5).unwrap();
        let mut rng = SimpleLcg::new(777);
        let pairs = b.sample_pair(30, &mut rng);
        assert_eq!(pairs.len(), 30);
        for (u, v) in &pairs {
            assert!(*u > 0.0 && *u < 1.0 && *v > 0.0 && *v < 1.0);
        }
    }

    // ---- Utilities ----

    #[test]
    fn test_kendall_tau_perfect_concordance() {
        let u = vec![0.1, 0.2, 0.3, 0.4, 0.5];
        let v = vec![0.1, 0.2, 0.3, 0.4, 0.5];
        assert!(approx_eq(compute_kendall_tau(&u, &v), 1.0, 1e-10));
    }

    #[test]
    fn test_kendall_tau_perfect_discordance() {
        let u = vec![0.1, 0.2, 0.3, 0.4, 0.5];
        let v = vec![0.5, 0.4, 0.3, 0.2, 0.1];
        assert!(approx_eq(compute_kendall_tau(&u, &v), -1.0, 1e-10));
    }

    #[test]
    fn test_debye_function_at_zero() {
        // D₁(0) should approach 1
        assert!((debye_1(0.0) - 1.0).abs() < 1e-6);
    }
}