russell_stat 1.13.0

use crate::{ProbabilityDistribution, StrError};
use rand::Rng;
use rand_distr::{Distribution, Frechet};
use russell_lab::math::gamma;

const FRECHET_MIN_DELTA_X: f64 = 1e-15;

/// Implements the Frechet distribution
///
/// This is a Type II Extreme Value Distribution (largest value)
///
/// See: <https://en.wikipedia.org/wiki/Fr%C3%A9chet_distribution>
///
/// ![Frechet](https://raw.githubusercontent.com/cpmech/russell/main/russell_stat/data/figures/plot_distribution_functions_frechet.svg)
pub struct DistributionFrechet {
    location: f64, // location of the minimum value
    scale: f64,    // scale parameter
    shape: f64,    // shape parameter

    sampler: Frechet<f64>, // sampler
}

impl DistributionFrechet {
    /// Allocates a new instance
    ///
    /// # Input
    ///
    /// * `location` -- location parameter
    /// * `shape` -- shape parameter
    pub fn new(location: f64, scale: f64, shape: f64) -> Result<Self, StrError> {
        Ok(DistributionFrechet {
            location,
            scale,
            shape,
            sampler: Frechet::new(location, scale, shape).map_err(|_| "invalid parameters")?,
        })
    }
}

impl ProbabilityDistribution for DistributionFrechet {
    /// Evaluates the Probability Density Function (PDF)
    ///
    /// # Examples
    ///
    /// ```
    /// use russell_lab::math::NAPIER;
    /// use russell_stat::*;
    ///
    /// fn main() -> Result<(), StrError> {
    ///     let (location, scale, shape) = (0.0, 1.0, 2.0);
    ///     let dist = DistributionFrechet::new(location, scale, shape)?;
    ///     assert_eq!(dist.pdf(1.0), 2.0 / NAPIER);
    ///     // Mathematica: PDF[FrechetDistribution[2, 1], 1] == 2/E
    ///     Ok(())
    /// }
    /// ```
    fn pdf(&self, x: f64) -> f64 {
        if x - self.location < FRECHET_MIN_DELTA_X {
            return 0.0;
        }
        let z = (x - self.location) / self.scale;
        f64::exp(-f64::powf(z, -self.shape)) * f64::powf(z, -1.0 - self.shape) * self.shape / self.scale
    }

    /// Evaluates the Cumulative Distribution Function (CDF)
    ///
    /// # Examples
    ///
    /// ```
    /// use russell_lab::approx_eq;
    /// use russell_stat::*;
    ///
    /// fn main() -> Result<(), StrError> {
    ///     // Example from <https://reference.wolfram.com/language/ref/FrechetDistribution.html>
    ///     let (location, scale, shape) = (0.0, 6.3, 0.71);
    ///     let dist = DistributionFrechet::new(location, scale, shape)?;
    ///     let p_xx_less_than_30 = dist.cdf(30.0);
    ///     let p_xx_greater_than_30 = 1.0 - p_xx_less_than_30;
    ///     println!("P(X > 30) = {}", p_xx_greater_than_30);
    ///     approx_eq(p_xx_greater_than_30, 0.2812192884182272, 1e-15);
    ///     Ok(())
    /// }
    /// ```
    fn cdf(&self, x: f64) -> f64 {
        if x - self.location < FRECHET_MIN_DELTA_X {
            return 0.0;
        }
        let z = (x - self.location) / self.scale;
        f64::exp(-f64::powf(z, -self.shape))
    }

    /// Returns the Mean
    ///
    /// # Examples
    ///
    /// ```
    /// use russell_lab::{approx_eq, math::SQRT_PI};
    /// use russell_stat::*;
    ///
    /// fn main() -> Result<(), StrError> {
    ///     let (location, scale, shape) = (0.0, 3.0, 2.0);
    ///     let dist = DistributionFrechet::new(location, scale, shape)?;
    ///     approx_eq(dist.mean(), 3.0 * SQRT_PI, 1e-15);
    ///     // Mathematica: Mean[FrechetDistribution[2, 3]]
    ///     Ok(())
    /// }
    /// ```
    fn mean(&self) -> f64 {
        if self.shape > 1.0 {
            return self.location + self.scale * gamma(1.0 - 1.0 / self.shape);
        }
        f64::INFINITY
    }

    /// Returns the Variance
    ///
    /// # Examples
    ///
    /// ```
    /// use russell_lab::approx_eq;
    /// use russell_stat::*;
    ///
    /// fn main() -> Result<(), StrError> {
    ///     let (location, scale, shape) = (0.0, 3.0, 4.0);
    ///     let dist = DistributionFrechet::new(location, scale, shape)?;
    ///     approx_eq(dist.variance(), 2.4372698060239768, 1e-14);
    ///     // Mathematica: N[Variance[FrechetDistribution[4, 3]], 17]
    ///     Ok(())
    /// }
    /// ```
    fn variance(&self) -> f64 {
        if self.shape > 2.0 {
            let ss = self.scale * self.scale;
            return ss * (gamma(1.0 - 2.0 / self.shape) - f64::powi(gamma(1.0 - 1.0 / self.shape), 2));
        }
        f64::INFINITY
    }

    /// Generates a pseudo-random number belonging to this probability distribution
    ///
    /// # Examples
    ///
    /// ```
    /// use russell_stat::*;
    ///
    /// fn main() -> Result<(), StrError> {
    ///     let (location, scale, shape) = (0.0, 3.0, 4.0);
    ///     let dist = DistributionFrechet::new(location, scale, shape)?;
    ///     let mut rng = get_rng();
    ///     println!("sample = {}", dist.sample(&mut rng));
    ///     Ok(())
    /// }
    /// ```
    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
        self.sampler.sample(rng)
    }
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#[cfg(test)]
mod tests {
    use crate::{get_rng, DistributionFrechet, ProbabilityDistribution};
    use russell_lab::approx_eq;

    // Data from the following R-code (run with Rscript frechet.R):
    /*
    # needs r-cran-evd
    library(evd)
    X <- seq(0, 4, 0.5)
    L <- c(0, 0.5)  # location
    C <- c(1, 2)    # scale
    A <- c(1, 2, 3) # shape
    Y <- matrix(ncol=5)
    first <- TRUE
    for (l in L) {
        for (c in C) {
            for (a in A) {
                pdf <- dfrechet(X, l, c, a)
                cdf <- pfrechet(X, l, c, a)
                for (i in 1:length(X)) {
                    if (first) {
                        Y <- rbind(c(X[i], l, c, a, pdf[i], cdf[i]))
                        first <- FALSE
                    } else {
                        Y <- rbind(Y, c(X[i], l, c, a, pdf[i], cdf[i]))
                    }
                }
            }
        }
    }
    write.table(format(Y, digits=15), "/tmp/frechet.dat", row.names=FALSE, col.names=c("x","location","scale","shape","pdf","cdf"), quote=FALSE)
    print("file </tmp/frechet.dat> written")
    */

    #[test]
    fn frechet_handles_errors() {
        assert_eq!(
            DistributionFrechet::new(2.0, 3.0, f64::INFINITY).err(),
            Some("invalid parameters")
        );
    }

    #[test]
    fn frechet_works() {
        #[rustfmt::skip]
        // x, location, scale, shape, pdf, cdf
        let data = [
            [0.00000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 1.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
            [5.00000000000000e-01, 0.00000000000000e+00, 1.00000000000000e+00, 1.00000000000000e+00, 5.41341132946451e-01, 1.35335283236613e-01],
            [1.00000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 1.00000000000000e+00, 3.67879441171442e-01, 3.67879441171442e-01],
            [1.50000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 1.00000000000000e+00, 2.28185386236708e-01, 5.13417119032592e-01],
            [2.00000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 1.00000000000000e+00, 1.51632664928158e-01, 6.06530659712633e-01],
            [2.50000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 1.00000000000000e+00, 1.07251207365702e-01, 6.70320046035639e-01],
            [3.00000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 1.00000000000000e+00, 7.96145900637543e-02, 7.16531310573789e-01],
            [3.50000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 1.00000000000000e+00, 6.13450851490029e-02, 7.51477293075286e-01],
            [4.00000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 1.00000000000000e+00, 4.86750489419628e-02, 7.78800783071405e-01],
            [0.00000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 2.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
            [5.00000000000000e-01, 0.00000000000000e+00, 1.00000000000000e+00, 2.00000000000000e+00, 2.93050222219747e-01, 1.83156388887342e-02],
            [1.00000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 2.00000000000000e+00, 7.35758882342885e-01, 3.67879441171442e-01],
            [1.50000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 2.00000000000000e+00, 3.79958748699232e-01, 6.41180388429955e-01],
            [2.00000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 2.00000000000000e+00, 1.94700195767851e-01, 7.78800783071405e-01],
            [2.50000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 2.00000000000000e+00, 1.09074404987675e-01, 8.52143788966211e-01],
            [3.00000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 2.00000000000000e+00, 6.62843938381015e-02, 8.94839316814370e-01],
            [3.50000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 2.00000000000000e+00, 4.29905748010600e-02, 9.21610447297725e-01],
            [4.00000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 2.00000000000000e+00, 2.93566582129211e-02, 9.39413062813476e-01],
            [0.00000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 3.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
            [5.00000000000000e-01, 0.00000000000000e+00, 1.00000000000000e+00, 3.00000000000000e+00, 1.61022061393206e-02, 3.35462627902512e-04],
            [1.00000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 3.00000000000000e+00, 1.10363832351433e+00, 3.67879441171442e-01],
            [1.50000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 3.00000000000000e+00, 4.40632343233130e-01, 7.43567079205906e-01],
            [2.00000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 3.00000000000000e+00, 1.65468169234612e-01, 8.82496902584595e-01],
            [2.50000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 3.00000000000000e+00, 7.20387839639600e-02, 9.38004999530729e-01],
            [3.00000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 3.00000000000000e+00, 3.56903868259736e-02, 9.63640444301286e-01],
            [3.50000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 3.00000000000000e+00, 1.95307877314814e-02, 9.76946277985140e-01],
            [4.00000000000000e+00, 0.00000000000000e+00, 1.00000000000000e+00, 3.00000000000000e+00, 1.15370676211571e-02, 9.84496437005408e-01],
            [0.00000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 1.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
            [5.00000000000000e-01, 0.00000000000000e+00, 2.00000000000000e+00, 1.00000000000000e+00, 1.46525111109873e-01, 1.83156388887342e-02],
            [1.00000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 1.00000000000000e+00, 2.70670566473225e-01, 1.35335283236613e-01],
            [1.50000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 1.00000000000000e+00, 2.34308567213979e-01, 2.63597138115727e-01],
            [2.00000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 1.00000000000000e+00, 1.83939720585721e-01, 3.67879441171442e-01],
            [2.50000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 1.00000000000000e+00, 1.43785268517511e-01, 4.49328964117222e-01],
            [3.00000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 1.00000000000000e+00, 1.14092693118354e-01, 5.13417119032592e-01],
            [3.50000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 1.00000000000000e+00, 9.21988770624913e-02, 5.64718122007759e-01],
            [4.00000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 1.00000000000000e+00, 7.58163324640792e-02, 6.06530659712633e-01],
            [0.00000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 2.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
            [5.00000000000000e-01, 0.00000000000000e+00, 2.00000000000000e+00, 2.00000000000000e+00, 7.20225118203259e-06, 1.12535174719259e-07],
            [1.00000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 2.00000000000000e+00, 1.46525111109873e-01, 1.83156388887342e-02],
            [1.50000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 2.00000000000000e+00, 4.00624155036601e-01, 1.69013315406066e-01],
            [2.00000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 2.00000000000000e+00, 3.67879441171442e-01, 3.67879441171442e-01],
            [2.50000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 2.00000000000000e+00, 2.69973721110041e-01, 5.27292424043049e-01],
            [3.00000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 2.00000000000000e+00, 1.89979374349616e-01, 6.41180388429955e-01],
            [3.50000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 2.00000000000000e+00, 1.34609406946660e-01, 7.21422290354756e-01],
            [4.00000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 2.00000000000000e+00, 9.73500978839256e-02, 7.78800783071405e-01],
            [0.00000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 3.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
            [5.00000000000000e-01, 0.00000000000000e+00, 2.00000000000000e+00, 3.00000000000000e+00, 6.15863381970675e-26, 1.60381089054864e-28],
            [1.00000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 3.00000000000000e+00, 8.05110306966028e-03, 3.35462627902512e-04],
            [1.50000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 3.00000000000000e+00, 4.43003781677631e-01, 9.34461101976254e-02],
            [2.00000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 3.00000000000000e+00, 5.51819161757164e-01, 3.67879441171442e-01],
            [2.50000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 3.00000000000000e+00, 3.68207332052299e-01, 5.99295787845538e-01],
            [3.00000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 3.00000000000000e+00, 2.20316171616565e-01, 7.43567079205906e-01],
            [3.50000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 3.00000000000000e+00, 1.32710267758176e-01, 8.29784773144222e-01],
            [4.00000000000000e+00, 0.00000000000000e+00, 2.00000000000000e+00, 3.00000000000000e+00, 8.27340846173058e-02, 8.82496902584595e-01],
            [0.00000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 1.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
            [5.00000000000000e-01, 5.00000000000000e-01, 1.00000000000000e+00, 1.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
            [1.00000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 1.00000000000000e+00, 5.41341132946451e-01, 1.35335283236613e-01],
            [1.50000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 1.00000000000000e+00, 3.67879441171442e-01, 3.67879441171442e-01],
            [2.00000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 1.00000000000000e+00, 2.28185386236708e-01, 5.13417119032592e-01],
            [2.50000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 1.00000000000000e+00, 1.51632664928158e-01, 6.06530659712633e-01],
            [3.00000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 1.00000000000000e+00, 1.07251207365702e-01, 6.70320046035639e-01],
            [3.50000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 1.00000000000000e+00, 7.96145900637543e-02, 7.16531310573789e-01],
            [4.00000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 1.00000000000000e+00, 6.13450851490029e-02, 7.51477293075286e-01],
            [0.00000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 2.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
            [5.00000000000000e-01, 5.00000000000000e-01, 1.00000000000000e+00, 2.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
            [1.00000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 2.00000000000000e+00, 2.93050222219747e-01, 1.83156388887342e-02],
            [1.50000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 2.00000000000000e+00, 7.35758882342885e-01, 3.67879441171442e-01],
            [2.00000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 2.00000000000000e+00, 3.79958748699232e-01, 6.41180388429955e-01],
            [2.50000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 2.00000000000000e+00, 1.94700195767851e-01, 7.78800783071405e-01],
            [3.00000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 2.00000000000000e+00, 1.09074404987675e-01, 8.52143788966211e-01],
            [3.50000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 2.00000000000000e+00, 6.62843938381015e-02, 8.94839316814370e-01],
            [4.00000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 2.00000000000000e+00, 4.29905748010600e-02, 9.21610447297725e-01],
            [0.00000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 3.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
            [5.00000000000000e-01, 5.00000000000000e-01, 1.00000000000000e+00, 3.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
            [1.00000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 3.00000000000000e+00, 1.61022061393206e-02, 3.35462627902512e-04],
            [1.50000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 3.00000000000000e+00, 1.10363832351433e+00, 3.67879441171442e-01],
            [2.00000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 3.00000000000000e+00, 4.40632343233130e-01, 7.43567079205906e-01],
            [2.50000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 3.00000000000000e+00, 1.65468169234612e-01, 8.82496902584595e-01],
            [3.00000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 3.00000000000000e+00, 7.20387839639600e-02, 9.38004999530729e-01],
            [3.50000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 3.00000000000000e+00, 3.56903868259736e-02, 9.63640444301286e-01],
            [4.00000000000000e+00, 5.00000000000000e-01, 1.00000000000000e+00, 3.00000000000000e+00, 1.95307877314814e-02, 9.76946277985140e-01],
            [0.00000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 1.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
            [5.00000000000000e-01, 5.00000000000000e-01, 2.00000000000000e+00, 1.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
            [1.00000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 1.00000000000000e+00, 1.46525111109873e-01, 1.83156388887342e-02],
            [1.50000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 1.00000000000000e+00, 2.70670566473225e-01, 1.35335283236613e-01],
            [2.00000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 1.00000000000000e+00, 2.34308567213979e-01, 2.63597138115727e-01],
            [2.50000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 1.00000000000000e+00, 1.83939720585721e-01, 3.67879441171442e-01],
            [3.00000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 1.00000000000000e+00, 1.43785268517511e-01, 4.49328964117222e-01],
            [3.50000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 1.00000000000000e+00, 1.14092693118354e-01, 5.13417119032592e-01],
            [4.00000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 1.00000000000000e+00, 9.21988770624913e-02, 5.64718122007759e-01],
            [0.00000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 2.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
            [5.00000000000000e-01, 5.00000000000000e-01, 2.00000000000000e+00, 2.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
            [1.00000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 2.00000000000000e+00, 7.20225118203259e-06, 1.12535174719259e-07],
            [1.50000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 2.00000000000000e+00, 1.46525111109873e-01, 1.83156388887342e-02],
            [2.00000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 2.00000000000000e+00, 4.00624155036601e-01, 1.69013315406066e-01],
            [2.50000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 2.00000000000000e+00, 3.67879441171442e-01, 3.67879441171442e-01],
            [3.00000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 2.00000000000000e+00, 2.69973721110041e-01, 5.27292424043049e-01],
            [3.50000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 2.00000000000000e+00, 1.89979374349616e-01, 6.41180388429955e-01],
            [4.00000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 2.00000000000000e+00, 1.34609406946660e-01, 7.21422290354756e-01],
            [0.00000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 3.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
            [5.00000000000000e-01, 5.00000000000000e-01, 2.00000000000000e+00, 3.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
            [1.00000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 3.00000000000000e+00, 6.15863381970675e-26, 1.60381089054864e-28],
            [1.50000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 3.00000000000000e+00, 8.05110306966028e-03, 3.35462627902512e-04],
            [2.00000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 3.00000000000000e+00, 4.43003781677631e-01, 9.34461101976254e-02],
            [2.50000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 3.00000000000000e+00, 5.51819161757164e-01, 3.67879441171442e-01],
            [3.00000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 3.00000000000000e+00, 3.68207332052299e-01, 5.99295787845538e-01],
            [3.50000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 3.00000000000000e+00, 2.20316171616565e-01, 7.43567079205906e-01],
            [4.00000000000000e+00, 5.00000000000000e-01, 2.00000000000000e+00, 3.00000000000000e+00, 1.32710267758176e-01, 8.29784773144222e-01],
        ];
        for row in data {
            let [x, location, scale, shape, pdf, cdf] = row;
            let d = DistributionFrechet::new(location, scale, shape).unwrap();
            approx_eq(d.pdf(x), pdf, 1e-14);
            approx_eq(d.cdf(x), cdf, 1e-14);
        }
    }

    #[test]
    fn mean_and_variance_work() {
        let location = 8.782275;
        let scale = 1.0;
        let shape = 4.095645;
        let d = DistributionFrechet::new(location, scale, shape).unwrap();
        approx_eq(d.mean(), 10.0, 1e-6);
        approx_eq(d.variance(), 0.25, 1e-6);

        let d = DistributionFrechet::new(location, scale, 1.0).unwrap();
        assert_eq!(d.mean(), f64::INFINITY);
        assert_eq!(d.variance(), f64::INFINITY);
    }

    #[test]
    fn sample_works() {
        let d = DistributionFrechet::new(1.0, 2.0, 3.0).unwrap();
        let mut rng = get_rng();
        d.sample(&mut rng);
    }
}