r2rs-nmath 0.1.1

use strafe_type::{FloatConstraint, LogProbability64, Positive64, Probability64, Real64};

use crate::{
    distribution::tukey::{p::ptukey, q::qtukey},
    traits::{Distribution, RNG},
};

mod p;
mod q;

/// # The Studentized Range Distribution
///
/// ## Description
///
/// Functions of the distribution of the studentized range, $ R/s $, where $ R $ is the range of a
/// standard normal sample and $ df*s^2 $ is independently distributed as chi-squared with df
/// degrees of freedom, see pchisq.
///
/// ## Arguments
///
/// * nmeans: sample size for range (same for each group).
/// * df: degrees of freedom for s (see below).
/// * nranges: number of groups whose maximum range is considered.
///
/// ## Details
///
/// If ng = nranges is greater than one, R is the maximum of ng groups of nmeans observations each.
///
/// ## Note
///
/// A Legendre 16-point formula is used for the integral of ptukey. The computations are relatively
/// expensive, especially for qtukey which uses a simple secant method for finding the inverse of
/// ptukey. qtukey will be accurate to the 4th decimal place.
///
/// ## Source
///
/// qtukey is in part adapted from Odeh and Evans (1974).
///
/// ## References
///
/// Copenhaver, Margaret Diponzio and Holland, Burt S. (1988). Computation of the distribution of
/// the maximum studentized range statistic with application to multiple significance testing of
/// simple effects. Journal of Statistical Computation and Simulation, 30, 1–15.
/// doi: 10.1080/00949658808811082.
///
/// Odeh, R. E. and Evans, J. O. (1974). Algorithm AS 70: Percentage Points of the Normal
/// Distribution. Applied Statistics, 23, 96–97. doi: 10.2307/2347061.
///
/// ## See Also
///
/// Distributions for standard distributions, including pnorm and qnorm for the corresponding
/// functions for the normal distribution.
///
/// ## Examples
///
/// ```rust
/// # use r2rs_base::traits::StatisticalSlice;
/// # use r2rs_nmath::{distribution::TukeyBuilder, traits::Distribution};
/// # use strafe_plot::prelude::{IntoDrawingArea, Line, Plot, PlotOptions, SVGBackend, BLACK};
/// # use strafe_type::FloatConstraint;
/// let tukey = TukeyBuilder::new().with_means(6).with_df(5).build();
/// let x = <[f64]>::sequence(-1.0, 8.0, 1000);
/// let y = x
///     .iter()
///     .map(|x| tukey.probability(x, true).unwrap())
///     .collect::<Vec<_>>();
///
/// let root = SVGBackend::new("prob.svg", (1024, 768)).into_drawing_area();
/// Plot::new()
///     .with_options(PlotOptions {
///         x_axis_label: "x".to_string(),
///         y_axis_label: "probability".to_string(),
///         ..Default::default()
///     })
///     .with_plottable(Line {
///         x,
///         y,
///         color: BLACK,
///         ..Default::default()
///     })
///     .plot(&root)
///     .unwrap();
/// # use std::fs::rename;
/// #     drop(root);
/// #     rename(
/// #             format!("prob.svg"),
/// #             format!("src/distribution/tukey/doctest_out/prob.svg"),
/// #     )
/// #     .unwrap();
/// ```
#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = embed_doc_image::embed_image!("prob", "src/distribution/tukey/doctest_out/prob.svg")))]
#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = "![Probability][prob]"))]
///
/// The precision may be not much more than about 8 digits
/// ```rust
/// # use r2rs_nmath::{distribution::TukeyBuilder, traits::Distribution};
/// # use strafe_type::FloatConstraint;
/// let tukey = TukeyBuilder::new().with_means(2).with_df(5).build();
/// let x = (0..=10).collect::<Vec<_>>();
/// let p = x
///     .iter()
///     .map(|x| tukey.probability(x, true).unwrap())
///     .collect::<Vec<_>>();
/// println!("{p:?}");
///
/// let df = (2..=11).collect::<Vec<_>>();
/// let q = df
///     .iter()
///     .map(|df| {
///         TukeyBuilder::new()
///             .with_means(2)
///             .with_df(df)
///             .build()
///             .quantile(0.95, true)
///             .unwrap()
///     })
///     .collect::<Vec<_>>();
/// println!("{q:?}");
///
/// let err = q
///     .iter()
///     .zip(df.iter())
///     .map(|(q, df)| {
///         0.95 - TukeyBuilder::new()
///             .with_means(2)
///             .with_df(df)
///             .build()
///             .probability(q, true)
///             .unwrap()
///     })
///     .collect::<Vec<_>>();
/// println!("{err:?}");
/// # use std::{fs::File, io::Write};
/// # let mut f = File::create("src/distribution/tukey/doctest_out/err.md").unwrap();
/// # writeln!(f, "```output").unwrap();
/// # writeln!(f, "{p:?}").unwrap();
/// # writeln!(f, "{q:?}").unwrap();
/// # writeln!(f, "{err:?}").unwrap();
/// # writeln!(f, "```").unwrap();
/// ```
#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = include_str!("doctest_out/err.md")))]
pub struct Tukey {
    ranges: Positive64,
    means: Positive64,
    df: Positive64,
}

impl Distribution for Tukey {
    fn density<R: Into<Real64>>(&self, _x: R) -> Real64 {
        unimplemented!()
    }

    fn log_density<R: Into<Real64>>(&self, _x: R) -> Real64 {
        unimplemented!()
    }

    fn probability<R: Into<Real64>>(&self, q: R, lower_tail: bool) -> Probability64 {
        let ret = ptukey(q, self.ranges, self.means, self.df, lower_tail, false);
        ret.unwrap().into()
    }

    fn log_probability<R: Into<Real64>>(&self, q: R, lower_tail: bool) -> LogProbability64 {
        let ret = ptukey(q, self.ranges, self.means, self.df, lower_tail, true);
        ret.unwrap().into()
    }

    fn quantile<P: Into<Probability64>>(&self, p: P, lower_tail: bool) -> Real64 {
        qtukey(
            p.into().unwrap(),
            self.ranges,
            self.means,
            self.df,
            lower_tail,
            false,
        )
    }

    fn log_quantile<LP: Into<LogProbability64>>(&self, p: LP, lower_tail: bool) -> Real64 {
        qtukey(
            p.into().unwrap(),
            self.ranges,
            self.means,
            self.df,
            lower_tail,
            true,
        )
    }

    fn random_sample<R: RNG>(&self, _rng: &mut R) -> Real64 {
        unimplemented!()
    }
}

pub struct TukeyBuilder {
    ranges: Option<Positive64>,
    means: Option<Positive64>,
    df: Option<Positive64>,
}

impl TukeyBuilder {
    pub fn new() -> Self {
        Self {
            ranges: None,
            means: None,
            df: None,
        }
    }

    pub fn with_ranges<P: Into<Positive64>>(&mut self, ranges: P) -> &mut Self {
        self.ranges = Some(ranges.into());
        self
    }

    pub fn with_means<P: Into<Positive64>>(&mut self, means: P) -> &mut Self {
        self.means = Some(means.into());
        self
    }

    pub fn with_df<P: Into<Positive64>>(&mut self, df: P) -> &mut Self {
        self.df = Some(df.into());
        self
    }

    pub fn build(&self) -> Tukey {
        let ranges = self.ranges.unwrap_or(1.0.into());
        let means = self.means.unwrap_or(2.0.into());
        let df = self.df.unwrap_or(2.0.into());

        Tukey { ranges, means, df }
    }
}

#[cfg(test)]
mod tests;

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

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