rs-stats 3.0.0

//! # Exponential Distribution
//!
//! This module implements the exponential distribution, a continuous probability distribution
//! that models the time between events in a Poisson point process.
//!
//! ## Key Characteristics
//! - Continuous probability distribution
//! - Memoryless property (future states depend only on the present, not the past)
//! - Used to model waiting times and inter-arrival times
//!
//! ## Common Applications
//! - Time between arrivals in a Poisson process
//! - Lifetime analysis (e.g., how long until a component fails)
//! - Queue theory and service times
//! - Radioactive decay
//!
//! ## Mathematical Formulation
//! The probability density function (PDF) is given by:
//!
//! f(x; λ) = λe^(-λx) for x ≥ 0
//!
//! where:
//! - λ (lambda) is the rate parameter (λ > 0)
//! - x is the random variable (time or space)
//!
//! The cumulative distribution function (CDF) is:
//!
//! F(x; λ) = 1 - e^(-λx) for x ≥ 0

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

/// Probability density function (PDF) for the Exponential distribution.
///
/// Calculates the probability density at point `x` for an exponential distribution
/// with rate parameter `lambda`.
///
/// # Arguments
/// * `x` - The point at which to evaluate the PDF (must be non-negative)
/// * `lambda` - The rate parameter (must be positive)
///
/// # Returns
/// The probability density at point `x`.
///
/// # Errors
/// Returns an error if:
/// - `x` is negative
/// - `lambda` is not positive
/// - Type conversion to f64 fails
///
#[inline]
fn exponential_pdf(x: f64, lambda: f64) -> StatsResult<f64> {
    if x < 0.0 {
        return Err(StatsError::InvalidInput {
            message: "exponential_pdf: x must be non-negative".to_string(),
        });
    }
    if lambda <= 0.0 {
        return Err(StatsError::InvalidInput {
            message: "exponential_pdf: lambda must be positive".to_string(),
        });
    }
    Ok(if x == 0.0 {
        lambda
    } else {
        lambda * (-lambda * x).exp()
    })
}

/// Cumulative distribution function (CDF) for the Exponential distribution.
///
/// Calculates the probability of a random variable being less than or equal to `x`
/// for an exponential distribution with rate parameter `lambda`.
///
/// # Arguments
/// * `x` - The point at which to evaluate the CDF (must be non-negative)
/// * `lambda` - The rate parameter (must be positive)
///
/// # Returns
/// The cumulative probability at point `x`.
///
/// # Errors
/// Returns an error if:
/// - `x` is negative
/// - `lambda` is not positive
/// - Type conversion to f64 fails
///
#[inline]
fn exponential_cdf(x: f64, lambda: f64) -> StatsResult<f64> {
    if x < 0.0 {
        return Err(StatsError::InvalidInput {
            message: "exponential_cdf: x must be non-negative".to_string(),
        });
    }
    if lambda <= 0.0 {
        return Err(StatsError::InvalidInput {
            message: "exponential_cdf: lambda must be positive".to_string(),
        });
    }
    Ok(1.0 - (-lambda * x).exp())
}

/// Inverse cumulative distribution function for the Exponential distribution.
///
/// Calculates the value of `x` for which the CDF equals the given probability `p`.
///
/// # Arguments
/// * `p` - The probability (must be between 0 and 1)
/// * `lambda` - The rate parameter (must be positive)
///
/// # Returns
/// The value `x` such that P(X ≤ x) = p.
///
/// # Errors
/// Returns an error if:
/// - `p` is not between 0 and 1
/// - `lambda` is not positive
/// - Type conversion to f64 fails
///
#[inline]
fn exponential_inverse_cdf(p: f64, lambda: f64) -> StatsResult<f64> {
    if !(0.0..=1.0).contains(&p) {
        return Err(StatsError::InvalidInput {
            message: "exponential_inverse_cdf: p must be between 0 and 1".to_string(),
        });
    }
    if lambda <= 0.0 {
        return Err(StatsError::InvalidInput {
            message: "exponential_inverse_cdf: lambda must be positive".to_string(),
        });
    }
    Ok(-((1.0 - p).ln()) / lambda)
}

/// Mean (expected value) of the Exponential distribution.
///
/// Calculates the mean of an exponential distribution with rate parameter `lambda`.
///
/// # Arguments
/// * `lambda` - The rate parameter (must be positive)
///
/// # Returns
/// The mean of the distribution.
///
/// # Errors
/// Returns an error if:
/// - `lambda` is not positive
/// - Type conversion to f64 fails
///

/// Variance of the Exponential distribution.
///
/// Calculates the variance of an exponential distribution with rate parameter `lambda`.
///
/// # Arguments
/// * `lambda` - The rate parameter (must be positive)
///
/// # Returns
/// The variance of the distribution.
///
/// # Errors
/// Returns an error if:
/// - `lambda` is not positive
/// - Type conversion to f64 fails
///

// ── Typed struct + Distribution impl ──────────────────────────────────────────

/// Exponential distribution Exp(λ) as a typed struct.
///
/// # Examples
/// ```
/// use rs_stats::distributions::exponential_distribution::Exponential;
/// use rs_stats::distributions::traits::Distribution;
///
/// let e = Exponential::new(2.0).unwrap();
/// assert!((e.mean() - 0.5).abs() < 1e-10);
/// ```
#[derive(Debug, Clone, Copy)]
pub struct Exponential {
    /// Rate parameter λ > 0
    pub lambda: f64,
}

impl Exponential {
    /// Creates an `Exponential` distribution with validation.
    pub fn new(lambda: f64) -> StatsResult<Self> {
        if lambda <= 0.0 {
            return Err(StatsError::InvalidInput {
                message: "Exponential::new: lambda must be positive".to_string(),
            });
        }
        Ok(Self { lambda })
    }

    /// MLE: λ = 1 / mean(data).  Data must be non-negative.
    pub fn fit(data: &[f64]) -> StatsResult<Self> {
        if data.is_empty() {
            return Err(StatsError::InvalidInput {
                message: "Exponential::fit: data must not be empty".to_string(),
            });
        }
        if data.iter().any(|&x| x < 0.0) {
            return Err(StatsError::InvalidInput {
                message: "Exponential::fit: all data values must be non-negative".to_string(),
            });
        }
        let mean = data.iter().sum::<f64>() / data.len() as f64;
        Self::new(1.0 / mean)
    }
}

impl crate::distributions::traits::Distribution for Exponential {
    type X = f64;
    fn name(&self) -> &str {
        "Exponential"
    }
    fn num_params(&self) -> usize {
        1
    }
    fn pdf(&self, x: f64) -> StatsResult<f64> {
        exponential_pdf(x, self.lambda)
    }
    fn logpdf(&self, x: f64) -> StatsResult<f64> {
        if x < 0.0 {
            return Err(StatsError::InvalidInput {
                message: "Exponential::logpdf: x must be non-negative".to_string(),
            });
        }
        Ok(self.lambda.ln() - self.lambda * x)
    }
    fn cdf(&self, x: f64) -> StatsResult<f64> {
        exponential_cdf(x, self.lambda)
    }
    fn inverse_cdf(&self, p: f64) -> StatsResult<f64> {
        exponential_inverse_cdf(p, self.lambda)
    }
    fn mean(&self) -> f64 {
        1.0 / self.lambda
    }
    fn variance(&self) -> f64 {
        1.0 / (self.lambda * self.lambda)
    }
}

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

    const EPSILON: f64 = 1e-10;

    #[test]
    fn test_exponential_pdf() {
        let lambda = 2.0;

        // PDF at x = 0
        let result = exponential_pdf(0.0, lambda).unwrap();
        assert_eq!(result, lambda);

        // PDF at x = 1
        let result = exponential_pdf(1.0, lambda).unwrap();
        let expected = lambda * (-lambda).exp();
        assert!((result - expected).abs() < EPSILON);

        // PDF at x = 0.5
        let result = exponential_pdf(0.5, lambda).unwrap();
        let expected = lambda * (-lambda * 0.5).exp();
        assert!((result - expected).abs() < EPSILON);
    }

    #[test]
    fn test_exponential_cdf() {
        let lambda = 2.0_f64;

        // CDF at x = 0
        let result = exponential_cdf(0.0, lambda).unwrap();
        assert!((result - 0.0).abs() < EPSILON);

        // CDF at x = 1
        let result = exponential_cdf(1.0, lambda).unwrap();
        let expected = 1.0 - (-lambda).exp();
        assert!((result - expected).abs() < EPSILON);

        // CDF at x = 0.5
        let result = exponential_cdf(0.5, lambda).unwrap();
        let expected = 1.0 - (-lambda * 0.5).exp();
        assert!((result - expected).abs() < EPSILON);
    }

    #[test]
    fn test_exponential_inverse_cdf() {
        let lambda = 2.0_f64;

        // Test inverse CDF with various probabilities
        let test_cases = vec![0.1, 0.25, 0.5, 0.75, 0.9];

        for p in test_cases {
            let x = exponential_inverse_cdf(p, lambda).unwrap();
            let cdf = exponential_cdf(x, lambda).unwrap();
            assert!(
                (cdf - p).abs() < EPSILON,
                "Inverse CDF failed for p = {}: got {}, expected {}",
                p,
                cdf,
                p
            );
        }
    }

    #[test]
    fn test_exponential_mean_via_trait() {
        use crate::distributions::traits::Distribution;
        let dist = Exponential::new(2.0).unwrap();
        assert!((dist.mean() - 0.5).abs() < EPSILON);
        assert!((dist.variance() - 0.25).abs() < EPSILON);
    }

    #[test]
    fn test_exponential_pdf_invalid_lambda() {
        let result = exponential_pdf(1.0, -2.0);
        assert!(result.is_err());
        match result {
            Err(StatsError::InvalidInput { message }) => {
                assert!(message.contains("lambda must be positive"));
            }
            _ => panic!("Expected InvalidInput error"),
        }
    }

    #[test]
    fn test_exponential_pdf_invalid_x() {
        let result = exponential_pdf(-1.0, 2.0);
        assert!(result.is_err());
        match result {
            Err(StatsError::InvalidInput { message }) => {
                assert!(message.contains("x must be non-negative"));
            }
            _ => panic!("Expected InvalidInput error"),
        }
    }

    #[test]
    fn test_exponential_inverse_cdf_p_negative() {
        let result = exponential_inverse_cdf(-0.1, 2.0);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            StatsError::InvalidInput { .. }
        ));
    }

    #[test]
    fn test_exponential_inverse_cdf_p_greater_than_one() {
        let result = exponential_inverse_cdf(1.5, 2.0);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            StatsError::InvalidInput { .. }
        ));
    }

    #[test]
    fn test_exponential_cdf_invalid_lambda() {
        let result = exponential_cdf(1.0, -2.0);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            StatsError::InvalidInput { .. }
        ));
    }

    #[test]
    fn test_exponential_cdf_invalid_x() {
        let result = exponential_cdf(-1.0, 2.0);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            StatsError::InvalidInput { .. }
        ));
    }

    #[test]
    fn test_exponential_pdf_x_positive() {
        // Test the branch where x > 0 (not x == 0)
        let result = exponential_pdf(0.5, 2.0).unwrap();
        let lambda: f64 = 2.0;
        let x: f64 = 0.5;
        let expected = lambda * (-lambda * x).exp();
        assert!((result - expected).abs() < EPSILON);
    }

    #[test]
    fn test_exponential_inverse_cdf_p_zero() {
        // When p = 0, inverse CDF should be 0 (or very close to 0)
        let result = exponential_inverse_cdf(0.0, 2.0).unwrap();
        assert_eq!(result, 0.0);
    }

    #[test]
    fn test_exponential_inverse_cdf_p_one() {
        // When p = 1, inverse CDF should approach infinity
        // But due to numerical limits, we check it's very large
        let result = exponential_inverse_cdf(1.0, 2.0).unwrap();
        assert!(result.is_infinite() || result > 1e10);
    }

    #[test]
    fn test_exponential_inverse_cdf_lambda_zero() {
        let result = exponential_inverse_cdf(0.5, 0.0);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            StatsError::InvalidInput { .. }
        ));
    }

    #[test]
    fn test_exponential_inverse_cdf_lambda_negative() {
        let result = exponential_inverse_cdf(0.5, -1.0);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            StatsError::InvalidInput { .. }
        ));
    }
}