r2rs-nmath 0.1.1

mod d;
mod p;
mod q;
mod r;

use strafe_type::{LogProbability64, Probability64, Rational64, Real64};

pub(crate) use self::{d::*, p::*, q::*, r::*};
use crate::traits::{Distribution, RNG};

/// # The Weibull Distribution
///
/// ## Description
///
/// Density, distribution function, quantile function and random generation for the Weibull
/// distribution with parameters shape and scale.
///
/// ## Arguments
///
/// * shape, scale: shape and scale parameters, the latter defaulting to 1.
///
/// ## Details
///
/// The Weibull distribution with shape parameter a and scale parameter b has density given by
///
/// $ f(x) = (a/b) (x/b)^{a-1} exp(- (x/b)^a) $
///
/// for x > 0. The cumulative distribution function is $F(x) = 1 - exp(- (x/b)^a)$ on x > 0, the mean
/// is $E(X) = b \Gamma(1 + 1/a)$, and the $Var(X) = b^2 * (\Gamma(1 + 2/a) - (\Gamma(1 + 1/a))^2)$.
///
/// ## Density Plot
///
/// ```rust
/// # use r2rs_base::traits::StatisticalSlice;
/// # use r2rs_nmath::{distribution::WeibullBuilder, traits::Distribution};
/// # use strafe_plot::prelude::{IntoDrawingArea, Line, Plot, PlotOptions, SVGBackend, BLACK};
/// # use strafe_type::FloatConstraint;
/// let weibull = WeibullBuilder::new().build();
/// let x = <[f64]>::sequence(-1.0, 5.0, 1000);
/// let y = x
///     .iter()
///     .map(|x| weibull.density(x).unwrap())
///     .collect::<Vec<_>>();
///
/// let root = SVGBackend::new("density.svg", (1024, 768)).into_drawing_area();
/// Plot::new()
///     .with_options(PlotOptions {
///         x_axis_label: "x".to_string(),
///         y_axis_label: "density".to_string(),
///         ..Default::default()
///     })
///     .with_plottable(Line {
///         x,
///         y,
///         color: BLACK,
///         ..Default::default()
///     })
///     .plot(&root)
///     .unwrap();
/// # use std::fs::rename;
/// #     drop(root);
/// #     rename(
/// #             format!("density.svg"),
/// #             format!("src/distribution/weibull/doctest_out/density.svg"),
/// #     )
/// #     .unwrap();
/// ```
#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = embed_doc_image::embed_image!("density", "src/distribution/weibull/doctest_out/density.svg")))]
#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = "![Density][density]"))]
///
/// ## Note
///
/// The cumulative hazard $H(t) = - log(1 - F(t))$ is
///
/// -pweibull(t, a, b, lower = FALSE, log = TRUE)
/// which is just $H(t) = (t/b)^a$.
///
/// ## Source
///
/// \[dpq\]weibull are calculated directly from the definitions. rweibull uses inversion.
///
/// References
/// Johnson, N. L., Kotz, S. and Balakrishnan, N. (1995) Continuous Univariate Distributions,
/// volume 1, chapter 21. Wiley, New York.
///
/// ## See Also
///
/// Distributions for other standard distributions, including the Exponential which is a special
/// case of the Weibull distribution.
///
/// ## Examples
///
/// Using this x for everything below
/// ```rust
/// # use num_traits::FloatConst;
/// # use r2rs_nmath::{
/// #     distribution::{ExponentialBuilder, LogNormalBuilder, WeibullBuilder},
/// #     rng::MersenneTwister,
/// #     traits::{Distribution, RNG},
/// # };
/// # use strafe_type::FloatConstraint;
/// let mut rng = MersenneTwister::new();
/// rng.set_seed(1);
///
/// let lnorm = LogNormalBuilder::new().build();
///
/// let mut x = (0..50)
///     .map(|_| lnorm.random_sample(&mut rng).unwrap())
///     .collect::<Vec<_>>();
/// x.insert(0, 0.0);
/// ```
///
/// ```rust
/// # use num_traits::FloatConst;
/// # use r2rs_nmath::{
/// #     distribution::{ExponentialBuilder, LogNormalBuilder, WeibullBuilder},
/// #     rng::MersenneTwister,
/// #     traits::{Distribution, RNG},
/// # };
/// # use strafe_type::FloatConstraint;
/// # let mut rng = MersenneTwister::new();
/// # rng.set_seed(1);
/// #
/// # let lnorm = LogNormalBuilder::new().build();
/// #
/// # let mut x = (0..50)
/// #     .map(|_| lnorm.random_sample(&mut rng).unwrap())
/// #     .collect::<Vec<_>>();
/// # x.insert(0, 0.0);
/// let weibull = WeibullBuilder::new().with_shape(1).build();
/// let r1 = x
///     .iter()
///     .map(|x| weibull.density(x).unwrap())
///     .collect::<Vec<_>>();
///
/// let exp = ExponentialBuilder::new().build();
/// let r2 = x
///     .iter()
///     .map(|x| exp.density(x).unwrap())
///     .collect::<Vec<_>>();
///
/// println!("{r1:?}");
/// println!("{r2:?}");
/// # use std::{fs::File, io::Write};
/// # let mut f = File::create("src/distribution/weibull/doctest_out/cmp1.md").unwrap();
/// # writeln!(f, "```output").unwrap();
/// # writeln!(f, "{r1:?}").unwrap();
/// # writeln!(f, "{r1:?}").unwrap();
/// # writeln!(f, "```").unwrap();
/// ```
#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = include_str!("doctest_out/cmp1.md")))]
///
/// ```rust
/// # use num_traits::FloatConst;
/// # use r2rs_nmath::{
/// #     distribution::{ExponentialBuilder, LogNormalBuilder, WeibullBuilder},
/// #     rng::MersenneTwister,
/// #     traits::{Distribution, RNG},
/// # };
/// # use strafe_type::FloatConstraint;
/// # let mut rng = MersenneTwister::new();
/// # rng.set_seed(1);
/// #
/// # let lnorm = LogNormalBuilder::new().build();
/// #
/// # let mut x = (0..50)
/// #     .map(|_| lnorm.random_sample(&mut rng).unwrap())
/// #     .collect::<Vec<_>>();
/// # x.insert(0, 0.0);
/// let weibull = WeibullBuilder::new()
///     .with_shape(1)
///     .with_scale(f64::PI())
///     .build();
/// let r1 = x
///     .iter()
///     .map(|x| weibull.probability(x, true).unwrap())
///     .collect::<Vec<_>>();
///
/// let exp = ExponentialBuilder::new().with_rate(1.0 / f64::PI()).build();
/// let r2 = x
///     .iter()
///     .map(|x| exp.probability(x, true).unwrap())
///     .collect::<Vec<_>>();
///
/// println!("{r1:?}");
/// println!("{r2:?}");
/// # use std::{fs::File, io::Write};
/// # let mut f = File::create("src/distribution/weibull/doctest_out/cmp2.md").unwrap();
/// # writeln!(f, "```output").unwrap();
/// # writeln!(f, "{r1:?}").unwrap();
/// # writeln!(f, "{r1:?}").unwrap();
/// # writeln!(f, "```").unwrap();
/// ```
#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = include_str!("doctest_out/cmp2.md")))]
///
/// Cumulative hazard H()
/// ```rust
/// # use num_traits::FloatConst;
/// # use r2rs_nmath::{
/// #     distribution::{ExponentialBuilder, LogNormalBuilder, WeibullBuilder},
/// #     rng::MersenneTwister,
/// #     traits::{Distribution, RNG},
/// # };
/// # use strafe_type::FloatConstraint;
/// # let mut rng = MersenneTwister::new();
/// # rng.set_seed(1);
/// #
/// # let lnorm = LogNormalBuilder::new().build();
/// #
/// # let mut x = (0..50)
/// #     .map(|_| lnorm.random_sample(&mut rng).unwrap())
/// #     .collect::<Vec<_>>();
/// # x.insert(0, 0.0);
/// let weibull = WeibullBuilder::new()
///     .with_shape(2.5)
///     .with_scale(f64::PI())
///     .build();
/// let r1 = x
///     .iter()
///     .map(|x| weibull.log_probability(x, false).unwrap())
///     .collect::<Vec<_>>();
///
/// let r2 = x
///     .iter()
///     .map(|x| -(x / f64::PI()).powf(2.5))
///     .collect::<Vec<_>>();
///
/// println!("{r1:?}");
/// println!("{r2:?}");
/// # use std::{fs::File, io::Write};
/// # let mut f = File::create("src/distribution/weibull/doctest_out/cmp3.md").unwrap();
/// # writeln!(f, "```output").unwrap();
/// # writeln!(f, "{r1:?}").unwrap();
/// # writeln!(f, "{r1:?}").unwrap();
/// # writeln!(f, "```").unwrap();
/// ```
#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = include_str!("doctest_out/cmp3.md")))]
///
/// ```rust
/// # use num_traits::FloatConst;
/// # use r2rs_nmath::{
/// #     distribution::{ExponentialBuilder, LogNormalBuilder, WeibullBuilder},
/// #     rng::MersenneTwister,
/// #     traits::{Distribution, RNG},
/// # };
/// # use strafe_type::FloatConstraint;
/// # let mut rng = MersenneTwister::new();
/// # rng.set_seed(1);
/// #
/// # let lnorm = LogNormalBuilder::new().build();
/// #
/// # let mut x = (0..50)
/// #     .map(|_| lnorm.random_sample(&mut rng).unwrap())
/// #     .collect::<Vec<_>>();
/// # x.insert(0, 0.0);
/// let weibull = WeibullBuilder::new()
///     .with_shape(1)
///     .with_scale(f64::PI())
///     .build();
/// let r1 = x
///     .iter()
///     .map(|x| weibull.quantile(x / 11.0, true).unwrap())
///     .collect::<Vec<_>>();
///
/// let exp = ExponentialBuilder::new().with_rate(1.0 / f64::PI()).build();
/// let r2 = x
///     .iter()
///     .map(|x| exp.quantile(x / 11.0, true).unwrap())
///     .collect::<Vec<_>>();
///
/// println!("{r1:?}");
/// println!("{r2:?}");
/// # use std::{fs::File, io::Write};
/// # let mut f = File::create("src/distribution/weibull/doctest_out/cmp4.md").unwrap();
/// # writeln!(f, "```output").unwrap();
/// # writeln!(f, "{r1:?}").unwrap();
/// # writeln!(f, "{r1:?}").unwrap();
/// # writeln!(f, "```").unwrap();
/// ```
#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = include_str!("doctest_out/cmp4.md")))]
pub struct Weibull {
    shape: Rational64,
    scale: Rational64,
}

impl Distribution for Weibull {
    fn density<R: Into<Real64>>(&self, x: R) -> Real64 {
        dweibull(x, self.shape, self.scale, false)
    }

    fn log_density<R: Into<Real64>>(&self, x: R) -> Real64 {
        dweibull(x, self.shape, self.scale, true)
    }

    fn probability<R: Into<Real64>>(&self, q: R, lower_tail: bool) -> Probability64 {
        pweibull(q, self.shape, self.scale, lower_tail)
    }

    fn log_probability<R: Into<Real64>>(&self, q: R, lower_tail: bool) -> LogProbability64 {
        log_pweibull(q, self.shape, self.scale, lower_tail)
    }

    fn quantile<P: Into<Probability64>>(&self, p: P, lower_tail: bool) -> Real64 {
        qweibull(p, self.shape, self.scale, lower_tail)
    }

    fn log_quantile<LP: Into<LogProbability64>>(&self, p: LP, lower_tail: bool) -> Real64 {
        log_qweibull(p, self.shape, self.scale, lower_tail)
    }

    fn random_sample<R: RNG>(&self, rng: &mut R) -> Real64 {
        rweibull(self.shape, self.scale, rng)
    }
}

pub struct WeibullBuilder {
    shape: Option<Rational64>,
    scale: Option<Rational64>,
}

impl WeibullBuilder {
    pub fn new() -> Self {
        Self {
            shape: None,
            scale: None,
        }
    }

    pub fn with_shape<R: Into<Rational64>>(&mut self, shape: R) -> &mut Self {
        self.shape = Some(shape.into());
        self
    }

    pub fn with_scale<R: Into<Rational64>>(&mut self, scale: R) -> &mut Self {
        self.scale = Some(scale.into());
        self
    }

    pub fn build(&self) -> Weibull {
        let shape = self.shape.unwrap_or(1.0.into());
        let scale = self.scale.unwrap_or(1.0.into());

        Weibull { shape, scale }
    }
}

#[cfg(test)]
mod tests;

#[cfg(all(test, feature = "enable_proptest"))]
mod proptests;

#[cfg(all(test, feature = "enable_covtest"))]
mod covtests;