u_analytics/weibull/

reliability.rs

//! Reliability analysis from fitted Weibull parameters.
//!
//! Provides reliability function, hazard rate, MTBF, B-life, and other
//! common reliability engineering metrics.

use u_numflow::special::gamma;

use super::mle::{two_sided_z, WeibullInterval};

/// Reliability analysis results from a fitted Weibull distribution.
///
/// Computes reliability engineering metrics (reliability function, hazard
/// rate, MTBF, B-life) from Weibull shape (beta) and scale (eta) parameters.
///
/// # Mathematical Background
///
/// Given a Weibull distribution with shape beta > 0 and scale eta > 0:
/// - Reliability: R(t) = exp(-(t/eta)^beta)
/// - Hazard rate: lambda(t) = (beta/eta) * (t/eta)^(beta-1)
/// - MTBF: eta * Gamma(1 + 1/beta)
///
/// # Examples
///
/// ```
/// use u_analytics::weibull::ReliabilityAnalysis;
/// let ra = ReliabilityAnalysis::new(2.0, 100.0).unwrap();
/// assert!((ra.reliability(0.0) - 1.0).abs() < 1e-10);
/// assert!(ra.hazard_rate(50.0) > 0.0);
/// assert!(ra.mtbf() > 0.0);
/// let b10 = ra.b_life(0.10).unwrap();
/// assert!(b10 > 0.0 && b10 < 100.0);
/// ```
///
/// # References
/// - Meeker & Escobar (1998), *Statistical Methods for Reliability Data*, Wiley.
/// - Nelson, W. (1982). *Applied Life Data Analysis*, Wiley. §2.
/// - Lawless, J.F. (2003). *Statistical Models and Methods for Lifetime Data*, 2nd ed.
#[derive(Debug, Clone)]
pub struct ReliabilityAnalysis {
    /// Shape parameter (beta).
    shape: f64,
    /// Scale parameter (eta).
    scale: f64,
    /// Covariance of `(shape, scale)`, when the parameters were estimated --
    /// what the confidence intervals are made from.
    covariance: Option<[[f64; 2]; 2]>,
}

impl ReliabilityAnalysis {
    /// Creates a new reliability analysis from Weibull parameters.
    ///
    /// # Arguments
    /// * `shape` - Shape parameter beta (must be positive and finite)
    /// * `scale` - Scale parameter eta (must be positive and finite)
    ///
    /// # Returns
    /// `None` if either parameter is non-positive or non-finite.
    ///
    /// # Examples
    ///
    /// ```
    /// use u_analytics::weibull::ReliabilityAnalysis;
    /// let ra = ReliabilityAnalysis::new(2.0, 100.0);
    /// assert!(ra.is_some());
    ///
    /// // Invalid parameters return None
    /// assert!(ReliabilityAnalysis::new(-1.0, 100.0).is_none());
    /// assert!(ReliabilityAnalysis::new(2.0, 0.0).is_none());
    /// ```
    pub fn new(shape: f64, scale: f64) -> Option<Self> {
        if !shape.is_finite() || !scale.is_finite() || shape <= 0.0 || scale <= 0.0 {
            return None;
        }
        Some(Self {
            shape,
            scale,
            covariance: None,
        })
    }

    /// The same analysis with the covariance of the estimated `(shape,
    /// scale)` -- row and column order shape, scale, as
    /// [`WeibullMleResult::covariance`](super::WeibullMleResult::covariance)
    /// gives it -- so the B-lives and reliabilities carry confidence
    /// intervals ([`Self::b_life_interval`], [`Self::reliability_interval`]).
    ///
    /// `None` unless the matrix is a covariance: finite, symmetric (to a
    /// relative 1e-9) and positive definite.
    ///
    /// # Examples
    ///
    /// ```
    /// use u_analytics::weibull::ReliabilityAnalysis;
    /// let ra = ReliabilityAnalysis::new(2.0, 100.0).unwrap();
    /// assert!(ra.clone().with_covariance([[0.04, 0.1], [0.1, 25.0]]).is_some());
    /// assert!(ra.with_covariance([[0.04, 2.0], [2.0, 25.0]]).is_none()); // not positive definite
    /// ```
    pub fn with_covariance(self, covariance: [[f64; 2]; 2]) -> Option<Self> {
        let [[a, b], [c, d]] = covariance;
        let finite = [a, b, c, d].iter().all(|v| v.is_finite());
        let symmetric = (b - c).abs() <= 1e-9 * b.abs().max(c.abs()).max(f64::MIN_POSITIVE);
        if !(finite && symmetric && a > 0.0 && d > 0.0 && a * d - b * c > 0.0) {
            return None;
        }
        Some(Self {
            covariance: Some(covariance),
            ..self
        })
    }

    /// The covariance of `(shape, scale)` the intervals are made from, when
    /// there is one.
    pub fn covariance(&self) -> Option<[[f64; 2]; 2]> {
        self.covariance
    }

    /// Creates a reliability analysis from an MLE fitting result.
    ///
    /// # Panics
    /// This function does not panic. MLE results always have valid parameters.
    ///
    /// The fit's covariance comes along, so the analysis gives confidence
    /// intervals.
    pub fn from_mle(result: &super::mle::WeibullMleResult) -> Self {
        Self {
            shape: result.shape,
            scale: result.scale,
            covariance: Some(result.covariance),
        }
    }

    /// Creates a reliability analysis from an MRR fitting result.
    ///
    /// # Panics
    /// This function does not panic. MRR results always have valid parameters.
    pub fn from_mrr(result: &super::mrr::WeibullMrrResult) -> Self {
        Self {
            shape: result.shape,
            scale: result.scale,
            covariance: None,
        }
    }

    /// Returns the shape parameter (beta).
    pub fn shape(&self) -> f64 {
        self.shape
    }

    /// Returns the scale parameter (eta).
    pub fn scale(&self) -> f64 {
        self.scale
    }

    /// Reliability (survival) function at time t.
    ///
    /// ```text
    /// R(t) = exp(-(t/eta)^beta)
    /// ```
    ///
    /// For t < 0, returns 1.0 (no failure before time 0).
    ///
    /// # Examples
    ///
    /// ```
    /// use u_analytics::weibull::ReliabilityAnalysis;
    /// let ra = ReliabilityAnalysis::new(2.0, 100.0).unwrap();
    ///
    /// // R(0) = 1.0 (full reliability at time zero)
    /// assert!((ra.reliability(0.0) - 1.0).abs() < 1e-10);
    ///
    /// // R(eta) = exp(-1) ≈ 0.368 for any shape parameter
    /// let expected = (-1.0_f64).exp();
    /// assert!((ra.reliability(100.0) - expected).abs() < 1e-10);
    /// ```
    ///
    /// # Reference
    /// Weibull (1951), *Journal of Applied Mechanics* 18(3), pp. 293-297.
    pub fn reliability(&self, t: f64) -> f64 {
        if t <= 0.0 {
            return 1.0;
        }
        let z = t / self.scale;
        (-z.powf(self.shape)).exp()
    }

    /// Failure rate (hazard function) at time t.
    ///
    /// ```text
    /// lambda(t) = (beta/eta) * (t/eta)^(beta-1)
    /// ```
    ///
    /// - beta < 1: Decreasing failure rate (infant mortality)
    /// - beta = 1: Constant failure rate (random/exponential failures)
    /// - beta > 1: Increasing failure rate (wear-out)
    ///
    /// For t <= 0, returns 0.0.
    ///
    /// # Examples
    ///
    /// ```
    /// use u_analytics::weibull::ReliabilityAnalysis;
    /// let ra = ReliabilityAnalysis::new(2.0, 100.0).unwrap();
    ///
    /// // Hazard rate is positive for t > 0
    /// assert!(ra.hazard_rate(50.0) > 0.0);
    ///
    /// // With beta > 1, hazard rate increases over time (wear-out)
    /// assert!(ra.hazard_rate(50.0) < ra.hazard_rate(80.0));
    /// ```
    ///
    /// # Reference
    /// Meeker & Escobar (1998), *Statistical Methods for Reliability Data*, Ch. 4.
    pub fn hazard_rate(&self, t: f64) -> f64 {
        if t <= 0.0 {
            return 0.0;
        }
        let z = t / self.scale;
        (self.shape / self.scale) * z.powf(self.shape - 1.0)
    }

    /// Mean Time Between Failures (MTBF).
    ///
    /// ```text
    /// MTBF = eta * Gamma(1 + 1/beta)
    /// ```
    ///
    /// # Examples
    ///
    /// ```
    /// use u_analytics::weibull::ReliabilityAnalysis;
    /// let ra = ReliabilityAnalysis::new(2.0, 100.0).unwrap();
    /// let mtbf = ra.mtbf();
    /// assert!(mtbf > 0.0);
    ///
    /// // For beta=1 (exponential), MTBF = eta
    /// let exp_ra = ReliabilityAnalysis::new(1.0, 50.0).unwrap();
    /// assert!((exp_ra.mtbf() - 50.0).abs() < 1e-8);
    /// ```
    ///
    /// # Reference
    /// Johnson, Kotz & Balakrishnan (1994), *Continuous Univariate Distributions*,
    /// Vol. 1, Chapter 21.
    pub fn mtbf(&self) -> f64 {
        self.scale * gamma(1.0 + 1.0 / self.shape)
    }

    /// Time at which reliability drops to a given level.
    ///
    /// Solves R(t) = p for t:
    ///
    /// ```text
    /// t = eta * (-ln(p))^(1/beta)
    /// ```
    ///
    /// # Arguments
    /// * `p` - Desired reliability level, must be in (0, 1).
    ///
    /// # Returns
    /// `None` if `p` is outside (0, 1).
    ///
    /// # Reference
    /// Abernethy (2006), *The New Weibull Handbook*, 5th ed.
    pub fn time_to_reliability(&self, p: f64) -> Option<f64> {
        if p <= 0.0 || p >= 1.0 {
            return None;
        }
        Some(self.scale * (-p.ln()).powf(1.0 / self.shape))
    }

    /// B-life: time at which a given fraction of the population has failed.
    ///
    /// B10 life (10% failed) = `b_life(0.10)`, which is equivalent to
    /// `time_to_reliability(0.90)`.
    ///
    /// # Arguments
    /// * `fraction_failed` - Fraction of population that has failed, must be in (0, 1).
    ///
    /// # Returns
    /// `None` if `fraction_failed` is outside (0, 1).
    ///
    /// # Examples
    ///
    /// ```
    /// use u_analytics::weibull::ReliabilityAnalysis;
    /// let ra = ReliabilityAnalysis::new(2.0, 100.0).unwrap();
    ///
    /// // B10: time when 10% of the population has failed
    /// let b10 = ra.b_life(0.10).unwrap();
    /// assert!(b10 > 0.0 && b10 < 100.0);
    ///
    /// // B-lives increase with failure fraction: B5 < B10 < B50
    /// let b5 = ra.b_life(0.05).unwrap();
    /// let b50 = ra.b_life(0.50).unwrap();
    /// assert!(b5 < b10 && b10 < b50);
    /// ```
    ///
    /// # Reference
    /// Abernethy (2006), *The New Weibull Handbook*, 5th ed., Chapter 2.
    pub fn b_life(&self, fraction_failed: f64) -> Option<f64> {
        if fraction_failed <= 0.0 || fraction_failed >= 1.0 {
            return None;
        }
        self.time_to_reliability(1.0 - fraction_failed)
    }

    /// `g' Σ g` for a gradient in `(shape, scale)`.
    fn variance_along(&self, g: [f64; 2]) -> Option<f64> {
        let c = self.covariance?;
        Some((g[0] * g[0] * c[0][0] + 2.0 * g[0] * g[1] * c[0][1] + g[1] * g[1] * c[1][1]).max(0.0))
    }

    /// The B-life at `fraction_failed` with its two-sided interval at
    /// `confidence`, by the delta method on
    /// `ln t_p = ln η + ln(−ln(1 − p)) / β` (normal on the log scale, so the
    /// limits stay positive).
    ///
    /// Returns `(estimate, interval)`; `None` without a covariance
    /// ([`Self::with_covariance`], [`Self::from_mle`]) or when either argument
    /// is outside (0, 1).
    ///
    /// # Examples
    ///
    /// ```
    /// use u_analytics::weibull::ReliabilityAnalysis;
    /// let ra = ReliabilityAnalysis::new(2.0, 100.0)
    ///     .unwrap()
    ///     .with_covariance([[0.04, 0.1], [0.1, 25.0]])
    ///     .unwrap();
    /// let (b10, ci) = ra.b_life_interval(0.10, 0.90).unwrap();
    /// assert!(ci.lower < b10 && b10 < ci.upper);
    /// ```
    ///
    /// # Reference
    /// Meeker & Escobar (1998), *Statistical Methods for Reliability Data*, §8.4.
    pub fn b_life_interval(
        &self,
        fraction_failed: f64,
        confidence: f64,
    ) -> Option<(f64, WeibullInterval)> {
        let t_p = self.b_life(fraction_failed)?;
        let z = two_sided_z(confidence)?;
        let (beta, eta) = (self.shape, self.scale);
        let w = (-(1.0 - fraction_failed).ln()).ln();
        // Gradient of ln t_p in (beta, eta).
        let half = z * self.variance_along([-w / (beta * beta), 1.0 / eta])?.sqrt();
        Some((
            t_p,
            WeibullInterval {
                lower: t_p * (-half).exp(),
                upper: t_p * half.exp(),
            },
        ))
    }

    /// The reliability `R(t)` with its two-sided interval at `confidence`, by
    /// the delta method on `u = ln(−ln R(t)) = β (ln t − ln η)`, which is
    /// unbounded -- so the limits `exp(−exp(u ± z·se))` stay inside (0, 1).
    /// At `t ≤ 0` the model has no failures and the interval is `[1, 1]`.
    ///
    /// Returns `(estimate, interval)`; `None` without a covariance or when
    /// `confidence` is outside (0, 1) or `t` is not finite.
    ///
    /// # Examples
    ///
    /// ```
    /// use u_analytics::weibull::ReliabilityAnalysis;
    /// let ra = ReliabilityAnalysis::new(2.0, 100.0)
    ///     .unwrap()
    ///     .with_covariance([[0.04, 0.1], [0.1, 25.0]])
    ///     .unwrap();
    /// let (r, ci) = ra.reliability_interval(50.0, 0.95).unwrap();
    /// assert!(0.0 < ci.lower && ci.lower < r && r < ci.upper && ci.upper < 1.0);
    /// ```
    ///
    /// # Reference
    /// Meeker & Escobar (1998), *Statistical Methods for Reliability Data*, §8.4.
    pub fn reliability_interval(&self, t: f64, confidence: f64) -> Option<(f64, WeibullInterval)> {
        let z = two_sided_z(confidence)?;
        if !t.is_finite() {
            return None;
        }
        self.covariance?;
        if t <= 0.0 {
            return Some((
                1.0,
                WeibullInterval {
                    lower: 1.0,
                    upper: 1.0,
                },
            ));
        }
        let (beta, eta) = (self.shape, self.scale);
        let u = beta * (t / eta).ln();
        // Gradient of u in (beta, eta).
        let half = z * self.variance_along([u / beta, -beta / eta])?.sqrt();
        let r = |u: f64| (-u.exp()).exp();
        Some((
            self.reliability(t),
            WeibullInterval {
                lower: r(u + half),
                upper: r(u - half),
            },
        ))
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::weibull::{weibull_mle, weibull_mrr};

    fn estimated() -> ReliabilityAnalysis {
        ReliabilityAnalysis::new(1.8, 120.0)
            .unwrap()
            .with_covariance([[0.05, 0.9], [0.9, 60.0]])
            .unwrap()
    }

    /// The delta-method variance with the gradient taken by central
    /// differences of the function itself, not the closed form.
    fn numeric_variance(ra: &ReliabilityAnalysis, f: impl Fn(f64, f64) -> f64) -> f64 {
        let (b, e) = (ra.shape(), ra.scale());
        let (hb, he) = (1e-6 * b, 1e-6 * e);
        let g = [
            (f(b + hb, e) - f(b - hb, e)) / (2.0 * hb),
            (f(b, e + he) - f(b, e - he)) / (2.0 * he),
        ];
        let c = ra.covariance().unwrap();
        g[0] * g[0] * c[0][0] + 2.0 * g[0] * g[1] * c[0][1] + g[1] * g[1] * c[1][1]
    }

    #[test]
    fn interval_widths_are_the_delta_method_of_the_function() {
        let ra = estimated();
        let z = u_numflow::special::inverse_normal_cdf(0.975);
        let p = 0.1;
        let (t_p, ci) = ra.b_life_interval(p, 0.95).unwrap();
        let ln_tp = |b: f64, e: f64| e.ln() + (-(1.0 - p).ln()).ln() / b;
        let half = z * numeric_variance(&ra, ln_tp).sqrt();
        assert!(((ci.upper / t_p).ln() - half).abs() < 1e-7, "{ci:?}");
        assert!(((t_p / ci.lower).ln() - half).abs() < 1e-7, "{ci:?}");

        let t = 75.0;
        let (r, ci) = ra.reliability_interval(t, 0.95).unwrap();
        let u = |b: f64, e: f64| b * (t / e).ln();
        let half = z * numeric_variance(&ra, u).sqrt();
        let u_of = |r: f64| (-r.ln()).ln();
        assert!((u_of(ci.lower) - u_of(r) - half).abs() < 1e-7, "{ci:?}");
        assert!((u_of(r) - u_of(ci.upper) - half).abs() < 1e-7, "{ci:?}");
    }

    /// At t = B_p the two intervals are one statement: u = β(ln t − ln η) is
    /// −β times (ln t_p − ln η − w/β) moved, so se(u) = β · se(ln t_p).
    #[test]
    fn the_reliability_interval_at_a_b_life_is_that_b_life_s_interval() {
        let ra = estimated();
        let p = 0.2;
        let (t_p, b) = ra.b_life_interval(p, 0.9).unwrap();
        let (r, ci) = ra.reliability_interval(t_p, 0.9).unwrap();
        assert!((r - (1.0 - p)).abs() < 1e-12);
        let half_u = (-ci.lower.ln()).ln() - (-r.ln()).ln();
        let half_t = (b.upper / t_p).ln();
        assert!(
            (half_u - ra.shape() * half_t).abs() < 1e-9,
            "{half_u} vs {half_t}"
        );
    }

    #[test]
    fn reliability_limits_stay_inside_zero_one_and_are_exact_at_time_zero() {
        let ra = estimated();
        for t in [1e-3, 1.0, 50.0, 120.0, 400.0, 2000.0] {
            let (r, ci) = ra.reliability_interval(t, 0.99).unwrap();
            assert!(
                0.0 <= ci.lower && ci.lower <= r && r <= ci.upper && ci.upper <= 1.0,
                "{t}: {ci:?}"
            );
        }
        let (r, ci) = ra.reliability_interval(0.0, 0.95).unwrap();
        assert_eq!((r, ci.lower, ci.upper), (1.0, 1.0, 1.0));
        assert!(ra.reliability_interval(f64::NAN, 0.95).is_none());
        assert!(ra.reliability_interval(10.0, 1.0).is_none());
    }

    #[test]
    fn intervals_need_a_covariance() {
        let ra = ReliabilityAnalysis::new(2.0, 100.0).unwrap();
        assert!(ra.b_life_interval(0.1, 0.95).is_none());
        assert!(ra.reliability_interval(10.0, 0.95).is_none());
        let mrr = weibull_mrr(&[12.0, 25.0, 33.0, 47.0, 58.0, 71.0]).unwrap();
        assert!(ReliabilityAnalysis::from_mrr(&mrr).covariance().is_none());
    }

    #[test]
    fn with_covariance_takes_only_a_covariance() {
        let ra = || ReliabilityAnalysis::new(2.0, 100.0).unwrap();
        assert!(ra().with_covariance([[0.04, 0.1], [0.1, 25.0]]).is_some());
        assert!(
            ra().with_covariance([[0.04, 0.1], [0.2, 25.0]]).is_none(),
            "asymmetric"
        );
        assert!(
            ra().with_covariance([[0.04, 1.0], [1.0, 25.0]]).is_none(),
            "singular"
        );
        assert!(
            ra().with_covariance([[-0.04, 0.0], [0.0, 25.0]]).is_none(),
            "negative variance"
        );
        assert!(ra()
            .with_covariance([[f64::NAN, 0.0], [0.0, 25.0]])
            .is_none());
    }

    #[test]
    fn the_mle_s_b_life_interval_is_this_one() {
        let fit = weibull_mle(&[55.0, 187.0, 216.0, 240.0, 244.0, 335.0, 361.0, 373.0]).unwrap();
        let ra = ReliabilityAnalysis::from_mle(&fit);
        assert_eq!(fit.b_life_interval(0.1, 0.9), ra.b_life_interval(0.1, 0.9));
    }

    /// Monte Carlo: n = 40 Weibull(2, 100) samples, MLE, 90 % interval for
    /// R(60) and for B10. A Wald interval at this n covers within a few
    /// points of nominal (Meeker & Escobar §8.4 — the log/loglog scales are
    /// chosen for that); 600 samples give a binomial s.e. of 1.2 %.
    #[test]
    fn the_intervals_cover_near_their_level() {
        use u_numflow::distributions::{Sample, Weibull};
        let truth = ReliabilityAnalysis::new(2.0, 100.0).unwrap();
        let (r_true, b10_true) = (truth.reliability(60.0), truth.b_life(0.1).unwrap());
        let law = Weibull::new(2.0, 100.0).unwrap();
        let mut rng = u_numflow::random::create_rng(2026);
        let runs = 600;
        let (mut r_in, mut b_in) = (0, 0);
        for _ in 0..runs {
            let data: Vec<f64> = (0..40).map(|_| law.sample(&mut rng)).collect();
            let ra = ReliabilityAnalysis::from_mle(&weibull_mle(&data).unwrap());
            let (_, r) = ra.reliability_interval(60.0, 0.9).unwrap();
            let (_, b) = ra.b_life_interval(0.1, 0.9).unwrap();
            r_in += usize::from(r.lower <= r_true && r_true <= r.upper);
            b_in += usize::from(b.lower <= b10_true && b10_true <= b.upper);
        }
        let (rc, bc) = (r_in as f64 / runs as f64, b_in as f64 / runs as f64);
        assert!((0.85..=0.95).contains(&rc), "R(60) coverage {rc}");
        assert!((0.85..=0.95).contains(&bc), "B10 coverage {bc}");
    }

    #[test]
    fn test_new_valid() {
        assert!(ReliabilityAnalysis::new(2.0, 50.0).is_some());
    }

    #[test]
    fn test_new_invalid() {
        assert!(ReliabilityAnalysis::new(0.0, 50.0).is_none());
        assert!(ReliabilityAnalysis::new(-1.0, 50.0).is_none());
        assert!(ReliabilityAnalysis::new(2.0, 0.0).is_none());
        assert!(ReliabilityAnalysis::new(2.0, -1.0).is_none());
        assert!(ReliabilityAnalysis::new(f64::NAN, 50.0).is_none());
        assert!(ReliabilityAnalysis::new(2.0, f64::INFINITY).is_none());
    }

    #[test]
    fn test_reliability_at_zero() {
        let ra = ReliabilityAnalysis::new(2.0, 50.0).expect("valid parameters");
        assert!((ra.reliability(0.0) - 1.0).abs() < 1e-15);
    }

    #[test]
    fn test_reliability_decreasing() {
        let ra = ReliabilityAnalysis::new(2.0, 50.0).expect("valid parameters");
        let mut prev = 1.0;
        for i in 1..=100 {
            let t = i as f64;
            let r = ra.reliability(t);
            assert!(
                r <= prev + 1e-15,
                "reliability should be non-increasing: R({}) = {} > R({}) = {}",
                t,
                r,
                t - 1.0,
                prev
            );
            prev = r;
        }
    }

    #[test]
    fn test_reliability_at_scale() {
        // R(eta) = exp(-1) for any beta
        let ra = ReliabilityAnalysis::new(2.0, 50.0).expect("valid parameters");
        let r_at_scale = ra.reliability(50.0);
        let expected = (-1.0_f64).exp();
        assert!(
            (r_at_scale - expected).abs() < 1e-10,
            "R(eta) = {}, expected exp(-1) = {}",
            r_at_scale,
            expected
        );
    }

    #[test]
    fn test_reliability_negative_time() {
        let ra = ReliabilityAnalysis::new(2.0, 50.0).expect("valid parameters");
        assert!((ra.reliability(-10.0) - 1.0).abs() < 1e-15);
    }

    #[test]
    fn test_hazard_rate_exponential() {
        // beta=1 => constant hazard rate = 1/eta
        let ra = ReliabilityAnalysis::new(1.0, 20.0).expect("valid parameters");
        for t in [5.0, 10.0, 20.0, 50.0] {
            assert!(
                (ra.hazard_rate(t) - 1.0 / 20.0).abs() < 1e-10,
                "hazard at t={} should be 1/eta=0.05",
                t
            );
        }
    }

    #[test]
    fn test_hazard_rate_increasing() {
        // beta > 1 => increasing hazard rate (wear-out)
        let ra = ReliabilityAnalysis::new(3.0, 50.0).expect("valid parameters");
        let h1 = ra.hazard_rate(10.0);
        let h2 = ra.hazard_rate(20.0);
        let h3 = ra.hazard_rate(30.0);
        assert!(
            h1 < h2,
            "hazard should increase: h(10)={} >= h(20)={}",
            h1,
            h2
        );
        assert!(
            h2 < h3,
            "hazard should increase: h(20)={} >= h(30)={}",
            h2,
            h3
        );
    }

    #[test]
    fn test_hazard_rate_zero_time() {
        let ra = ReliabilityAnalysis::new(2.0, 50.0).expect("valid parameters");
        assert!((ra.hazard_rate(0.0)).abs() < 1e-15);
        assert!((ra.hazard_rate(-5.0)).abs() < 1e-15);
    }

    #[test]
    fn test_mtbf_exponential() {
        // beta=1 => MTBF = eta * Gamma(2) = eta * 1 = eta
        let ra = ReliabilityAnalysis::new(1.0, 100.0).expect("valid parameters");
        assert!(
            (ra.mtbf() - 100.0).abs() < 1e-8,
            "MTBF = {}, expected 100.0 for exponential",
            ra.mtbf()
        );
    }

    #[test]
    fn test_mtbf_rayleigh() {
        // beta=2, eta=1 => MTBF = Gamma(1.5) = sqrt(pi)/2
        let ra = ReliabilityAnalysis::new(2.0, 1.0).expect("valid parameters");
        let expected = std::f64::consts::PI.sqrt() / 2.0;
        assert!(
            (ra.mtbf() - expected).abs() < 1e-10,
            "MTBF = {}, expected sqrt(pi)/2 = {}",
            ra.mtbf(),
            expected
        );
    }

    #[test]
    fn test_time_to_reliability_median() {
        // Median life: R(t) = 0.5 => t = eta * (ln(2))^(1/beta)
        let ra = ReliabilityAnalysis::new(2.0, 50.0).expect("valid parameters");
        let t_median = ra.time_to_reliability(0.5).expect("p=0.5 is valid");
        let expected = 50.0 * (2.0_f64.ln()).powf(0.5);
        assert!(
            (t_median - expected).abs() < 1e-10,
            "median life = {}, expected {}",
            t_median,
            expected
        );
    }

    #[test]
    fn test_time_to_reliability_invalid() {
        let ra = ReliabilityAnalysis::new(2.0, 50.0).expect("valid parameters");
        assert!(ra.time_to_reliability(0.0).is_none());
        assert!(ra.time_to_reliability(1.0).is_none());
        assert!(ra.time_to_reliability(-0.1).is_none());
        assert!(ra.time_to_reliability(1.5).is_none());
    }

    #[test]
    fn test_time_to_reliability_roundtrip() {
        let ra = ReliabilityAnalysis::new(2.5, 100.0).expect("valid parameters");
        for p in [0.1, 0.25, 0.5, 0.75, 0.9] {
            let t = ra.time_to_reliability(p).expect("valid p");
            let r = ra.reliability(t);
            assert!(
                (r - p).abs() < 1e-10,
                "roundtrip: time_to_reliability({}) = {}, reliability({}) = {}",
                p,
                t,
                t,
                r
            );
        }
    }

    #[test]
    fn test_b_life_b10() {
        // B10 = time when 10% have failed = time_to_reliability(0.90)
        let ra = ReliabilityAnalysis::new(2.0, 50.0).expect("valid parameters");
        let b10 = ra.b_life(0.10).expect("valid fraction");
        let t_r90 = ra.time_to_reliability(0.90).expect("valid p");
        assert!(
            (b10 - t_r90).abs() < 1e-10,
            "B10 = {}, time_to_reliability(0.90) = {}",
            b10,
            t_r90
        );
    }

    #[test]
    fn test_b_life_invalid() {
        let ra = ReliabilityAnalysis::new(2.0, 50.0).expect("valid parameters");
        assert!(ra.b_life(0.0).is_none());
        assert!(ra.b_life(1.0).is_none());
        assert!(ra.b_life(-0.1).is_none());
    }

    #[test]
    fn test_from_mle() {
        let data = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0];
        let mle = weibull_mle(&data).expect("MLE should converge");
        let ra = ReliabilityAnalysis::from_mle(&mle);

        assert!((ra.shape() - mle.shape).abs() < 1e-15);
        assert!((ra.scale() - mle.scale).abs() < 1e-15);
        assert!(ra.mtbf() > 0.0);
    }

    #[test]
    fn test_from_mrr() {
        let data = [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0];
        let mrr = weibull_mrr(&data).expect("MRR should succeed");
        let ra = ReliabilityAnalysis::from_mrr(&mrr);

        assert!((ra.shape() - mrr.shape).abs() < 1e-15);
        assert!((ra.scale() - mrr.scale).abs() < 1e-15);
        assert!(ra.mtbf() > 0.0);
    }

    #[test]
    fn test_b_life_ordering() {
        // B5 < B10 < B50 (more failures => more time)
        let ra = ReliabilityAnalysis::new(2.0, 50.0).expect("valid parameters");
        let b5 = ra.b_life(0.05).expect("valid");
        let b10 = ra.b_life(0.10).expect("valid");
        let b50 = ra.b_life(0.50).expect("valid");
        assert!(b5 < b10, "B5={} should < B10={}", b5, b10);
        assert!(b10 < b50, "B10={} should < B50={}", b10, b50);
    }

    #[test]
    fn test_hazard_pdf_reliability_relation() {
        // h(t) = f(t) / R(t) where f(t) is the Weibull PDF
        // f(t) = (beta/eta) * (t/eta)^(beta-1) * exp(-(t/eta)^beta)
        let ra = ReliabilityAnalysis::new(2.5, 50.0).expect("valid parameters");
        for t in [5.0, 20.0, 50.0, 80.0] {
            let z = t / ra.scale();
            let pdf =
                (ra.shape() / ra.scale()) * z.powf(ra.shape() - 1.0) * (-z.powf(ra.shape())).exp();
            let h = ra.hazard_rate(t);
            let r = ra.reliability(t);
            assert!(
                (h * r - pdf).abs() < 1e-10,
                "h(t)*R(t) = {} should equal f(t) = {} at t={}",
                h * r,
                pdf,
                t
            );
        }
    }

    #[test]
    fn test_accessors() {
        let ra = ReliabilityAnalysis::new(2.5, 100.0).expect("valid parameters");
        assert!((ra.shape() - 2.5).abs() < 1e-15);
        assert!((ra.scale() - 100.0).abs() < 1e-15);
    }

    /// Reference value tests: beta=2, eta=100.
    ///
    /// Expected values derived analytically:
    /// - Nelson (1982) §2; Lawless (2003) Ch. 2.
    #[test]
    fn test_reference_reliability_beta2_eta100() {
        let ra = ReliabilityAnalysis::new(2.0, 100.0).expect("valid parameters");

        // R(0) = 1.0 exactly
        assert!(
            (ra.reliability(0.0) - 1.0).abs() < 1e-15,
            "R(0) = {}, expected 1.0",
            ra.reliability(0.0)
        );

        // R(100) = exp(-1) ≈ 0.367879
        let r100 = ra.reliability(100.0);
        let expected_r100 = (-1.0_f64).exp(); // 0.36787944117144233
        assert!(
            (r100 - expected_r100).abs() < 1e-5,
            "R(100) = {:.6}, expected {:.6}",
            r100,
            expected_r100
        );

        // R(50) = exp(-0.25) ≈ 0.778801
        let r50 = ra.reliability(50.0);
        let expected_r50 = (-0.25_f64).exp(); // 0.7788007830714049
        assert!(
            (r50 - expected_r50).abs() < 1e-5,
            "R(50) = {:.6}, expected {:.6}",
            r50,
            expected_r50
        );
    }

    #[test]
    fn test_reference_cdf_beta2_eta100() {
        let ra = ReliabilityAnalysis::new(2.0, 100.0).expect("valid parameters");

        // F(100) = 1 - R(100) = 1 - exp(-1) ≈ 0.632121
        let f100 = 1.0 - ra.reliability(100.0);
        let expected_f100 = 1.0 - (-1.0_f64).exp(); // 0.6321205588285578
        assert!(
            (f100 - expected_f100).abs() < 1e-5,
            "F(100) = {:.6}, expected {:.6}",
            f100,
            expected_f100
        );
    }

    #[test]
    fn test_reference_hazard_rate_beta2_eta100() {
        let ra = ReliabilityAnalysis::new(2.0, 100.0).expect("valid parameters");

        // h(100) = (2/100)*(100/100)^1 = 0.02
        let h100 = ra.hazard_rate(100.0);
        assert!(
            (h100 - 0.02).abs() < 1e-6,
            "h(100) = {:.8}, expected 0.02",
            h100
        );

        // h(50) = (2/100)*(50/100)^1 = 0.01
        let h50 = ra.hazard_rate(50.0);
        assert!(
            (h50 - 0.01).abs() < 1e-6,
            "h(50) = {:.8}, expected 0.01",
            h50
        );
    }

    #[test]
    fn test_reference_mttf_beta2_eta100() {
        let ra = ReliabilityAnalysis::new(2.0, 100.0).expect("valid parameters");

        // MTTF = 100 * Gamma(1.5) = 100 * sqrt(pi)/2 ≈ 88.6227
        let mttf = ra.mtbf();
        let expected_mttf = 100.0 * std::f64::consts::PI.sqrt() / 2.0; // 88.62269254527580...
        assert!(
            (mttf - expected_mttf).abs() < 0.001,
            "MTTF = {:.6}, expected {:.6}",
            mttf,
            expected_mttf
        );
    }

    #[test]
    fn test_reference_b10_beta2_eta100() {
        let ra = ReliabilityAnalysis::new(2.0, 100.0).expect("valid parameters");

        // B10: F(t)=0.10 => t = 100*(-ln(0.9))^0.5 ≈ 32.458
        let b10 = ra.b_life(0.10).expect("valid fraction");
        let expected_b10 = 100.0 * (-0.9_f64.ln()).powf(0.5); // 100 * sqrt(0.10536) ≈ 32.4576
        assert!(
            (b10 - expected_b10).abs() < 0.01,
            "B10 = {:.5}, expected {:.5}",
            b10,
            expected_b10
        );
    }
}