commit 749184191d72900f7ea49400b1b3b2c48b5e04d8 Author: Stephen Waits Date: Tue Aug 27 12:48:10 2024 -0600 feat: initial commit diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..96ef6c0 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +/target +Cargo.lock diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..6f5802b --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,23 @@ +[package] +name = "distimate" +version = "0.1.0" +edition = "2021" +authors = ["Stephen Waits "] +description = "A crate providing probability distributions for estimation and risk analysis" +license = "MIT" +repository = "https://github.com/swaits/distimate" +documentation = "https://docs.rs/distimate" +readme = "README.md" +keywords = ["probability", "estimation", "risk-analysis", "monte-carlo", "pert"] +categories = ["mathematics", "science", "simulation"] + +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html + +[dependencies] +rand = "0.8.5" +statrs = "0.17.1" +thiserror = "1.0.63" + +[dev-dependencies] +approx = "0.5.1" +rand = { version = "0.8.5", features = ["std_rng"] } diff --git a/README.md b/README.md new file mode 100644 index 0000000..fc3a726 --- /dev/null +++ b/README.md @@ -0,0 +1,94 @@ +# distimate + +`distimate` is a Rust crate that provides probability distributions specifically +designed for estimation and risk analysis scenarios. It offers implementations +of various distributions that are commonly used in project management, cost +estimation, and other fields where uncertainty needs to be modeled. + +## Features + +- Implementations of PERT, Triangular, and Normal distributions tailored for + estimation purposes +- Consistent API across all distributions +- Support for common statistical operations (mean, variance, skewness, etc.) +- Sampling capabilities for Monte Carlo simulations +- Range-respecting implementations that work with min, likely (mode), and max estimates + +## Installation + +Add this to your `Cargo.toml`: + +```toml +[dependencies] +distimate = "" +``` + +## Usage + +Here's a quick example of how to use the PERT distribution: + +```rust +use distimate::prelude::*; +use distimate::Pert; + +fn main() { + let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + println!("Mean: {}", pert.mean().unwrap()); + println!("Variance: {}", pert.variance().unwrap()); + println!("95th percentile: {}", pert.inverse_cdf(0.95)); +} +``` + +## Distributions + +### PERT Distribution + +The PERT (Program Evaluation and Review Technique) distribution is commonly used +in project management and risk analysis. It's defined by minimum, most likely +(mode), and maximum values. + +```rust +let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); +``` + +### Triangular Distribution + +The Triangular distribution is often used when limited sample data is available. +It's also defined by minimum, most likely (mode), and maximum values. + +```rust +let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); +``` + +### Normal Distribution + +The Normal (Gaussian) distribution is implemented with modifications to respect +a given range. It's suitable for natural phenomena and when values are expected +to be symmetrically distributed around a mean. + +```rust +let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); +``` + +## Common Operations + +All distributions implement common traits for statistical operations: + +- `Distribution`: Provides `mean()`, `variance()`, `skewness()`, and `entropy()` +- `Median`: Provides `median()` +- `Mode`: Provides `mode()` +- `Continuous`: Provides `pdf()` and `ln_pdf()` +- `ContinuousCDF`: Provides `cdf()` and `inverse_cdf()` +- `Min` and `Max`: Provide `min()` and `max()` + +Additionally, all distributions can be sampled using the `sample()` method, +which is useful for Monte Carlo simulations. + +## Contributing + +Contributions are welcome! Please feel free to submit a Pull Request. + +## License + +This project is licensed under the MIT License - see the [LICENSE](LICENSE) file +for details. diff --git a/TODO.md b/TODO.md new file mode 100644 index 0000000..87fba22 --- /dev/null +++ b/TODO.md @@ -0,0 +1,52 @@ +Distributions to include: + +1. PERT (Program Evaluation and Review Technique) Distribution + + - Useful when you have estimates for minimum, most likely, and maximum + values, with a tendency to weight the most likely value more heavily. + +2. Triangular Distribution + + - Ideal when you have estimates for minimum, most likely, and maximum values, + but less certainty about the weight of the most likely value compared to + PERT. + +3. Normal (Gaussian) Distribution + + - Appropriate for natural phenomena and when you expect values to be + symmetrically distributed around a mean. + +4. Log-Normal Distribution + + - Suitable for modeling quantities that are always positive and have a skewed + distribution, such as project durations or costs. + +5. Uniform Distribution + + - Used when all values between a minimum and maximum are equally likely, + often in the absence of more specific information. + +6. Beta Distribution + + - Flexible for modeling various shapes of uncertainty, especially useful when + historical data is available to fit the distribution. + +7. Gamma Distribution + + - Often used for modeling waiting times or project durations, especially when + the distribution is skewed to the right. + +8. Weibull Distribution + + - Commonly used in reliability analysis and for modeling time-to-failure in + projects or systems. + +9. Exponential Distribution + + - Useful for modeling the time between independent events, such as time + between failures in a system. + +10. Discrete Uniform Distribution + + - Applicable when dealing with a finite number of equally likely outcomes, + such as selecting from a set of discrete options. diff --git a/src/distributions/mod.rs b/src/distributions/mod.rs new file mode 100644 index 0000000..5f37343 --- /dev/null +++ b/src/distributions/mod.rs @@ -0,0 +1,7 @@ +mod normal; +mod pert; +mod triangular; + +pub use self::normal::Normal; +pub use self::pert::Pert; +pub use self::triangular::Triangular; diff --git a/src/distributions/normal.rs b/src/distributions/normal.rs new file mode 100644 index 0000000..41f1cf8 --- /dev/null +++ b/src/distributions/normal.rs @@ -0,0 +1,742 @@ +//! Normal (Gaussian) distribution implementation for estimation. +//! +//! This module provides an implementation of a modified Normal distribution, +//! which is defined by minimum, most likely (mode), and maximum values. +//! It's designed for use in estimation scenarios where a range and most likely +//! value are known, but the underlying distribution is assumed to be normal. +//! +//! # Examples +//! +//! ``` +//! use distimate::prelude::*; +//! use distimate::Normal; +//! use approx::assert_relative_eq; +//! +//! let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); +//! assert_eq!(normal.mode(), 2.0); +//! assert_relative_eq!(normal.mean().unwrap(), 2.0, epsilon = 1e-6); +//! ``` + +use crate::prelude::*; +use statrs::distribution::Normal as StatrsNormal; + +/// Represents a modified Normal distribution for estimation purposes. +/// +/// This Normal distribution is defined by minimum, most likely (mode), and maximum values. +/// It uses these to create an underlying standard Normal distribution, but provides +/// methods to work within the specified range. +/// +/// # Examples +/// +/// ``` +/// use distimate::prelude::*; +/// use distimate::Normal; +/// use approx::assert_relative_eq; +/// +/// let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); +/// assert_eq!(normal.min(), 1.0); +/// assert_eq!(normal.mode(), 2.0); +/// assert_eq!(normal.max(), 3.0); +/// ``` +#[derive(Debug, Clone)] +pub struct Normal { + inner: StatrsNormal, + min: f64, + likely: f64, + max: f64, +} + +impl Normal { + /// Creates a new Normal distribution with the given minimum, most likely (mode), and maximum values. + /// + /// # Arguments + /// + /// * `min` - The minimum value of the distribution + /// * `likely` - The most likely value (mode) of the distribution + /// * `max` - The maximum value of the distribution + /// + /// # Returns + /// + /// Returns a `Result` containing the new `Normal` instance if the parameters are valid, + /// or an `Error` if the parameters are invalid. + /// + /// # Errors + /// + /// Returns an error if: + /// - `min` is greater than or equal to `max` + /// - `likely` is less than `min` or greater than `max` + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Normal; + /// + /// let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + /// assert_eq!(normal.min(), 1.0); + /// assert_eq!(normal.mode(), 2.0); + /// assert_eq!(normal.max(), 3.0); + /// ``` + pub fn new(min: f64, likely: f64, max: f64) -> Result { + if min >= max { + return Err(Error::InvalidRange { min, max }); + } + if likely < min || likely > max { + return Err(Error::InvalidMode { min, likely, max }); + } + + let mean = likely; + let std_dev = (max - min) / 6.0; // Assuming 99.7% of values fall within the range + + Ok(Normal { + inner: StatrsNormal::new(mean, std_dev)?, + min, + likely, + max, + }) + } + + /// Clamps a value to the [min, max] range of the distribution. + fn clamp(&self, x: f64) -> f64 { + x.max(self.min).min(self.max) + } +} + +/// Implementation of the `EstimationDistribution` trait for the Normal distribution. +/// +/// This trait marks the Normal distribution as suitable for use in estimation contexts. +/// It doesn't add any new methods, but signals that this distribution is appropriate +/// for modeling uncertain estimates or measurements within a specified range. +/// +/// The Normal distribution in this context is particularly useful when: +/// - There's a most likely value (mode) with symmetric uncertainty around it. +/// - The minimum and maximum values represent extreme scenarios. +/// - The underlying process is believed to follow a bell-shaped curve. +/// +/// # Examples +/// +/// ``` +/// use distimate::prelude::*; +/// use distimate::Normal; +/// +/// let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); +/// // The Normal distribution can now be used in contexts requiring an EstimationDistribution +/// ``` +impl EstimationDistribution for Normal {} + +/// Implementation of the `Distribution` trait for the Normal distribution. +/// +/// This implementation provides methods to calculate various statistical properties +/// of the Normal distribution, adapted for estimation scenarios with specified +/// minimum, most likely (mode), and maximum values. +/// +/// # Note +/// +/// The calculated values are based on the underlying Normal distribution and may +/// theoretically extend beyond the specified [min, max] range. For practical +/// estimation purposes, you may want to consider clamping extreme values. +impl Distribution for Normal { + /// Calculates the mean of the Normal distribution. + /// + /// For this adapted Normal distribution, the mean is equal to the most likely value (mode). + /// + /// # Returns + /// + /// Returns `Some(mean)` where `mean` is the most likely value specified during construction. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Normal; + /// + /// let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + /// assert_eq!(normal.mean(), Some(2.0)); + /// ``` + fn mean(&self) -> Option { + self.inner.mean() + } + + /// Calculates the variance of the Normal distribution. + /// + /// The variance is calculated based on the standard deviation derived from + /// the specified min and max values during construction. + /// + /// # Returns + /// + /// Returns `Some(variance)` where `variance` is the square of the standard deviation. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Normal; + /// + /// let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + /// assert_eq!(normal.variance().unwrap(), 0.1111111111111111); // (1/3)^2 + /// ``` + fn variance(&self) -> Option { + self.inner.variance() + } + + /// Calculates the skewness of the Normal distribution. + /// + /// The Normal distribution is always symmetric, so the skewness is always 0. + /// + /// # Returns + /// + /// Always returns `Some(0.0)`. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Normal; + /// + /// let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + /// assert_eq!(normal.skewness(), Some(0.0)); + /// ``` + fn skewness(&self) -> Option { + Some(0.0) // Normal distribution is always symmetric + } + + /// Calculates the entropy of the Normal distribution. + /// + /// The entropy is a measure of the average amount of information contained in + /// the distribution. For a Normal distribution, it depends on the standard deviation. + /// + /// # Returns + /// + /// Returns `Some(entropy)` where `entropy` is calculated based on the + /// standard deviation of the distribution. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Normal; + /// + /// let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + /// assert!(normal.entropy().is_some()); + /// ``` + fn entropy(&self) -> Option { + self.inner.entropy() + } +} + +/// Implementation of the `Median` trait for the Normal distribution. +/// +/// For a Normal distribution adapted for estimation purposes, the median is equal +/// to the most likely value (mode) specified during construction. This is because +/// the Normal distribution is symmetric around its mean, which in this case is set +/// to the most likely value. +impl Median for Normal { + /// Calculates the median of the Normal distribution. + /// + /// For this adapted Normal distribution, the median is equal to the most likely + /// value (mode) specified during construction. + /// + /// # Returns + /// + /// Returns the median of the distribution as an f64 value. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Normal; + /// + /// let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + /// assert_eq!(normal.median(), 2.0); + /// ``` + /// + /// # Note + /// + /// In a standard Normal distribution, the median, mean, and mode are all equal. + /// In this adapted version for estimation, we maintain this property by setting + /// all of these to the most likely value provided during construction. + fn median(&self) -> f64 { + self.likely + } +} + +/// Implementation of the `Mode` trait for the Normal distribution. +/// +/// For a Normal distribution adapted for estimation purposes, the mode is set to +/// the most likely value specified during construction. This aligns with the +/// typical use case in estimation scenarios where the most likely value is known. +impl Mode for Normal { + /// Returns the mode of the Normal distribution. + /// + /// The mode is the most likely value of the distribution, which in this case + /// is explicitly set during the construction of the Normal distribution instance. + /// + /// # Returns + /// + /// Returns the mode of the distribution as an f64 value, which is equal to + /// the 'likely' value provided when creating the distribution. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Normal; + /// + /// let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + /// assert_eq!(normal.mode(), 2.0); + /// ``` + /// + /// # Note + /// + /// In a standard Normal distribution, the mode, median, and mean are all equal. + /// In this adapted version for estimation, we maintain this property by setting + /// all of these to the most likely value provided during construction. + fn mode(&self) -> f64 { + self.likely + } +} + +/// Implementation of the `Continuous` trait for the Normal distribution. +/// +/// This implementation provides methods to calculate the probability density function (PDF) +/// and its natural logarithm for the Normal distribution, adapted for estimation scenarios +/// with specified minimum and maximum values. +impl Continuous for Normal { + /// Calculates the probability density function (PDF) of the Normal distribution at a given point. + /// + /// The PDF represents the relative likelihood of the distribution taking on a specific value. + /// For this adapted Normal distribution, the PDF is set to 0 for values outside the [min, max] range. + /// + /// # Arguments + /// + /// * `x` - The point at which to calculate the PDF + /// + /// # Returns + /// + /// The value of the PDF at the given point as an f64. Returns 0 if `x` is outside the [min, max] range. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Normal; + /// + /// let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + /// assert!(normal.pdf(2.0) > 0.0); + /// assert_eq!(normal.pdf(0.0), 0.0); // Outside the distribution's range + /// ``` + fn pdf(&self, x: f64) -> f64 { + if x < self.min || x > self.max { + 0.0 + } else { + self.inner.pdf(x) + } + } + + /// Calculates the natural logarithm of the probability density function (PDF) of the Normal distribution at a given point. + /// + /// This method is particularly useful for numerical stability in calculations involving very small probability densities. + /// For this adapted Normal distribution, the ln_pdf is set to negative infinity for values outside the [min, max] range. + /// + /// # Arguments + /// + /// * `x` - The point at which to calculate the log PDF + /// + /// # Returns + /// + /// The natural logarithm of the PDF at the given point as an f64. + /// Returns negative infinity if `x` is outside the [min, max] range. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Normal; + /// + /// let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + /// assert!(normal.ln_pdf(2.0).is_finite()); + /// assert_eq!(normal.ln_pdf(0.0), f64::NEG_INFINITY); // Outside the distribution's range + /// ``` + fn ln_pdf(&self, x: f64) -> f64 { + if x < self.min || x > self.max { + f64::NEG_INFINITY + } else { + self.inner.ln_pdf(x) + } + } +} + +/// Implementation of the `ContinuousCDF` trait for the Normal distribution. +/// +/// This implementation provides methods to calculate the cumulative distribution function (CDF) +/// and its inverse for the Normal distribution, adapted for estimation scenarios with +/// specified minimum and maximum values. +impl ContinuousCDF for Normal { + /// Calculates the cumulative distribution function (CDF) of the Normal distribution at a given point. + /// + /// The CDF represents the probability that a random variable from this distribution + /// will be less than or equal to the given value. For this adapted Normal distribution, + /// the CDF is clamped to 0 for values below the minimum and to 1 for values above the maximum. + /// + /// # Arguments + /// + /// * `x` - The point at which to calculate the CDF + /// + /// # Returns + /// + /// The value of the CDF at the given point as an f64, in the range [0, 1]. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Normal; + /// use approx::assert_relative_eq; + /// + /// let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(normal.cdf(2.0), 0.5, epsilon = 1e-6); + /// assert_eq!(normal.cdf(0.0), 0.0); // Below minimum + /// assert_eq!(normal.cdf(4.0), 1.0); // Above maximum + /// ``` + fn cdf(&self, x: f64) -> f64 { + if x <= self.min { + 0.0 + } else if x >= self.max { + 1.0 + } else { + self.inner.cdf(x) + } + } + + /// Calculates the inverse of the cumulative distribution function (inverse CDF) + /// of the Normal distribution for a given probability. + /// + /// This function is also known as the quantile function. It returns the value x + /// for which P(X ≤ x) = p, where p is the given probability. For this adapted + /// Normal distribution, the result is clamped to the [min, max] range. + /// + /// # Arguments + /// + /// * `p` - The probability for which to calculate the inverse CDF, must be in the range [0, 1] + /// + /// # Returns + /// + /// The value x for which the CDF of the distribution equals the given probability, + /// clamped to the [min, max] range of the distribution. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Normal; + /// use approx::assert_relative_eq; + /// + /// let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(normal.inverse_cdf(0.5), 2.0, epsilon = 1e-6); + /// assert_eq!(normal.inverse_cdf(0.0), 1.0); // Clamped to minimum + /// assert_eq!(normal.inverse_cdf(1.0), 3.0); // Clamped to maximum + /// ``` + fn inverse_cdf(&self, p: f64) -> f64 { + self.clamp(self.inner.inverse_cdf(p)) + } +} + +/// Implementation of the `Min` trait for the Normal distribution. +/// +/// This implementation provides a method to retrieve the minimum value +/// of the Normal distribution adapted for estimation scenarios. +impl Min for Normal { + /// Returns the minimum value of the Normal distribution. + /// + /// In the context of this estimation-focused Normal distribution, + /// the minimum value represents the lower bound of the possible outcomes. + /// It's one of the three key parameters (along with the most likely value + /// and the maximum) used to define this adapted Normal distribution. + /// + /// # Returns + /// + /// The minimum value of the distribution as an f64. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Normal; + /// + /// let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + /// assert_eq!(normal.min(), 1.0); + /// ``` + /// + /// # Note + /// + /// While a standard Normal distribution extends infinitely in both directions, + /// this adapted version for estimation purposes has a defined minimum value. + /// This allows for more practical modeling in scenarios where there's a known + /// lower limit to possible outcomes. + fn min(&self) -> f64 { + self.min + } +} + +/// Implementation of the `Max` trait for the Normal distribution. +/// +/// This implementation provides a method to retrieve the maximum value +/// of the Normal distribution adapted for estimation scenarios. +impl Max for Normal { + /// Returns the maximum value of the Normal distribution. + /// + /// In the context of this estimation-focused Normal distribution, + /// the maximum value represents the upper bound of the possible outcomes. + /// It's one of the three key parameters (along with the minimum value + /// and the most likely value) used to define this adapted Normal distribution. + /// + /// # Returns + /// + /// The maximum value of the distribution as an f64. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Normal; + /// + /// let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + /// assert_eq!(normal.max(), 3.0); + /// ``` + /// + /// # Note + /// + /// While a standard Normal distribution extends infinitely in both directions, + /// this adapted version for estimation purposes has a defined maximum value. + /// This allows for more practical modeling in scenarios where there's a known + /// upper limit to possible outcomes. + fn max(&self) -> f64 { + self.max + } +} + +/// Implementation of the `rand::distributions::Distribution` trait for the Normal distribution. +/// +/// This implementation allows for random sampling from the Normal distribution +/// adapted for estimation scenarios, using any random number generator that +/// implements the `rand::Rng` trait. +impl RandDistribution for Normal { + /// Generates a random sample from the Normal distribution. + /// + /// This method uses the underlying Normal distribution to generate a sample, + /// then clamps the result to ensure it falls within the specified [min, max] range. + /// This approach maintains the shape of the Normal distribution within the + /// specified range while ensuring all samples are valid for the estimation scenario. + /// + /// # Arguments + /// + /// * `rng` - A mutable reference to a random number generator that implements the `rand::Rng` trait + /// + /// # Returns + /// + /// A random sample from the Normal distribution as an f64, guaranteed to be + /// within the [min, max] range specified during the distribution's construction. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Normal; + /// use rand::prelude::*; + /// + /// let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + /// let mut rng = StdRng::seed_from_u64(42); // For reproducibility + /// + /// let sample = normal.sample(&mut rng); + /// assert!(sample >= 1.0 && sample <= 3.0); + /// ``` + /// + /// # Note + /// + /// While this method ensures all samples fall within the specified range, + /// it may lead to a slight overrepresentation of the minimum and maximum values + /// compared to a true Normal distribution. This trade-off is often acceptable + /// in practical estimation scenarios where staying within the specified range is crucial. + fn sample(&self, rng: &mut R) -> f64 { + self.clamp(self.inner.sample(rng)) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + use rand::distributions::Distribution as RandDistribution; + use rand::rngs::StdRng; + use rand::SeedableRng; + use statrs::statistics::{Data, Distribution, OrderStatistics}; + + #[test] + fn test_new() { + let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + assert_eq!(normal.min(), 1.0); + assert_eq!(normal.mode(), 2.0); + assert_eq!(normal.max(), 3.0); + assert_eq!(normal.mean().unwrap(), 2.0); + assert_relative_eq!(normal.std_dev().unwrap(), 1.0 / 3.0, epsilon = 1e-6); + } + + #[test] + fn test_invalid_parameters() { + assert!(Normal::new(3.0, 2.0, 1.0).is_err()); // min > max + assert!(Normal::new(1.0, 0.5, 3.0).is_err()); // likely < min + assert!(Normal::new(1.0, 3.5, 3.0).is_err()); // likely > max + } + + #[test] + fn test_pdf() { + let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + assert_relative_eq!(normal.pdf(2.0), normal.inner.pdf(2.0), epsilon = 1e-6); + assert_eq!(normal.pdf(0.5), 0.0); // Below min + assert_eq!(normal.pdf(3.5), 0.0); // Above max + } + + #[test] + fn test_cdf() { + let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + assert_relative_eq!(normal.cdf(2.0), 0.5, epsilon = 1e-6); + assert_eq!(normal.cdf(0.5), 0.0); // Below min + assert_eq!(normal.cdf(3.5), 1.0); // Above max + } + + #[test] + fn test_inverse_cdf() { + let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + assert_relative_eq!(normal.inverse_cdf(0.5), 2.0, epsilon = 1e-6); + assert_eq!(normal.inverse_cdf(0.0), 1.0); // Clamped to min + assert_eq!(normal.inverse_cdf(1.0), 3.0); // Clamped to max + } + + #[test] + fn test_sampling() { + let normal = Normal::new(1.0, 2.0, 3.0).unwrap(); + let mut rng = StdRng::seed_from_u64(42); + let samples: Vec = (0..1000).map(|_| normal.sample(&mut rng)).collect(); + + assert!(samples.iter().all(|&x| (1.0..=3.0).contains(&x))); + let mean = samples.iter().sum::() / samples.len() as f64; + assert_relative_eq!(mean, 2.0, epsilon = 0.1); + } + + #[test] + fn test_distribution_statistics() { + let test_cases = vec![ + (1.0, 2.0, 3.0), // Standard case + (0.0, 5.0, 10.0), // Wider range + (-5.0, 0.0, 5.0), // Range including negative values + (100.0, 200.0, 300.0), // Larger values + ]; + + for (min, likely, max) in test_cases { + let normal = Normal::new(min, likely, max).unwrap(); + let n = 100_000; // Number of samples + let mut rng = StdRng::seed_from_u64(42); // Fixed seed for reproducibility + + let samples: Vec = (0..n).map(|_| normal.sample(&mut rng)).collect(); + let mut data = Data::new(samples.clone()); + + // Mean + let sample_mean = data.mean().unwrap(); + let theoretical_mean = normal.mean().unwrap(); + assert_relative_eq!( + sample_mean, + theoretical_mean, + epsilon = 0.01, + max_relative = 0.01 + ); + + // Variance + let sample_variance = data.variance().unwrap(); + let theoretical_variance = normal.variance().unwrap(); + assert_relative_eq!( + sample_variance, + theoretical_variance, + epsilon = 0.05, + max_relative = 0.05 + ); + + // Standard deviation + let sample_std_dev = data.std_dev().unwrap(); + let theoretical_std_dev = normal.std_dev().unwrap(); + assert_relative_eq!( + sample_std_dev, + theoretical_std_dev, + epsilon = 0.05, + max_relative = 0.05 + ); + + // Skewness + let sample_skewness = 0.0; // normal distribution has no skewness + let theoretical_skewness = normal.skewness().unwrap(); + assert_relative_eq!( + sample_skewness, + theoretical_skewness, + epsilon = 0.1, + max_relative = 0.1 + ); + + // Median + let sample_median = data.median(); + let theoretical_median = normal.median(); + assert_relative_eq!( + sample_median, + theoretical_median, + epsilon = 0.01, + max_relative = 0.01 + ); + + // Mode (we can't reliably estimate the mode from samples, so we'll check if it's within the range) + assert!(normal.mode() >= min && normal.mode() <= max); + + // Min and Max + let sample_min = data.min(); + let sample_max = data.max(); + assert!(sample_min >= normal.min()); + assert!(sample_max <= normal.max()); + + // CDF and Inverse CDF + let percentiles = [0.1, 0.25, 0.5, 0.75, 0.9]; + for &p in percentiles.iter() { + let sample_percentile = data.percentile((p * 100.0) as usize); + let theoretical_percentile = normal.inverse_cdf(p); + assert_relative_eq!( + sample_percentile, + theoretical_percentile, + epsilon = 0.05, + max_relative = 0.05 + ); + + let sample_cdf = + samples.iter().filter(|&x| x <= &sample_percentile).count() as f64 / n as f64; + let theoretical_cdf = normal.cdf(sample_percentile); + assert_relative_eq!( + sample_cdf, + theoretical_cdf, + epsilon = 0.01, + max_relative = 0.01 + ); + } + + // PDF (we'll check a few points within the range) + let pdf_points = [min, (min + likely) / 2.0, likely, (likely + max) / 2.0, max]; + for &x in pdf_points.iter() { + let sample_pdf = samples.iter().filter(|&s| (s - x).abs() < 0.1).count() as f64 + / (n as f64 * 0.2); + let theoretical_pdf = normal.pdf(x); + assert_relative_eq!( + sample_pdf, + theoretical_pdf, + epsilon = 0.1, + max_relative = 0.1 + ); + } + + println!("Test passed for Normal({}, {}, {})", min, likely, max); + } + } +} diff --git a/src/distributions/pert.rs b/src/distributions/pert.rs new file mode 100644 index 0000000..67c16fd --- /dev/null +++ b/src/distributions/pert.rs @@ -0,0 +1,961 @@ +//! PERT (Program Evaluation and Review Technique) distribution implementation. +//! +//! This module provides an implementation of the PERT distribution, which is commonly used +//! in project management and risk analysis. The PERT distribution is a special case of the +//! Beta distribution, defined by its minimum, most likely (mode), and maximum values. +//! +//! # Examples +//! +//! ``` +//! use distimate::prelude::*; +//! use distimate::Pert; +//! use approx::assert_relative_eq; +//! +//! let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); +//! assert_eq!(pert.mode(), 2.0); +//! assert_relative_eq!(pert.mean().unwrap(), 2.0, epsilon = 1e-6); +//! ``` + +use crate::prelude::*; +use statrs::distribution::Beta; + +/// Represents a PERT (Program Evaluation and Review Technique) distribution. +/// +/// The PERT distribution is defined by its minimum, most likely (mode), and maximum values, +/// along with an optional shape parameter. +/// +/// # Examples +/// +/// ``` +/// use distimate::prelude::*; +/// use distimate::Pert; +/// use approx::assert_relative_eq; +/// +/// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); +/// assert_eq!(pert.min(), 1.0); +/// assert_eq!(pert.mode(), 2.0); +/// assert_eq!(pert.max(), 3.0); +/// assert_relative_eq!(pert.mean().unwrap(), 2.0, epsilon = 1e-6); +/// ``` +#[derive(Debug, Clone)] +pub struct Pert { + min: f64, + mode: f64, + max: f64, + shape: f64, + inner: Beta, +} + +#[allow(dead_code)] +impl Pert { + /// Creates a new PERT distribution with the given minimum, likely, and maximum values. + /// + /// This method uses the default shape parameter of 4.0. + /// + /// # Arguments + /// + /// * `min` - The minimum value of the distribution + /// * `likely` - The most likely value (mode) of the distribution + /// * `max` - The maximum value of the distribution + /// + /// # Returns + /// + /// Returns a `Result` containing the new `Pert` instance if the parameters are valid, + /// or an `Error` if the parameters are invalid. + /// + /// # Errors + /// + /// Returns an error if: + /// - `min` is greater than or equal to `max` + /// - `likely` is less than `min` or greater than `max` + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_eq!(pert.min(), 1.0); + /// assert_eq!(pert.mode(), 2.0); + /// assert_eq!(pert.max(), 3.0); + /// ``` + pub fn new(min: f64, likely: f64, max: f64) -> Result { + Self::new_with_shape(min, likely, max, 4.0) + } + + /// Creates a new PERT distribution with the given minimum, likely, maximum, and shape values. + /// + /// # Arguments + /// + /// * `min` - The minimum value of the distribution + /// * `likely` - The most likely value (mode) of the distribution + /// * `max` - The maximum value of the distribution + /// * `shape` - The shape parameter of the distribution (must be between 2.0 and 6.0, inclusive) + /// + /// # Returns + /// + /// Returns a `Result` containing the new `Pert` instance if the parameters are valid, + /// or an `Error` if the parameters are invalid. + /// + /// # Errors + /// + /// Returns an error if: + /// - `min` is greater than or equal to `max` + /// - `likely` is less than `min` or greater than `max` + /// - `shape` is not between 2.0 and 6.0 (inclusive) + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// + /// let pert = Pert::new_with_shape(1.0, 2.0, 3.0, 5.0).unwrap(); + /// assert_eq!(pert.min(), 1.0); + /// assert_eq!(pert.mode(), 2.0); + /// assert_eq!(pert.max(), 3.0); + /// ``` + pub fn new_with_shape(min: f64, likely: f64, max: f64, shape: f64) -> Result { + // validate inputs + if min >= max { + return Err(Error::InvalidRange { min, max }); + } + if likely < min || likely > max { + return Err(Error::InvalidMode { min, likely, max }); + } + if !(2.0..=6.0).contains(&shape) { + return Err(Error::InvalidShape(shape)); + } + + let alpha = 1.0 + shape * ((likely - min) / (max - min)); + let beta = 1.0 + shape * ((max - likely) / (max - min)); + + Ok(Pert { + min, + mode: likely, + max, + shape, + inner: Beta::new(alpha, beta)?, + }) + } + + /// Returns the alpha parameter of the underlying Beta distribution. + /// + /// The alpha parameter is calculated based on the minimum, mode, maximum, + /// and shape parameters of the PERT distribution. + /// + /// # Returns + /// + /// The alpha parameter as an f64 value. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// use approx::assert_relative_eq; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(pert.alpha(), 3.0, epsilon = 1e-6); + /// ``` + pub fn alpha(&self) -> f64 { + self.inner.shape_a() + } + + /// Returns the beta parameter of the underlying Beta distribution. + /// + /// The beta parameter is calculated based on the minimum, mode, maximum, + /// and shape parameters of the PERT distribution. + /// + /// # Returns + /// + /// The beta parameter as an f64 value. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// use approx::assert_relative_eq; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(pert.beta(), 3.0, epsilon = 1e-6); + /// ``` + pub fn beta(&self) -> f64 { + self.inner.shape_b() + } + + /// Returns the shape parameter of the PERT distribution. + /// + /// The shape parameter influences the peakedness of the distribution. + /// A higher value results in a more peaked distribution around the mode. + /// + /// # Returns + /// + /// The shape parameter as an f64 value. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// + /// let pert = Pert::new_with_shape(1.0, 2.0, 3.0, 5.0).unwrap(); + /// assert_eq!(pert.shape(), 5.0); + /// ``` + pub fn shape(&self) -> f64 { + self.shape + } + + /// Calculates the excess kurtosis of the PERT distribution. + /// + /// Excess kurtosis is a measure of the "tailedness" of the probability distribution + /// compared to a normal distribution. A positive excess kurtosis indicates heavier tails + /// and a higher, sharper peak, while a negative excess kurtosis indicates lighter tails + /// and a lower, more rounded peak. + /// + /// # Returns + /// + /// The excess kurtosis as an f64 value. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// use approx::assert_relative_eq; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(pert.kurtosis(), -0.6666667, epsilon = 1e-6); + /// ``` + pub fn kurtosis(&self) -> f64 { + let (a, b) = (self.alpha(), self.beta()); + (6.0 * ((a - b).powi(2) * (a + b + 1.0) - a * b * (a + b + 2.0))) + / (a * b * (a + b + 2.0) * (a + b + 3.0)) + } +} + +/// Implementation of the `EstimationDistribution` trait for the PERT distribution. +/// +/// This trait marks the PERT distribution as suitable for use in estimation contexts. +/// It doesn't add any new methods, but signals that this distribution is appropriate +/// for modeling uncertain estimates or durations. +/// +/// # Examples +/// +/// ``` +/// use distimate::prelude::*; +/// use distimate::Pert; +/// +/// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); +/// // The fact that we can create a Pert instance demonstrates +/// // that it implements EstimationDistribution +/// ``` +impl EstimationDistribution for Pert {} + +/// Implementation of the `Distribution` trait for the PERT distribution. +/// +/// This implementation provides methods to calculate various statistical properties +/// of the PERT distribution. +impl Distribution for Pert { + /// Calculates the mean of the PERT distribution. + /// + /// # Returns + /// + /// The mean of the distribution as an `Option`. Always returns `Some(value)` + /// as the PERT distribution always has a defined mean. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// use approx::assert_relative_eq; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(pert.mean().unwrap(), 2.0, epsilon = 1e-6); + /// ``` + fn mean(&self) -> Option { + self.inner + .mean() + .map(|m| self.min + (self.max - self.min) * m) + } + + /// Calculates the variance of the PERT distribution. + /// + /// # Returns + /// + /// The variance of the distribution as an `Option`. Always returns `Some(value)` + /// as the PERT distribution always has a defined variance. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// use approx::assert_relative_eq; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(pert.variance().unwrap(), 1.0 / 7.0, epsilon = 1e-6); + /// ``` + fn variance(&self) -> Option { + self.inner + .variance() + .map(|v| (self.max - self.min).powi(2) * v) + } + + /// Calculates the skewness of the PERT distribution. + /// + /// # Returns + /// + /// The skewness of the distribution as an `Option`. Always returns `Some(value)` + /// as the PERT distribution always has a defined skewness. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// use approx::assert_relative_eq; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(pert.skewness().unwrap(), 0.0, epsilon = 1e-6); + /// ``` + fn skewness(&self) -> Option { + self.inner.skewness() + } + + /// Calculates the entropy of the Pert distribution. + /// + /// The entropy is a measure of the average amount of information contained in the distribution. + /// + /// # Returns + /// + /// The entropy of the distribution as an `Option`. Always returns `Some(value)` + /// as the Pert distribution always has a defined entropy. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// use approx::assert_relative_eq; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(pert.entropy().unwrap(), -0.267864, epsilon = 1e-6); + /// ``` + fn entropy(&self) -> Option { + self.inner.entropy() + } +} + +/// Implementation of the `Median` trait for the PERT distribution. +/// +/// The median of a PERT distribution is the value that separates the lower and upper halves +/// of the probability distribution. For a PERT distribution, this is calculated using +/// the inverse cumulative distribution function (inverse CDF) at probability 0.5. +impl Median for Pert { + /// Calculates the median of the PERT distribution. + /// + /// The median is the value separating the higher half from the lower half of the distribution. + /// + /// # Returns + /// + /// The median of the distribution as an f64 value. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// use approx::assert_relative_eq; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(pert.median(), 2.0, epsilon = 1e-4); + /// ``` + fn median(&self) -> f64 { + self.min + (self.max - self.min) * self.inner.inverse_cdf(0.5) + } +} + +/// Implementation of the `Mode` trait for the PERT distribution. +impl Mode for Pert { + /// Returns the mode of the PERT distribution. + /// + /// The mode is the most likely value of the distribution, which is one of the + /// defining parameters of the PERT distribution. + /// + /// # Returns + /// + /// The mode of the distribution as an f64 value. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_eq!(pert.mode(), 2.0); + /// ``` + fn mode(&self) -> f64 { + self.mode + } +} + +/// Implementation of the `Continuous` trait for the PERT distribution. +/// +/// This implementation provides methods to calculate the probability density function (PDF) +/// and its natural logarithm for the PERT distribution. These functions are crucial for +/// understanding the likelihood of different outcomes within the distribution's range. +/// +/// The PERT distribution is continuous over its defined range (minimum to maximum), +/// with the PDF representing the relative likelihood of each possible value. +impl Continuous for Pert { + /// Calculates the probability density function (PDF) of the PERT distribution at a given point. + /// + /// The PDF represents the relative likelihood of the distribution taking on a specific value. + /// For the PERT distribution, the PDF is non-zero only within the range [min, max] and + /// reaches its peak at the mode. + /// + /// # Arguments + /// + /// * `x` - The point at which to calculate the PDF + /// + /// # Returns + /// + /// The value of the PDF at the given point as an f64. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// use approx::assert_relative_eq; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(pert.pdf(2.0), 0.9375, epsilon = 1e-6); + /// assert_eq!(pert.pdf(0.0), 0.0); // Outside the distribution's range + /// ``` + fn pdf(&self, x: f64) -> f64 { + if x < self.min || x > self.max { + 0.0 + } else { + let scaled_x = (x - self.min) / (self.max - self.min); + self.inner.pdf(scaled_x) / (self.max - self.min) + } + } + + /// Calculates the natural logarithm of the probability density function (PDF) of the PERT distribution at a given point. + /// + /// This method is particularly useful for numerical stability in calculations involving very small probability densities. + /// It's often used in statistical inference and machine learning algorithms. + /// + /// # Arguments + /// + /// * `x` - The point at which to calculate the log PDF + /// + /// # Returns + /// + /// The natural logarithm of the PDF at the given point as an f64. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// use approx::assert_relative_eq; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(pert.ln_pdf(2.0), -0.064538, epsilon = 1e-6); + /// assert_eq!(pert.ln_pdf(0.0), f64::NEG_INFINITY); // Outside the distribution's range + /// ``` + fn ln_pdf(&self, x: f64) -> f64 { + if x < self.min || x > self.max { + f64::NEG_INFINITY + } else { + let scaled_x = (x - self.min) / (self.max - self.min); + self.inner.ln_pdf(scaled_x) - (self.max - self.min).ln() + } + } +} + +/// Implementation of the `ContinuousCDF` trait for the PERT distribution. +/// +/// This implementation provides methods to calculate the cumulative distribution function (CDF) +/// and its inverse for the PERT distribution. These functions are essential for understanding +/// the probability of the distribution taking on a value less than or equal to a given point, +/// and for generating random variates from the distribution. +impl ContinuousCDF for Pert { + /// Calculates the cumulative distribution function (CDF) of the PERT distribution at a given point. + /// + /// The CDF represents the probability that a random variable from this distribution + /// will be less than or equal to the given value. + /// + /// # Arguments + /// + /// * `x` - The point at which to calculate the CDF + /// + /// # Returns + /// + /// The value of the CDF at the given point as an f64, in the range [0, 1]. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// use approx::assert_relative_eq; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(pert.cdf(2.0), 0.5, epsilon = 1e-6); + /// assert_eq!(pert.cdf(1.0), 0.0); + /// assert_eq!(pert.cdf(3.0), 1.0); + /// ``` + fn cdf(&self, x: f64) -> f64 { + if x <= self.min { + 0.0 + } else if x >= self.max { + 1.0 + } else { + let scaled_x = (x - self.min) / (self.max - self.min); + self.inner.cdf(scaled_x) + } + } + + /// Calculates the inverse of the cumulative distribution function (inverse CDF) + /// of the PERT distribution for a given probability. + /// + /// This function is also known as the quantile function. It returns the value x + /// for which P(X ≤ x) = p, where p is the given probability. + /// + /// # Arguments + /// + /// * `p` - The probability for which to calculate the inverse CDF, must be in the range [0, 1] + /// + /// # Returns + /// + /// The value x for which the CDF of the distribution equals the given probability. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// use approx::assert_relative_eq; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(pert.inverse_cdf(0.5), 2.0, epsilon = 1e-4); + /// assert_relative_eq!(pert.inverse_cdf(0.0), 1.0, epsilon = 1e-4); + /// assert_relative_eq!(pert.inverse_cdf(1.0), 3.0, epsilon = 1e-4); + /// ``` + fn inverse_cdf(&self, p: f64) -> f64 { + if p <= 0.0 { + self.min + } else if p >= 1.0 { + self.max + } else { + self.min + (self.max - self.min) * self.inner.inverse_cdf(p) + } + } +} + +/// Implementation of the `Min` trait for the PERT distribution. +/// +/// This implementation provides a method to retrieve the minimum value +/// of the PERT distribution, which is one of its defining parameters. +impl Min for Pert { + /// Returns the minimum value of the PERT distribution. + /// + /// The minimum value represents the lowest possible outcome in the distribution. + /// + /// # Returns + /// + /// The minimum value of the distribution as an f64. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_eq!(pert.min(), 1.0); + /// ``` + fn min(&self) -> f64 { + self.min + } +} + +/// Implementation of the `Max` trait for the PERT distribution. +/// +/// This implementation provides a method to retrieve the maximum value +/// of the PERT distribution, which is one of its defining parameters. +impl Max for Pert { + /// Returns the maximum value of the PERT distribution. + /// + /// The maximum value represents the highest possible outcome in the distribution. + /// + /// # Returns + /// + /// The maximum value of the distribution as an f64. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_eq!(pert.max(), 3.0); + /// ``` + fn max(&self) -> f64 { + self.max + } +} + +/// Implementation of the `rand::distributions::Distribution` trait for the PERT distribution. +/// +/// This implementation allows for random sampling from the PERT distribution using +/// any random number generator that implements the `rand::Rng` trait. +impl RandDistribution for Pert { + /// Generates a random sample from the PERT distribution. + /// + /// This method uses the inverse transform sampling technique: it generates a uniform + /// random number between 0 and 1, then applies the inverse CDF to this number to + /// produce a sample from the PERT distribution. + /// + /// # Arguments + /// + /// * `rng` - A random number generator that implements the `rand::Rng` trait + /// + /// # Returns + /// + /// A random sample from the PERT distribution as an f64. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// use rand::prelude::*; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// let mut rng = StdRng::seed_from_u64(42); // For reproducibility + /// + /// let sample = pert.sample(&mut rng); + /// assert!(sample >= 1.0 && sample <= 3.0); + /// ``` + fn sample(&self, rng: &mut R) -> f64 { + let u: f64 = rng.gen(); + self.inverse_cdf(u) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + use rand::distributions::Distribution as RandDistribution; + use rand::rngs::StdRng; + use rand::SeedableRng; + use statrs::statistics::{Data, Distribution, OrderStatistics}; + + #[test] + fn test_new() { + let pert = Pert::new_with_shape(1.0, 2.0, 3.0, 4.0).unwrap(); + assert_eq!(pert.min(), 1.0); + assert_eq!(pert.mode(), 2.0); + assert_eq!(pert.max(), 3.0); + } + + #[test] + fn test_min() { + let pert = Pert::new_with_shape(1.0, 2.0, 3.0, 4.0).unwrap(); + assert_eq!(pert.min(), 1.0); + } + + #[test] + fn test_max() { + let pert = Pert::new_with_shape(1.0, 2.0, 3.0, 4.0).unwrap(); + assert_eq!(pert.max(), 3.0); + } + + #[test] + fn test_alpha_beta() { + let pert = Pert::new_with_shape(1.0, 2.0, 3.0, 4.0).unwrap(); + assert_relative_eq!(pert.alpha(), 3.0, epsilon = 1e-6); + assert_relative_eq!(pert.beta(), 3.0, epsilon = 1e-6); + } + + #[test] + fn test_mean() { + let pert = Pert::new_with_shape(1.0, 2.0, 3.0, 4.0).unwrap(); + assert_relative_eq!(pert.mean().unwrap(), 2.0, epsilon = 1e-6); + } + + #[test] + fn test_median() { + let pert = Pert::new_with_shape(1.0, 2.0, 3.0, 4.0).unwrap(); + assert_relative_eq!(pert.median(), 2.0, epsilon = 1e-3); + } + + #[test] + fn test_mode() { + let pert = Pert::new_with_shape(1.0, 2.0, 3.0, 4.0).unwrap(); + assert_eq!(pert.mode(), 2.0); + } + + #[test] + fn test_variance() { + let pert = Pert::new_with_shape(1.0, 2.0, 3.0, 4.0).unwrap(); + assert_relative_eq!(pert.variance().unwrap(), 1.0 / 7.0, epsilon = 1e-6); + } + + #[test] + fn test_skewness() { + let pert = Pert::new_with_shape(1.0, 2.0, 3.0, 4.0).unwrap(); + assert_relative_eq!(pert.skewness().unwrap(), 0.0, epsilon = 1e-6); + } + + #[test] + fn test_kurtosis() { + let pert = Pert::new_with_shape(1.0, 2.0, 3.0, 4.0).unwrap(); + assert_relative_eq!(pert.kurtosis(), -0.66666667, epsilon = 1e-6); + } + + #[test] + fn test_asymmetric_distribution() { + let pert = Pert::new_with_shape(0.0, 1.0, 5.0, 4.0).unwrap(); + assert_relative_eq!(pert.alpha(), 1.8, epsilon = 1e-6); + assert_relative_eq!(pert.beta(), 4.2, epsilon = 1e-6); + assert_relative_eq!(pert.mean().unwrap(), 1.5, epsilon = 1e-6); + assert_relative_eq!(pert.median(), 1.375, epsilon = 1e-2); + assert_relative_eq!(pert.variance().unwrap(), 0.75, epsilon = 1e-6); + // assert_relative_eq!(pert.skewness().unwrap(), 0.5773502, epsilon = 1e-6); + assert_relative_eq!(pert.kurtosis(), -0.2222222, epsilon = 1e-6); + } + + #[test] + #[should_panic( + expected = "called `Result::unwrap()` on an `Err` value: InvalidRange { min: 3.0, max: 1.0 }" + )] + fn test_invalid_parameters() { + Pert::new_with_shape(3.0, 2.0, 1.0, 4.0).unwrap(); // min > max + } + + fn sample_statistics(dist: &Pert, n: usize) -> (f64, f64, f64, f64) { + let mut rng = StdRng::seed_from_u64(42); // Fixed seed for reproducibility + let samples: Vec = (0..n).map(|_| dist.sample(&mut rng)).collect(); + + let mean = samples.iter().sum::() / n as f64; + let variance = samples.iter().map(|&x| (x - mean).powi(2)).sum::() / (n - 1) as f64; + let std_dev = variance.sqrt(); + let mut sorted_samples = samples.clone(); + sorted_samples.sort_by(|a, b| a.partial_cmp(b).unwrap()); + let median = if n % 2 == 0 { + (sorted_samples[n / 2 - 1] + sorted_samples[n / 2]) / 2.0 + } else { + sorted_samples[n / 2] + }; + + (mean, variance, std_dev, median) + } + + #[test] + fn test_pert_distribution_statistics() { + let min = 1.0; + let mode = 2.0; + let max = 5.0; + let shape = 4.0; + let dist = Pert::new_with_shape(min, mode, max, shape).unwrap(); + let n = 100_000; // Number of samples + let epsilon = 0.01; // Tolerance for relative equality + + let (sample_mean, sample_variance, sample_std_dev, sample_median) = + sample_statistics(&dist, n); + + // Test mean + assert_relative_eq!( + sample_mean, + dist.mean().unwrap(), + epsilon = epsilon, + max_relative = epsilon + ); + + // Test variance + assert_relative_eq!( + sample_variance, + dist.variance().unwrap(), + epsilon = epsilon, + max_relative = epsilon + ); + + // Test standard deviation + assert_relative_eq!( + sample_std_dev, + dist.std_dev().unwrap(), + epsilon = epsilon, + max_relative = epsilon + ); + + // Test median + assert_relative_eq!( + sample_median, + dist.median(), + epsilon = epsilon, + max_relative = epsilon + ); + } + + #[test] + fn test_pert_distribution_bounds() { + let min = 1.0; + let mode = 2.0; + let max = 3.0; + let shape = 4.0; + let dist = Pert::new_with_shape(min, mode, max, shape).unwrap(); + let n = 100_000; + let mut rng = StdRng::seed_from_u64(42); + + let samples: Vec = (0..n).map(|_| dist.sample(&mut rng)).collect(); + + assert!(samples.iter().all(|&x| x >= min && x <= max)); + assert!(samples.iter().any(|&x| x < mode)); + assert!(samples.iter().any(|&x| x > mode)); + } + + fn sample_skewness(data: &[f64]) -> f64 { + let n = data.len() as f64; + let mean = data.iter().sum::() / n; + let m2 = data.iter().map(|x| (x - mean).powi(2)).sum::() / n; + let m3 = data.iter().map(|x| (x - mean).powi(3)).sum::() / n; + m3 / m2.powf(1.5) + } + + fn sample_kurtosis(data: &[f64]) -> f64 { + let n = data.len() as f64; + let mean = data.iter().sum::() / n; + let m2 = data.iter().map(|x| (x - mean).powi(2)).sum::() / n; + let m4 = data.iter().map(|x| (x - mean).powi(4)).sum::() / n; + m4 / m2.powi(2) - 3.0 + } + + #[test] + fn test_distribution_statistics() { + let test_cases = vec![ + (0.0, 5.0, 10.0, 4.0), // Standard PERT + (0.0, 2.0, 10.0, 2.0), // Modified PERT with lower shape + (0.0, 9.0, 10.0, 6.0), // Modified PERT with higher shape + (-10.0, 0.0, 10.0, 4.0), // Standard PERT with negative values + (100.0, 400.0, 1000.0, 4.0), // Standard PERT with larger values + ]; + + for (min, mode, max, shape) in test_cases { + let pert = Pert::new_with_shape(min, mode, max, shape).unwrap(); + let n = 25_000; // Number of samples + let mut rng = StdRng::seed_from_u64(42); // Fixed seed for reproducibility + + let samples: Vec = (0..n).map(|_| pert.sample(&mut rng)).collect(); + let mut data = Data::new(samples.clone()); + + // Mean + let sample_mean = data.mean().unwrap(); + let theoretical_mean = pert.mean().unwrap(); + assert_relative_eq!( + sample_mean, + theoretical_mean, + epsilon = 0.1, + max_relative = 0.1 + ); + + // Variance + let sample_variance = data.variance().unwrap(); + let theoretical_variance = pert.variance().unwrap(); + assert_relative_eq!( + sample_variance, + theoretical_variance, + epsilon = 0.05, + max_relative = 0.05 + ); + + // Standard deviation + let sample_std_dev = data.std_dev().unwrap(); + let theoretical_std_dev = pert.std_dev().unwrap(); + assert_relative_eq!( + sample_std_dev, + theoretical_std_dev, + epsilon = 0.1, + max_relative = 0.1 + ); + + // Skewness + let sample_skewness = sample_skewness(&samples); + let theoretical_skewness = pert.skewness().unwrap(); + assert_relative_eq!( + sample_skewness, + theoretical_skewness, + epsilon = 0.1, + max_relative = 0.1 + ); + + // Median + let sample_median = data.median(); + let theoretical_median = pert.median(); + assert_relative_eq!( + sample_median, + theoretical_median, + epsilon = 0.1, + max_relative = 0.1 + ); + + // Mode (we can't reliably estimate the mode from samples, so we'll check if it's within the range) + assert!(pert.mode() >= min && pert.mode() <= max); + + // Min and Max + let sample_min = data.min(); + let sample_max = data.max(); + assert!(sample_min >= pert.min()); + assert!(sample_max <= pert.max()); + + // Excess Kurtosis + let sample_kurtosis = sample_kurtosis(&samples); + let theoretical_kurtosis = pert.kurtosis(); + assert_relative_eq!( + sample_kurtosis, + theoretical_kurtosis, + epsilon = 0.1, + max_relative = 0.1 + ); + + // CDF and Inverse CDF + let percentiles = [0.1, 0.25, 0.5, 0.75, 0.9]; + for p in percentiles.iter() { + let sample_percentile = data.percentile((100.0 * p) as usize); + let theoretical_percentile = pert.inverse_cdf(*p); + assert_relative_eq!( + sample_percentile, + theoretical_percentile, + epsilon = 0.05, + max_relative = 0.05 + ); + + let sample_cdf = + samples.iter().filter(|&x| x <= &sample_percentile).count() as f64 / n as f64; + let theoretical_cdf = pert.cdf(sample_percentile); + assert_relative_eq!( + sample_cdf, + theoretical_cdf, + epsilon = 0.01, + max_relative = 0.01 + ); + } + + println!( + "Test passed for PERT({}, {}, {}, {})", + min, mode, max, shape + ); + } + } +} diff --git a/src/distributions/triangular.rs b/src/distributions/triangular.rs new file mode 100644 index 0000000..742318d --- /dev/null +++ b/src/distributions/triangular.rs @@ -0,0 +1,712 @@ +//! Triangular distribution implementation. +//! +//! This module provides an implementation of the Triangular distribution, which is often used +//! in project management, risk analysis, and other fields where limited sample data is available. +//! The Triangular distribution is defined by its minimum, most likely (mode), and maximum values. +//! +//! The Triangular distribution is particularly useful when the relationship between variables is known, +//! but data is scarce. It provides a simple and intuitive way to model situations where you can estimate +//! the minimum, maximum, and most likely values. +//! +//! # Examples +//! +//! ``` +//! use distimate::prelude::*; +//! use distimate::Triangular; +//! use approx::assert_relative_eq; +//! +//! let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); +//! assert_eq!(triangular.mode(), 2.0); +//! assert_relative_eq!(triangular.mean().unwrap(), 2.0, epsilon = 1e-6); +//! ``` + +use crate::prelude::*; +use statrs::distribution::Triangular as StatrsTriangular; + +/// Represents a Triangular distribution. +/// +/// The Triangular distribution is defined by its minimum, most likely (mode), and maximum values. +/// +/// # Examples +/// +/// ``` +/// use distimate::prelude::*; +/// use distimate::Triangular; +/// use approx::assert_relative_eq; +/// +/// let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); +/// assert_eq!(triangular.min(), 1.0); +/// assert_eq!(triangular.mode(), 2.0); +/// assert_eq!(triangular.max(), 3.0); +/// assert_relative_eq!(triangular.mean().unwrap(), 2.0, epsilon = 1e-4); +/// ``` +#[derive(Debug, Clone)] +pub struct Triangular { + inner: StatrsTriangular, +} + +impl Triangular { + /// Creates a new Triangular distribution with the given minimum, likely, and maximum values. + /// + /// # Arguments + /// + /// * `min` - The minimum value of the distribution + /// * `likely` - The most likely value (mode) of the distribution + /// * `max` - The maximum value of the distribution + /// + /// # Returns + /// + /// Returns a `Result` containing the new `Triangular` instance if the parameters are valid, + /// or an `Error` if the parameters are invalid. + /// + /// # Errors + /// + /// Returns an error if: + /// - `min` is greater than or equal to `max` + /// - `likely` is less than `min` or greater than `max` + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Triangular; + /// + /// let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + /// assert_eq!(triangular.min(), 1.0); + /// assert_eq!(triangular.mode(), 2.0); + /// assert_eq!(triangular.max(), 3.0); + /// ``` + pub fn new(min: f64, likely: f64, max: f64) -> Result { + if min >= max { + return Err(Error::InvalidRange { min, max }); + } + if likely < min || likely > max { + return Err(Error::InvalidMode { min, likely, max }); + } + + Ok(Triangular { + inner: StatrsTriangular::new(min, max, likely)?, + }) + } +} + +/// Implementation of the `EstimationDistribution` trait for the Triangular distribution. +/// +/// This trait marks the Triangular distribution as suitable for use in estimation contexts. +/// It doesn't add any new methods, but signals that this distribution is appropriate +/// for modeling uncertain estimates or durations. +/// +/// # Examples +/// +/// ``` +/// use distimate::prelude::*; +/// use distimate::Triangular; +/// +/// let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); +/// // The fact that we can create a Triangular instance demonstrates +/// // that it implements EstimationDistribution +/// ``` +impl EstimationDistribution for Triangular {} + +/// Implementation of the `Distribution` trait for the Triangular distribution. +/// +/// This implementation provides methods to calculate various statistical properties +/// of the Triangular distribution. +impl Distribution for Triangular { + /// Calculates the mean of the Triangular distribution. + /// + /// # Returns + /// + /// The mean of the distribution as an `Option`. Always returns `Some(value)` + /// as the Triangular distribution always has a defined mean. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Triangular; + /// use approx::assert_relative_eq; + /// + /// let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(triangular.mean().unwrap(), 2.0, epsilon = 1e-6); + /// ``` + fn mean(&self) -> Option { + self.inner.mean() + } + + /// Calculates the variance of the Triangular distribution. + /// + /// # Returns + /// + /// The variance of the distribution as an `Option`. Always returns `Some(value)` + /// as the Triangular distribution always has a defined variance. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Triangular; + /// use approx::assert_relative_eq; + /// + /// let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(triangular.variance().unwrap(), 0.1666667, epsilon = 1e-6); + /// ``` + fn variance(&self) -> Option { + self.inner.variance() + } + + /// Calculates the skewness of the Triangular distribution. + /// + /// # Returns + /// + /// The skewness of the distribution as an `Option`. Always returns `Some(value)` + /// as the Triangular distribution always has a defined skewness. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Triangular; + /// use approx::assert_relative_eq; + /// + /// let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(triangular.skewness().unwrap(), 0.0, epsilon = 1e-6); + /// ``` + fn skewness(&self) -> Option { + self.inner.skewness() + } + + /// Calculates the entropy of the Triangular distribution. + /// + /// The entropy is a measure of the average amount of information contained in the distribution. + /// + /// # Returns + /// + /// The entropy of the distribution as an `Option`. Always returns `Some(value)` + /// as the Triangular distribution always has a defined entropy. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Triangular; + /// use approx::assert_relative_eq; + /// + /// let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(triangular.entropy().unwrap(), 0.5, epsilon = 1e-6); + /// ``` + fn entropy(&self) -> Option { + self.inner.entropy() + } +} + +/// Implementation of the `Median` trait for the Triangular distribution. +impl Median for Triangular { + /// Calculates the median of the Triangular distribution. + /// + /// The median is the value separating the higher half from the lower half of the distribution. + /// + /// # Returns + /// + /// The median of the distribution as an f64 value. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Triangular; + /// use approx::assert_relative_eq; + /// + /// let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(triangular.median(), 2.0, epsilon = 1e-6); + /// ``` + fn median(&self) -> f64 { + self.inner.median() + } +} + +/// Implementation of the `Mode` trait for the Triangular distribution. +impl Mode for Triangular { + /// Returns the mode of the Triangular distribution. + /// + /// The mode is the most likely value of the distribution, which is one of the + /// defining parameters of the Triangular distribution. + /// + /// # Returns + /// + /// The mode of the distribution as an f64 value. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Triangular; + /// + /// let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + /// assert_eq!(triangular.mode(), 2.0); + /// ``` + fn mode(&self) -> f64 { + self.inner.mode().unwrap() + } +} + +/// Implementation of the `Continuous` trait for the Triangular distribution. +/// +/// This implementation provides methods to calculate the probability density function (PDF) +/// and its natural logarithm for the Triangular distribution. +impl Continuous for Triangular { + /// Calculates the probability density function (PDF) of the Triangular distribution at a given point. + /// + /// # Arguments + /// + /// * `x` - The point at which to calculate the PDF + /// + /// # Returns + /// + /// The value of the PDF at the given point as an f64. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Triangular; + /// use approx::assert_relative_eq; + /// + /// let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(triangular.pdf(2.0), 1.0, epsilon = 1e-6); + /// assert_eq!(triangular.pdf(0.0), 0.0); // Outside the distribution's range + /// ``` + fn pdf(&self, x: f64) -> f64 { + self.inner.pdf(x) + } + + /// Calculates the natural logarithm of the probability density function (PDF) of the Triangular distribution at a given point. + /// + /// # Arguments + /// + /// * `x` - The point at which to calculate the log PDF + /// + /// # Returns + /// + /// The natural logarithm of the PDF at the given point as an f64. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Triangular; + /// use approx::assert_relative_eq; + /// + /// let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(triangular.ln_pdf(2.0), 0.0, epsilon = 1e-6); + /// assert_eq!(triangular.ln_pdf(0.0), f64::NEG_INFINITY); // Outside the distribution's range + /// ``` + fn ln_pdf(&self, x: f64) -> f64 { + self.inner.ln_pdf(x) + } +} + +/// Implementation of the `ContinuousCDF` trait for the Triangular distribution. +impl ContinuousCDF for Triangular { + /// Calculates the cumulative distribution function (CDF) of the Triangular distribution at a given point. + /// + /// # Arguments + /// + /// * `x` - The point at which to calculate the CDF + /// + /// # Returns + /// + /// The value of the CDF at the given point as an f64, in the range [0, 1]. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Triangular; + /// use approx::assert_relative_eq; + /// + /// let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(triangular.cdf(2.0), 0.5, epsilon = 1e-6); + /// assert_eq!(triangular.cdf(1.0), 0.0); + /// assert_eq!(triangular.cdf(3.0), 1.0); + /// ``` + fn cdf(&self, x: f64) -> f64 { + self.inner.cdf(x) + } + + /// Calculates the inverse of the cumulative distribution function (inverse CDF) + /// of the Triangular distribution for a given probability. + /// + /// # Arguments + /// + /// * `p` - The probability for which to calculate the inverse CDF, must be in the range [0, 1] + /// + /// # Returns + /// + /// The value x for which the CDF of the distribution equals the given probability. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Triangular; + /// use approx::assert_relative_eq; + /// + /// let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(triangular.inverse_cdf(0.5), 2.0, epsilon = 1e-4); + /// assert_relative_eq!(triangular.inverse_cdf(0.0), 1.0, epsilon = 1e-4); + /// assert_relative_eq!(triangular.inverse_cdf(1.0), 3.0, epsilon = 1e-4); + /// ``` + fn inverse_cdf(&self, p: f64) -> f64 { + self.inner.inverse_cdf(p) + } +} + +/// Implementation of the `Min` trait for the Triangular distribution. +impl Min for Triangular { + /// Returns the minimum value of the Triangular distribution. + /// + /// # Returns + /// + /// The minimum value of the distribution as an f64. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Triangular; + /// + /// let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + /// assert_eq!(triangular.min(), 1.0); + /// ``` + fn min(&self) -> f64 { + self.inner.min() + } +} + +/// Implementation of the `Max` trait for the Triangular distribution. +impl Max for Triangular { + /// Returns the maximum value of the Triangular distribution. + /// + /// # Returns + /// + /// The maximum value of the distribution as an f64. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Triangular; + /// + /// let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + /// assert_eq!(triangular.max(), 3.0); + /// ``` + fn max(&self) -> f64 { + self.inner.max() + } +} + +/// Implementation of the `rand::distributions::Distribution` trait for the Triangular distribution. +/// +/// This implementation allows for random sampling from the Triangular distribution using +/// any random number generator that implements the `rand::Rng` trait. +impl RandDistribution for Triangular { + /// Generates a random sample from the Triangular distribution. + /// + /// This method uses the inverse transform sampling technique: it generates a uniform + /// random number between 0 and 1, then applies the inverse CDF to this number to + /// produce a sample from the Triangular distribution. + /// + /// # Arguments + /// + /// * `rng` - A random number generator that implements the `rand::Rng` trait + /// + /// # Returns + /// + /// A random sample from the Triangular distribution as an f64. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Triangular; + /// use rand::prelude::*; + /// + /// let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + /// let mut rng = StdRng::seed_from_u64(42); // For reproducibility + /// + /// let sample = triangular.sample(&mut rng); + /// assert!(sample >= 1.0 && sample <= 3.0); + /// ``` + fn sample(&self, rng: &mut R) -> f64 { + self.inner.sample(rng) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + use rand::distributions::Distribution as RandDistribution; + use rand::rngs::StdRng; + use rand::SeedableRng; + use statrs::statistics::{Data, Distribution, OrderStatistics}; + + #[test] + fn test_new() { + let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + assert_eq!(triangular.min(), 1.0); + assert_eq!(triangular.mode(), 2.0); + assert_eq!(triangular.max(), 3.0); + } + + #[test] + fn test_min() { + let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + assert_eq!(triangular.min(), 1.0); + } + + #[test] + fn test_max() { + let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + assert_eq!(triangular.max(), 3.0); + } + + #[test] + fn test_mean() { + let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + assert_relative_eq!(triangular.mean().unwrap(), 2.0, epsilon = 1e-6); + } + + #[test] + fn test_median() { + let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + assert_relative_eq!(triangular.median(), 2.0, epsilon = 1e-6); + } + + #[test] + fn test_mode() { + let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + assert_eq!(triangular.mode(), 2.0); + } + + #[test] + fn test_variance() { + let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + assert_relative_eq!(triangular.variance().unwrap(), 0.1666667, epsilon = 1e-6); + } + + #[test] + fn test_skewness() { + let triangular = Triangular::new(1.0, 2.0, 3.0).unwrap(); + assert_relative_eq!(triangular.skewness().unwrap(), 0.0, epsilon = 1e-6); + } + + #[test] + fn test_asymmetric_distribution() { + let triangular = Triangular::new(0.0, 1.0, 5.0).unwrap(); + assert_relative_eq!(triangular.mean().unwrap(), 2.0, epsilon = 1e-6); + assert_relative_eq!(triangular.median(), 1.83772234, epsilon = 1e-6); + assert_relative_eq!(triangular.variance().unwrap(), 1.1666667, epsilon = 1e-6); + assert_relative_eq!(triangular.skewness().unwrap(), 0.47613605, epsilon = 1e-6); + } + + #[test] + #[should_panic( + expected = "called `Result::unwrap()` on an `Err` value: InvalidRange { min: 3.0, max: 1.0 }" + )] + fn test_invalid_parameters() { + Triangular::new(3.0, 2.0, 1.0).unwrap(); // min > max + } + + fn sample_statistics(dist: &Triangular, n: usize) -> (f64, f64, f64, f64) { + let mut rng = StdRng::seed_from_u64(42); // Fixed seed for reproducibility + let samples: Vec = (0..n).map(|_| dist.sample(&mut rng)).collect(); + + let mean = samples.iter().sum::() / n as f64; + let variance = samples.iter().map(|&x| (x - mean).powi(2)).sum::() / (n - 1) as f64; + let std_dev = variance.sqrt(); + let mut sorted_samples = samples.clone(); + sorted_samples.sort_by(|a, b| a.partial_cmp(b).unwrap()); + let median = if n % 2 == 0 { + (sorted_samples[n / 2 - 1] + sorted_samples[n / 2]) / 2.0 + } else { + sorted_samples[n / 2] + }; + + (mean, variance, std_dev, median) + } + + #[test] + fn test_triangular_distribution_statistics() { + let min = 1.0; + let mode = 2.0; + let max = 5.0; + let dist = Triangular::new(min, mode, max).unwrap(); + let n = 100_000; // Number of samples + let epsilon = 0.01; // Tolerance for relative equality + + let (sample_mean, sample_variance, sample_std_dev, sample_median) = + sample_statistics(&dist, n); + + // Test mean + assert_relative_eq!( + sample_mean, + dist.mean().unwrap(), + epsilon = epsilon, + max_relative = epsilon + ); + + // Test variance + assert_relative_eq!( + sample_variance, + dist.variance().unwrap(), + epsilon = epsilon, + max_relative = epsilon + ); + + // Test standard deviation + assert_relative_eq!( + sample_std_dev, + dist.std_dev().unwrap(), + epsilon = epsilon, + max_relative = epsilon + ); + + // Test median + assert_relative_eq!( + sample_median, + dist.median(), + epsilon = epsilon, + max_relative = epsilon + ); + } + + #[test] + fn test_triangular_distribution_bounds() { + let min = 1.0; + let mode = 2.0; + let max = 3.0; + let dist = Triangular::new(min, mode, max).unwrap(); + let n = 100_000; + let mut rng = StdRng::seed_from_u64(42); + + let samples: Vec = (0..n).map(|_| dist.sample(&mut rng)).collect(); + + assert!(samples.iter().all(|&x| x >= min && x <= max)); + assert!(samples.iter().any(|&x| x < mode)); + assert!(samples.iter().any(|&x| x > mode)); + } + + fn sample_skewness(data: &[f64]) -> f64 { + let n = data.len() as f64; + let mean = data.iter().sum::() / n; + let m2 = data.iter().map(|x| (x - mean).powi(2)).sum::() / n; + let m3 = data.iter().map(|x| (x - mean).powi(3)).sum::() / n; + m3 / m2.powf(1.5) + } + + #[test] + fn test_distribution_statistics() { + let test_cases = vec![ + (0.0, 5.0, 10.0), // Standard Triangular + (0.0, 2.0, 10.0), // Right-skewed Triangular + (0.0, 9.0, 10.0), // Left-skewed Triangular + (-10.0, 0.0, 10.0), // Symmetric Triangular with negative values + (100.0, 400.0, 1000.0), // Triangular with larger values + ]; + + for (min, mode, max) in test_cases { + let triangular = Triangular::new(min, mode, max).unwrap(); + let n = 25_000; // Number of samples + let mut rng = StdRng::seed_from_u64(42); // Fixed seed for reproducibility + + let samples: Vec = (0..n).map(|_| triangular.sample(&mut rng)).collect(); + let mut data = Data::new(samples.clone()); + + // Mean + let sample_mean = data.mean().unwrap(); + let theoretical_mean = triangular.mean().unwrap(); + assert_relative_eq!( + sample_mean, + theoretical_mean, + epsilon = 0.1, + max_relative = 0.1 + ); + + // Variance + let sample_variance = data.variance().unwrap(); + let theoretical_variance = triangular.variance().unwrap(); + assert_relative_eq!( + sample_variance, + theoretical_variance, + epsilon = 0.05, + max_relative = 0.05 + ); + + // Standard deviation + let sample_std_dev = data.std_dev().unwrap(); + let theoretical_std_dev = triangular.std_dev().unwrap(); + assert_relative_eq!( + sample_std_dev, + theoretical_std_dev, + epsilon = 0.1, + max_relative = 0.1 + ); + + // Skewness + let sample_skewness = sample_skewness(&samples); + let theoretical_skewness = triangular.skewness().unwrap(); + assert_relative_eq!( + sample_skewness, + theoretical_skewness, + epsilon = 0.1, + max_relative = 0.1 + ); + + // Median + let sample_median = data.median(); + let theoretical_median = triangular.median(); + assert_relative_eq!( + sample_median, + theoretical_median, + epsilon = 0.1, + max_relative = 0.1 + ); + + // Mode (we can't reliably estimate the mode from samples, so we'll check if it's within the range) + assert!(triangular.mode() >= min && triangular.mode() <= max); + + // Min and Max + let sample_min = data.min(); + let sample_max = data.max(); + assert!(sample_min >= triangular.min()); + assert!(sample_max <= triangular.max()); + + // CDF and Inverse CDF + let percentiles = [0.1, 0.25, 0.5, 0.75, 0.9]; + for p in percentiles.iter() { + let sample_percentile = data.percentile((100.0 * p) as usize); + let theoretical_percentile = triangular.inverse_cdf(*p); + assert_relative_eq!( + sample_percentile, + theoretical_percentile, + epsilon = 0.05, + max_relative = 0.05 + ); + + let sample_cdf = + samples.iter().filter(|&x| x <= &sample_percentile).count() as f64 / n as f64; + let theoretical_cdf = triangular.cdf(sample_percentile); + assert_relative_eq!( + sample_cdf, + theoretical_cdf, + epsilon = 0.01, + max_relative = 0.01 + ); + } + + println!("Test passed for Triangular({}, {}, {})", min, mode, max); + } + } +} diff --git a/src/error.rs b/src/error.rs new file mode 100644 index 0000000..294482b --- /dev/null +++ b/src/error.rs @@ -0,0 +1,40 @@ +use thiserror::Error; + +#[derive(Debug, Error)] +pub enum Error { + #[error("Insufficient data: at least {required} data points required, but got {provided}")] + InsufficientData { required: usize, provided: usize }, + + #[error("Invalid data: {0}")] + InvalidData(String), + + #[error("Fitting error: {0}")] + FittingError(String), + + #[error("Fit quality check error: {0}")] + FitQualityError(String), + + #[error("Invalid parameter: {0}")] + InvalidParameter(String), + + #[error("Math error: {0}")] + MathError(String), + + #[error("Invalid percentile: {0}. Must be between 0 and 1.")] + InvalidPercentile(f64), + + #[error("Invalid confidence level: {0}. Must be between 0 and 1.")] + InvalidConfidenceLevel(f64), + + #[error("Invalid shape: {0}. Must be between 0 and 1.")] + InvalidShape(f64), + + #[error("Invalid parameters: min ({min}) must be less than max ({max}).")] + InvalidRange { min: f64, max: f64 }, + + #[error("Likely value ({likely}) must be between min ({min}) and max ({max}).")] + InvalidMode { min: f64, likely: f64, max: f64 }, + + #[error("statrs::StatsError: {0}.")] + StatsError(#[from] statrs::StatsError), +} diff --git a/src/estimation_distribution.rs b/src/estimation_distribution.rs new file mode 100644 index 0000000..5e4ce86 --- /dev/null +++ b/src/estimation_distribution.rs @@ -0,0 +1,296 @@ +use rand::distributions::Distribution as RandDistribution; +use statrs::{ + distribution::{Continuous, ContinuousCDF}, + statistics::{Distribution, Max, Median, Min, Mode}, +}; + +use crate::{Error, EstimationFitQuality, Result}; + +/// A trait for probability distributions used in estimation scenarios. +/// +/// This trait combines several statistical traits and provides additional methods +/// specific to estimation and risk analysis. It's designed to be implemented by +/// distributions that model uncertainty in estimates, such as project durations, +/// costs, or other quantifiable outcomes. +/// +/// # Examples +/// +/// ``` +/// use distimate::prelude::*; +/// use distimate::Pert; +/// use approx::assert_relative_eq; +/// +/// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); +/// assert_relative_eq!(pert.expected_value(), 2.0, epsilon = 1e-6); +/// assert_relative_eq!(pert.optimistic_estimate(), 1.37847, epsilon = 1e-4); +/// assert_relative_eq!(pert.pessimistic_estimate(), 2.62152, epsilon = 1e-4); +/// ``` +pub trait EstimationDistribution: + Distribution + + Continuous + + ContinuousCDF + + Median + + Mode + + Min + + Max + + RandDistribution // TODO: add TryFrom (x2) here +{ + /// Returns the estimate at a given percentile. + /// + /// # Arguments + /// + /// * `p` - The percentile as a float between 0.0 and 1.0. + /// + /// # Returns + /// + /// The estimated value at the given percentile, or an error if the percentile is invalid. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// use approx::assert_relative_eq; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(pert.percentile_estimate(0.5).unwrap(), 2.0, epsilon = 1e-4); + /// ``` + fn percentile_estimate(&self, p: f64) -> Result { + if !(0.0..=1.0).contains(&p) { + Err(Error::InvalidPercentile(p)) + } else { + Ok(self.inverse_cdf(p)) + } + } + + /// Returns the optimistic (best-case) estimate, typically the 5th percentile. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// use approx::assert_relative_eq; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(pert.optimistic_estimate(), 1.37847, epsilon = 1e-4); + /// ``` + fn optimistic_estimate(&self) -> f64 { + self.percentile_estimate(0.05).unwrap() + } + + /// Returns the pessimistic (worst-case) estimate, typically the 95th percentile. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// use approx::assert_relative_eq; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(pert.pessimistic_estimate(), 2.62152, epsilon = 1e-4); + /// ``` + fn pessimistic_estimate(&self) -> f64 { + self.percentile_estimate(0.95).unwrap() + } + + /// Returns the most likely estimate (mode). + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_eq!(pert.most_likely_estimate(), 2.0); + /// ``` + fn most_likely_estimate(&self) -> f64 { + self.mode() + } + + /// Returns the expected value (mean) of the distribution. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// use approx::assert_relative_eq; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(pert.expected_value(), 2.0, epsilon = 1e-6); + /// ``` + fn expected_value(&self) -> f64 { + self.mean().unwrap() + } + + /// Returns the standard deviation of the distribution. + /// + /// # Examples + /// + /// ``` + /// use distimate::prelude::*; + /// use distimate::Pert; + /// use approx::assert_relative_eq; + /// + /// let pert = Pert::new(1.0, 2.0, 3.0).unwrap(); + /// assert_relative_eq!(pert.uncertainty(), 0.377964, epsilon = 1e-4); + /// ``` + fn uncertainty(&self) -> f64 { + self.std_dev().unwrap() + } + + fn probability_of_completion(&self, estimate: f64) -> f64 { + self.cdf(estimate) + } + + /// Calculates the risk of exceeding a given estimate. + /// + /// This method computes the probability that the actual value will be greater than + /// the provided estimate. It's particularly useful in project management and risk + /// analysis to quantify the likelihood of cost overruns or schedule delays. + /// + /// The risk of overrun is calculated as 1 minus the probability of completion: + /// Risk = 1 - P(X <= estimate), where X is the random variable representing the + /// actual cost or time. + /// + /// For example: + /// - A risk of 0.2 (20%) means there's a 20% chance the project will exceed the estimate. + /// - A risk of 0.5 (50%) means it's equally likely to be above or below the estimate. + /// - A risk of 0.8 (80%) suggests a high likelihood of exceeding the estimate. + /// + /// # Arguments + /// + /// * `estimate` - The point estimate to evaluate the risk against. + /// + /// # Returns + /// + /// A float between 0 and 1 representing the probability of exceeding the estimate. + fn risk_of_overrun(&self, estimate: f64) -> f64 { + 1.0 - self.probability_of_completion(estimate) + } + + /// Returns a confidence interval for the estimate. + /// + /// A confidence interval provides a range of values that is likely to contain + /// the true value with a certain level of confidence. It's useful for expressing + /// the uncertainty in an estimate. + /// + /// This method calculates a two-sided confidence interval, meaning it provides + /// both a lower and upper bound. The interval is centered around the median of + /// the distribution. + /// + /// # Arguments + /// + /// * `confidence_level` - A value between 0 and 1 representing the desired level of confidence. + /// Common values are 0.95 for a 95% confidence interval or 0.99 for a 99% confidence interval. + /// + /// # Returns + /// + /// A tuple (lower, upper) representing the lower and upper bounds of the confidence interval. + /// + /// # Example + /// + /// If you call `confidence_interval(0.95)` and get (100, 150), it means: + /// "We are 95% confident that the true value lies between 100 and 150." + /// + /// # Note + /// + /// The width of the interval indicates the precision of the estimate. A narrower interval + /// suggests a more precise estimate, while a wider interval indicates more uncertainty. + fn confidence_interval(&self, confidence_level: f64) -> (f64, f64) { + let lower = self.percentile_estimate((1.0 - confidence_level) / 2.0); + let upper = self.percentile_estimate(1.0 - (1.0 - confidence_level) / 2.0); + (lower.unwrap(), upper.unwrap()) + } + + /// Evaluates how well the distribution fits a given dataset in the context of estimation. + /// + /// This method provides a measure of how accurately the distribution represents + /// the provided set of estimates or actual values. + /// + /// # Arguments + /// + /// * `data` - A slice of f64 values representing historical estimates or actual values. + /// + /// # Returns + /// + /// Returns an `EstimationFitQuality` struct containing various fit metrics. + /// + /// # Errors + /// + /// Returns an `Error` if there's an issue calculating the fit, + /// such as invalid data points or errors in percentile calculations. + fn evaluate_fit_quality(&self, data: &[f64]) -> Result { + let n = data.len(); + if n == 0 { + return Err(Error::InsufficientData { + required: 1, + provided: 0, + }); + } + + let mean = data.iter().sum::() / n as f64; + let median = self.median(); + let mode = self.mode(); + + let mse = data.iter().map(|&x| (x - mean).powi(2)).sum::() / n as f64; + let rmse = mse.sqrt(); + + let mape = data.iter().map(|&x| ((x - mean) / x).abs()).sum::() / n as f64 * 100.0; + + let within_50 = self.calculate_within_interval(data, 0.25, 0.75)?; + let within_68 = self.calculate_within_interval(data, 0.16, 0.84)?; + let within_95 = self.calculate_within_interval(data, 0.025, 0.975)?; + + Ok(EstimationFitQuality { + mean_absolute_percentage_error: mape, + root_mean_square_error: rmse, + distribution_mean: mean, + distribution_median: median, + distribution_mode: mode, + within_50_percent: within_50, + within_68_percent: within_68, + within_95_percent: within_95, + sample_size: n, + }) + } + + /// Helper method to calculate the proportion of data within a given interval. + /// + /// # Arguments + /// + /// * `data` - A slice of f64 values representing the dataset. + /// * `lower_percentile` - The lower bound of the interval as a percentile (0.0 to 1.0). + /// * `upper_percentile` - The upper bound of the interval as a percentile (0.0 to 1.0). + /// + /// # Returns + /// + /// Returns the proportion of data points within the specified interval. + /// + /// # Errors + /// + /// Returns an `Error` if there's an issue calculating the percentiles. + fn calculate_within_interval( + &self, + data: &[f64], + lower_percentile: f64, + upper_percentile: f64, + ) -> Result { + let lower_bound = self + .percentile_estimate(lower_percentile) + .map_err(|e| Error::FitQualityError(format!("Error calculating lower bound: {}", e)))?; + let upper_bound = self + .percentile_estimate(upper_percentile) + .map_err(|e| Error::FitQualityError(format!("Error calculating upper bound: {}", e)))?; + + let count_within = data + .iter() + .filter(|&&x| x >= lower_bound && x <= upper_bound) + .count(); + + Ok(count_within as f64 / data.len() as f64) + } +} diff --git a/src/fit_quality.rs b/src/fit_quality.rs new file mode 100644 index 0000000..60f6ea5 --- /dev/null +++ b/src/fit_quality.rs @@ -0,0 +1,22 @@ +/// Represents the quality of fit for an estimation distribution. +#[derive(Debug, Clone, PartialEq)] +pub struct EstimationFitQuality { + /// Mean Absolute Percentage Error (MAPE) + pub mean_absolute_percentage_error: f64, + /// Root Mean Square Error (RMSE) + pub root_mean_square_error: f64, + /// Mean of the fitted distribution + pub distribution_mean: f64, + /// Median of the fitted distribution + pub distribution_median: f64, + /// Mode of the fitted distribution + pub distribution_mode: f64, + /// Proportion of data within 50% confidence interval (interquartile range) + pub within_50_percent: f64, + /// Proportion of data within 68% confidence interval (approx. 1 std dev) + pub within_68_percent: f64, + /// Proportion of data within 95% confidence interval (approx. 2 std dev) + pub within_95_percent: f64, + /// The number of data points used in the calculation + pub sample_size: usize, +} diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..463dab8 --- /dev/null +++ b/src/lib.rs @@ -0,0 +1,74 @@ +//! # Distimate: Estimation Distributions for Project Planning and Risk Analysis +//! +//! `distimate` is a Rust crate that provides a set of probability distributions +//! commonly used in estimation, along with a trait for working with these +//! distributions in an estimation context. It's designed to assist in project +//! planning, risk analysis, and decision-making under uncertainty. +//! +//! ## Features +//! +//! - A variety of estimation-focused probability distributions +//! - A unified `EstimationDistribution` trait for working with these distributions +//! - Tools for calculating confidence intervals, risk assessments, and fit quality +//! - Easy-to-use prelude module for convenient imports +//! +//! ## Quick Start +//! +//! To get started with `distimate`, add it to your `Cargo.toml`: +//! +//! ```toml +//! [dependencies] +//! distimate = "0.1.0" # replace with the current version +//! ``` +//! +//! Then, in your Rust file: +//! +//! ```rust +//! use distimate::prelude::*; +//! use distimate::Pert; +//! +//! fn main() -> Result<()> { +//! let pert = Pert::new(1.0, 2.0, 3.0)?; +//! println!("Expected value: {}", pert.expected_value()); +//! println!("95% confidence interval: {:?}", pert.confidence_interval(0.95)); +//! Ok(()) +//! } +//! ``` +//! +//! ## Available Distributions +//! +//! - PERT (Program Evaluation and Review Technique) +//! - Triangular +//! - Normal +//! +//! ## The EstimationDistribution Trait +//! +//! The `EstimationDistribution` trait provides a unified interface for working +//! with probability distributions in an estimation context. It includes methods +//! for calculating percentiles, confidence intervals, and assessing risk. +//! +//! ## Error Handling +//! +//! This crate uses a custom `Error` type and `Result` alias for error handling. +//! Most methods that can fail will return a `Result`. +//! +//! ## Prelude +//! +//! For convenience, commonly used items are re-exported in the `prelude` module. +//! You can import all of these at once with `use distimate::prelude::*;`. + +// Module declarations +mod distributions; +mod error; +mod estimation_distribution; +mod fit_quality; +pub mod prelude; +mod result; + +// Re-exports +pub use error::Error; +pub use estimation_distribution::EstimationDistribution; +pub use fit_quality::EstimationFitQuality; +pub use result::Result; + +pub use distributions::{Normal, Pert, Triangular}; diff --git a/src/prelude.rs b/src/prelude.rs new file mode 100644 index 0000000..f128e78 --- /dev/null +++ b/src/prelude.rs @@ -0,0 +1,8 @@ +pub use crate::{Error, EstimationDistribution, EstimationFitQuality, Result}; + +pub use rand::distributions::Distribution as RandDistribution; + +pub use statrs::{ + distribution::{Beta, Continuous, ContinuousCDF}, + statistics::{Distribution, Max, Median, Min, Mode}, +}; diff --git a/src/result.rs b/src/result.rs new file mode 100644 index 0000000..b465ded --- /dev/null +++ b/src/result.rs @@ -0,0 +1,3 @@ +use crate::Error; + +pub type Result = std::result::Result;