r2rs_nmath/distribution/f/
mod.rs

1mod d;
2mod non_central;
3mod p;
4mod q;
5mod r;
6
7use strafe_type::{LogProbability64, Positive64, Probability64, Rational64, Real64};
8
9pub(crate) use self::{d::*, non_central::*, p::*, q::*, r::*};
10use crate::traits::{Distribution, RNG};
11
12/// # The F Distribution
13///
14/// ## Description:
15///
16/// Density, distribution function, quantile function and random
17/// generation for the F distribution with ‘df1’ and ‘df2’ degrees of
18/// freedom (and optional non-centrality parameter ‘ncp’).
19///
20/// ## Arguments:
21///
22/// * df1, df2: degrees of freedom.  ‘Inf’ is allowed.
23/// * ncp: non-centrality parameter. If omitted the central F is
24/// assumed.
25///
26/// ## Details:
27///
28/// The F distribution with ‘df1 =’ n1 and ‘df2 =’ n2 degrees of
29/// freedom has density
30///
31/// $f(x) = \frac{\Gamma(\frac{n1 + n2}{2})}{\Gamma(\frac{n1}{2}) \Gamma(\frac{n2}{2})}
32/// (\frac{n1}{n2})^{\frac{n1}{2}} x^{\frac{n1}{2} - 1}
33/// (1 + \frac{n1}{n2} x)^{-\frac{n1 + n2}{2}}$
34///
35/// for $x > 0$.
36///
37/// It is the distribution of the ratio of the mean squares of n1 and
38/// n2 independent standard normals, and hence of the ratio of two
39/// independent chi-squared variates each divided by its degrees of
40/// freedom.  Since the ratio of a normal and the root mean-square of
41/// m independent normals has a Student's t_m distribution, the square
42/// of a t_m variate has a F distribution on 1 and m degrees of
43/// freedom.
44///
45/// The non-central F distribution is again the ratio of mean squares
46/// of independent normals of unit variance, but those in the
47/// numerator are allowed to have non-zero means and ‘ncp’ is the sum
48/// of squares of the means.  See Chisquare for further details on
49/// non-central distributions.
50///
51/// ## Value:
52///
53/// ‘df’ gives the density, ‘pf’ gives the distribution function ‘qf’
54/// gives the quantile function, and ‘rf’ generates random deviates.
55///
56/// Invalid arguments will result in return value ‘NaN’, with a
57/// warning.
58///
59/// The length of the result is determined by ‘n’ for ‘rf’, and is the
60/// maximum of the lengths of the numerical arguments for the other
61/// functions.
62///
63/// The numerical arguments other than ‘n’ are recycled to the length
64/// of the result.  Only the first elements of the logical arguments
65/// are used.
66///
67/// ## Density Plot
68///
69/// ```rust
70/// # use r2rs_base::traits::StatisticalSlice;
71/// # use r2rs_nmath::{distribution::FBuilder, traits::Distribution};
72/// # use strafe_plot::prelude::{IntoDrawingArea, Line, Plot, PlotOptions, SVGBackend, BLACK};
73/// # use strafe_type::FloatConstraint;
74/// let f = FBuilder::new().build();
75/// let x = <[f64]>::sequence(-0.5, 4.0, 1000);
76/// let y = x
77///     .iter()
78///     .map(|x| f.density(x).unwrap())
79///     .collect::<Vec<_>>();
80///
81/// let root = SVGBackend::new("density.svg", (1024, 768)).into_drawing_area();
82/// Plot::new()
83///     .with_options(PlotOptions {
84///         x_axis_label: "x".to_string(),
85///         y_axis_label: "density".to_string(),
86///         ..Default::default()
87///     })
88///     .with_plottable(Line {
89///         x,
90///         y,
91///         color: BLACK,
92///         ..Default::default()
93///     })
94///     .plot(&root)
95///     .unwrap();
96/// # use std::fs::rename;
97/// #     drop(root);
98/// #     rename(
99/// #             format!("density.svg"),
100/// #             format!("src/distribution/f/doctest_out/density.svg"),
101/// #     )
102/// #     .unwrap();
103/// ```
104#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = embed_doc_image::embed_image!("density", "src/distribution/f/doctest_out/density.svg")))]
105#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = "![Density][density]"))]
106///
107/// ## Note:
108///
109/// Supplying ‘ncp = 0’ uses the algorithm for the non-central
110/// distribution, which is not the same algorithm used if ‘ncp’ is
111/// omitted.  This is to give consistent behaviour in extreme cases
112/// with values of ‘ncp’ very near zero.
113///
114/// The code for non-zero ‘ncp’ is principally intended to be used for
115/// moderate values of ‘ncp’: it will not be highly accurate,
116/// especially in the tails, for large values.
117///
118/// ## Source:
119///
120/// For the central case of ‘df’, computed _via_ a binomial
121/// probability, code contributed by Catherine Loader (see ‘dbinom’);
122/// for the non-central case computed _via_ ‘dbeta’, code contributed
123/// by Peter Ruckdeschel.
124///
125/// For ‘pf’, _via_ ‘pbeta’ (or for large ‘df2’, _via_ ‘pchisq’).
126///
127/// For ‘qf’, _via_ ‘qchisq’ for large ‘df2’, else _via_ ‘qbeta’.
128///
129/// ## References:
130///
131/// Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) _The New S
132/// Language_.  Wadsworth & Brooks/Cole.
133///
134/// Johnson, N. L., Kotz, S. and Balakrishnan, N. (1995) _Continuous
135/// Univariate Distributions_, volume 2, chapters 27 and 30.  Wiley,
136/// New York.
137///
138/// ## See Also:
139///
140/// Distributions for other standard distributions, including ‘dchisq’
141/// for chi-squared and ‘dt’ for Student's t distributions.
142///
143/// ## Examples:
144///
145/// Lower Tails
146/// ```rust
147/// # use r2rs_base::traits::StatisticalSlice;
148/// # use r2rs_nmath::{
149/// #     distribution::{FBuilder, TBuilder},
150/// #     traits::Distribution,
151/// # };
152/// # use strafe_type::FloatConstraint;
153/// let x = <[f64]>::sequence(0.001, 5.0, 100);
154/// let nu = 4.0;
155///
156/// let t = TBuilder::new().with_df(nu).unwrap().build();
157/// let f = FBuilder::new().with_df1(1).with_df2(nu).build();
158/// let r1 = x
159///     .iter()
160///     .map(|x| 2.0 * t.probability(x, true).unwrap() - 1.0)
161///     .collect::<Vec<_>>();
162/// let r2 = x
163///     .iter()
164///     .map(|x| f.probability(x.powi(2), true).unwrap())
165///     .collect::<Vec<_>>();
166/// println!("{r1:?}");
167/// println!("{r2:?}");
168/// # use std::{fs::File, io::Write};
169/// # let mut f = File::create("src/distribution/f/doctest_out/lower_tail.md").unwrap();
170/// # writeln!(f, "```output").unwrap();
171/// # writeln!(f, "{r1:?}").unwrap();
172/// # writeln!(f, "{r2:?}").unwrap();
173/// # writeln!(f, "```").unwrap();
174/// ```
175#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = include_str!("doctest_out/lower_tail.md")))]
176///
177/// Upper Tails
178/// ```rust
179/// # use r2rs_base::traits::StatisticalSlice;
180/// # use r2rs_nmath::{
181/// #     distribution::{FBuilder, TBuilder},
182/// #     traits::Distribution,
183/// # };
184/// # use strafe_type::FloatConstraint;
185/// let x = <[f64]>::sequence(0.001, 5.0, 100);
186/// let nu = 4.0;
187///
188/// let t = TBuilder::new().with_df(nu).unwrap().build();
189/// let f = FBuilder::new().with_df1(1).with_df2(nu).build();
190/// let r1 = x
191///     .iter()
192///     .map(|x| 2.0 * t.probability(x, false).unwrap())
193///     .collect::<Vec<_>>();
194/// let r2 = x
195///     .iter()
196///     .map(|x| f.probability(x.powi(2), false).unwrap())
197///     .collect::<Vec<_>>();
198/// println!("{r1:?}");
199/// println!("{r2:?}");
200/// # use std::{fs::File, io::Write};
201/// # let mut f = File::create("src/distribution/f/doctest_out/upper_tail.md").unwrap();
202/// # writeln!(f, "```output").unwrap();
203/// # writeln!(f, "{r1:?}").unwrap();
204/// # writeln!(f, "{r2:?}").unwrap();
205/// # writeln!(f, "```").unwrap();
206/// ```
207#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = include_str!("doctest_out/upper_tail.md")))]
208///
209/// The density of the square of a $t_m$ is $2*\frac{\{t-density}(x, m)}{2*x}$
210///
211/// Check this is the same as the density of $F_{1,m}$
212/// ```rust
213/// # use r2rs_base::traits::StatisticalSlice;
214/// # use r2rs_nmath::{
215/// #     distribution::{FBuilder, TBuilder},
216/// #     traits::Distribution,
217/// # };
218/// # use strafe_type::FloatConstraint;
219/// let x = <[f64]>::sequence(0.001, 5.0, 100);
220/// let nu = 5.0;
221///
222/// let t = TBuilder::new().with_df(nu).unwrap().build();
223/// let f = FBuilder::new().with_df1(1).with_df2(nu).build();
224/// let r1 = x
225///     .iter()
226///     .map(|x| t.density(x).unwrap() / x)
227///     .collect::<Vec<_>>();
228/// let r2 = x
229///     .iter()
230///     .map(|x| f.density(x.powi(2)).unwrap())
231///     .collect::<Vec<_>>();
232/// println!("{r1:?}");
233/// println!("{r2:?}");
234/// # use std::{fs::File, io::Write};
235/// # let mut f = File::create("src/distribution/f/doctest_out/dens.md").unwrap();
236/// # writeln!(f, "```output").unwrap();
237/// # writeln!(f, "{r1:?}").unwrap();
238/// # writeln!(f, "{r2:?}").unwrap();
239/// # writeln!(f, "```").unwrap();
240/// ```
241#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = include_str!("doctest_out/dens.md")))]
242///
243/// Identity:  qf(2*p - 1, 1, df) == qt(p, df)^2  for  p >= 1/2
244/// ```rust
245/// use r2rs_base::traits::{QuantileType, StatisticalSlice};
246/// use r2rs_nmath::{
247///     distribution::{FBuilder, TBuilder},
248///     traits::Distribution,
249/// };
250/// use strafe_type::FloatConstraint;
251///
252/// let p = <[f64]>::sequence(0.5, 0.99, 50);
253/// let df = 10.0;
254///
255/// let rel_err = |x: &[f64], y: &[f64]| {
256///     let mean = x
257///         .iter()
258///         .chain(y.iter())
259///         .map(|f| f.abs())
260///         .collect::<Vec<_>>()
261///         .mean();
262///     x.iter()
263///         .zip(y.iter())
264///         .map(|(x, y)| if x == y { 0.0 } else { (x - y).abs() / mean })
265///         .collect::<Vec<_>>()
266/// };
267///
268/// let t = TBuilder::new().with_df(df).unwrap().build();
269/// let f = FBuilder::new().with_df1(1).with_df2(df).build();
270/// let r1 = p
271///     .iter()
272///     .map(|p| t.quantile(p, true).unwrap().powi(2))
273///     .collect::<Vec<_>>();
274/// let r2 = p
275///     .iter()
276///     .map(|p| f.quantile(2.0 * p - 1.0, true).unwrap())
277///     .collect::<Vec<_>>();
278///
279/// let error = rel_err(&r1, &r2).quantile(&[0.9], QuantileType::S);
280///
281/// println!("{error:?}");
282/// # use std::{fs::File, io::Write};
283/// # let mut f = File::create("src/distribution/f/doctest_out/quan.md").unwrap();
284/// # writeln!(f, "```output").unwrap();
285/// # writeln!(f, "{error:?}").unwrap();
286/// # writeln!(f, "```").unwrap();
287/// ```
288#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = include_str!("doctest_out/quan.md")))]
289pub struct F {
290    df1: Rational64,
291    df2: Rational64,
292    ncp: Option<Positive64>,
293}
294
295impl Distribution for F {
296    fn density<R: Into<Real64>>(&self, x: R) -> Real64 {
297        if let Some(ncp) = self.ncp {
298            dnf(x, self.df1, self.df2, ncp, false)
299        } else {
300            df(x, self.df1, self.df2, false)
301        }
302    }
303
304    fn log_density<R: Into<Real64>>(&self, x: R) -> Real64 {
305        if let Some(ncp) = self.ncp {
306            dnf(x, self.df1, self.df2, ncp, true)
307        } else {
308            df(x, self.df1, self.df2, true)
309        }
310    }
311
312    fn probability<R: Into<Real64>>(&self, q: R, lower_tail: bool) -> Probability64 {
313        if let Some(ncp) = self.ncp {
314            pnf(q, self.df1, self.df2, ncp, lower_tail)
315        } else {
316            pf(q, self.df1, self.df2, lower_tail)
317        }
318    }
319
320    fn log_probability<R: Into<Real64>>(&self, q: R, lower_tail: bool) -> LogProbability64 {
321        if let Some(ncp) = self.ncp {
322            log_pnf(q, self.df1, self.df2, ncp, lower_tail)
323        } else {
324            log_pf(q, self.df1, self.df2, lower_tail)
325        }
326    }
327
328    fn quantile<P: Into<Probability64>>(&self, p: P, lower_tail: bool) -> Real64 {
329        if let Some(ncp) = self.ncp {
330            qnf(p, self.df1, self.df2, ncp, lower_tail)
331        } else {
332            qf(p, self.df1, self.df2, lower_tail)
333        }
334    }
335
336    fn log_quantile<LP: Into<LogProbability64>>(&self, p: LP, lower_tail: bool) -> Real64 {
337        if let Some(ncp) = self.ncp {
338            log_qnf(p, self.df1, self.df2, ncp, lower_tail)
339        } else {
340            log_qf(p, self.df1, self.df2, lower_tail)
341        }
342    }
343
344    fn random_sample<R: RNG>(&self, rng: &mut R) -> Real64 {
345        if let Some(ncp) = self.ncp {
346            rnf(self.df1, self.df2, ncp, rng)
347        } else {
348            rf(self.df1, self.df2, rng)
349        }
350    }
351}
352
353pub struct FBuilder {
354    df1: Option<Rational64>,
355    df2: Option<Rational64>,
356    ncp: Option<Positive64>,
357}
358
359impl FBuilder {
360    pub fn new() -> Self {
361        Self {
362            df1: None,
363            df2: None,
364            ncp: None,
365        }
366    }
367
368    pub fn with_df1<R: Into<Rational64>>(&mut self, df1: R) -> &mut Self {
369        self.df1 = Some(df1.into());
370        self
371    }
372
373    pub fn with_df2<R: Into<Rational64>>(&mut self, df2: R) -> &mut Self {
374        self.df2 = Some(df2.into());
375        self
376    }
377
378    pub fn with_ncp<P: Into<Positive64>>(&mut self, ncp: P) -> &mut Self {
379        self.ncp = Some(ncp.into());
380        self
381    }
382
383    pub fn build(&self) -> F {
384        let df1 = self.df1.unwrap_or(1.0.into());
385        let df2 = self.df2.unwrap_or(1.0.into());
386
387        F {
388            df1,
389            df2,
390            ncp: self.ncp,
391        }
392    }
393}
394
395#[cfg(test)]
396mod tests;
397
398#[cfg(all(test, feature = "enable_proptest"))]
399mod proptests;
400
401#[cfg(all(test, feature = "enable_covtest"))]
402mod covtests;
r2rs_nmath/distribution/f/mod.rs

r2rs_nmath/distribution/f/
mod.rs
