statrs/distribution/
fisher_snedecor.rs

1use crate::distribution::{Continuous, ContinuousCDF};
2use crate::function::beta;
3use crate::statistics::*;
4use crate::{Result, StatsError};
5use rand::Rng;
6use std::f64;
7
8/// Implements the
9/// [Fisher-Snedecor](https://en.wikipedia.org/wiki/F-distribution) distribution
10/// also commonly known as the F-distribution
11///
12/// # Examples
13///
14/// ```
15/// use statrs::distribution::{FisherSnedecor, Continuous};
16/// use statrs::statistics::Distribution;
17/// use statrs::prec;
18///
19/// let n = FisherSnedecor::new(3.0, 3.0).unwrap();
20/// assert_eq!(n.mean().unwrap(), 3.0);
21/// assert!(prec::almost_eq(n.pdf(1.0), 0.318309886183790671538, 1e-15));
22/// ```
23#[derive(Debug, Copy, Clone, PartialEq)]
24pub struct FisherSnedecor {
25    freedom_1: f64,
26    freedom_2: f64,
27}
28
29impl FisherSnedecor {
30    /// Constructs a new fisher-snedecor distribution with
31    /// degrees of freedom `freedom_1` and `freedom_2`
32    ///
33    /// # Errors
34    ///
35    /// Returns an error if `freedom_1` or `freedom_2` are `NaN`.
36    /// Also returns an error if `freedom_1 <= 0.0` or `freedom_2 <= 0.0`
37    ///
38    /// # Examples
39    ///
40    /// ```
41    /// use statrs::distribution::FisherSnedecor;
42    ///
43    /// let mut result = FisherSnedecor::new(1.0, 1.0);
44    /// assert!(result.is_ok());
45    ///
46    /// result = FisherSnedecor::new(0.0, 0.0);
47    /// assert!(result.is_err());
48    /// ```
49    pub fn new(freedom_1: f64, freedom_2: f64) -> Result<FisherSnedecor> {
50        if !freedom_1.is_finite() || freedom_1 <= 0.0 || !freedom_2.is_finite() || freedom_2 <= 0.0
51        {
52            Err(StatsError::BadParams)
53        } else {
54            Ok(FisherSnedecor {
55                freedom_1,
56                freedom_2,
57            })
58        }
59    }
60
61    /// Returns the first degree of freedom for the
62    /// fisher-snedecor distribution
63    ///
64    /// # Examples
65    ///
66    /// ```
67    /// use statrs::distribution::FisherSnedecor;
68    ///
69    /// let n = FisherSnedecor::new(2.0, 3.0).unwrap();
70    /// assert_eq!(n.freedom_1(), 2.0);
71    /// ```
72    pub fn freedom_1(&self) -> f64 {
73        self.freedom_1
74    }
75
76    /// Returns the second degree of freedom for the
77    /// fisher-snedecor distribution
78    ///
79    /// # Examples
80    ///
81    /// ```
82    /// use statrs::distribution::FisherSnedecor;
83    ///
84    /// let n = FisherSnedecor::new(2.0, 3.0).unwrap();
85    /// assert_eq!(n.freedom_2(), 3.0);
86    /// ```
87    pub fn freedom_2(&self) -> f64 {
88        self.freedom_2
89    }
90}
91
92impl ::rand::distributions::Distribution<f64> for FisherSnedecor {
93    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
94        (super::gamma::sample_unchecked(rng, self.freedom_1 / 2.0, 0.5) * self.freedom_2)
95            / (super::gamma::sample_unchecked(rng, self.freedom_2 / 2.0, 0.5) * self.freedom_1)
96    }
97}
98
99impl ContinuousCDF<f64, f64> for FisherSnedecor {
100    /// Calculates the cumulative distribution function for the fisher-snedecor
101    /// distribution
102    /// at `x`
103    ///
104    /// # Formula
105    ///
106    /// ```ignore
107    /// I_((d1 * x) / (d1 * x + d2))(d1 / 2, d2 / 2)
108    /// ```
109    ///
110    /// where `d1` is the first degree of freedom, `d2` is
111    /// the second degree of freedom, and `I` is the regularized incomplete
112    /// beta function
113    fn cdf(&self, x: f64) -> f64 {
114        if x < 0.0 {
115            0.0
116        } else if x.is_infinite() {
117            1.0
118        } else {
119            beta::beta_reg(
120                self.freedom_1 / 2.0,
121                self.freedom_2 / 2.0,
122                self.freedom_1 * x / (self.freedom_1 * x + self.freedom_2),
123            )
124        }
125    }
126
127    /// Calculates the survival function for the fisher-snedecor
128    /// distribution at `x`
129    ///
130    /// # Formula
131    ///
132    /// ```ignore
133    /// I_(1 - ((d1 * x) / (d1 * x + d2))(d2 / 2, d1 / 2)
134    /// ```
135    ///
136    /// where `d1` is the first degree of freedom, `d2` is
137    /// the second degree of freedom, and `I` is the regularized incomplete
138    /// beta function
139    fn sf(&self, x: f64) -> f64 {
140        if x < 0.0 {
141            1.0
142        } else if x.is_infinite() {
143            0.0
144        } else {
145            beta::beta_reg(
146                self.freedom_2 / 2.0,
147                self.freedom_1 / 2.0, 
148                1. - ((self.freedom_1 * x) / (self.freedom_1 * x + self.freedom_2))
149            )
150        }
151    }
152}
153
154impl Min<f64> for FisherSnedecor {
155    /// Returns the minimum value in the domain of the
156    /// fisher-snedecor distribution representable by a double precision
157    /// float
158    ///
159    /// # Formula
160    ///
161    /// ```ignore
162    /// 0
163    /// ```
164    fn min(&self) -> f64 {
165        0.0
166    }
167}
168
169impl Max<f64> for FisherSnedecor {
170    /// Returns the maximum value in the domain of the
171    /// fisher-snedecor distribution representable by a double precision
172    /// float
173    ///
174    /// # Formula
175    ///
176    /// ```ignore
177    /// INF
178    /// ```
179    fn max(&self) -> f64 {
180        f64::INFINITY
181    }
182}
183
184impl Distribution<f64> for FisherSnedecor {
185    /// Returns the mean of the fisher-snedecor distribution
186    ///
187    /// # Panics
188    ///
189    /// If `freedom_2 <= 2.0`
190    ///
191    /// # Remarks
192    ///
193    /// Returns `NaN` if `freedom_2` is `INF`
194    ///
195    /// # Formula
196    ///
197    /// ```ignore
198    /// d2 / (d2 - 2)
199    /// ```
200    ///
201    /// where `d2` is the second degree of freedom
202    fn mean(&self) -> Option<f64> {
203        if self.freedom_2 <= 2.0 {
204            None
205        } else {
206            Some(self.freedom_2 / (self.freedom_2 - 2.0))
207        }
208    }
209    /// Returns the variance of the fisher-snedecor distribution
210    ///
211    /// # Panics
212    ///
213    /// If `freedom_2 <= 4.0`
214    ///
215    /// # Remarks
216    ///
217    /// Returns `NaN` if `freedom_1` or `freedom_2` is `INF`
218    ///
219    /// # Formula
220    ///
221    /// ```ignore
222    /// (2 * d2^2 * (d1 + d2 - 2)) / (d1 * (d2 - 2)^2 * (d2 - 4))
223    /// ```
224    ///
225    /// where `d1` is the first degree of freedom and `d2` is
226    /// the second degree of freedom
227    fn variance(&self) -> Option<f64> {
228        if self.freedom_2 <= 4.0 {
229            None
230        } else {
231            let val =
232                (2.0 * self.freedom_2 * self.freedom_2 * (self.freedom_1 + self.freedom_2 - 2.0))
233                    / (self.freedom_1
234                        * (self.freedom_2 - 2.0)
235                        * (self.freedom_2 - 2.0)
236                        * (self.freedom_2 - 4.0));
237            Some(val)
238        }
239    }
240    /// Returns the skewness of the fisher-snedecor distribution
241    ///
242    /// # Panics
243    ///
244    /// If `freedom_2 <= 6.0`
245    ///
246    /// # Remarks
247    ///
248    /// Returns `NaN` if `freedom_1` or `freedom_2` is `INF`
249    ///
250    /// # Formula
251    ///
252    /// ```ignore
253    /// ((2d1 + d2 - 2) * sqrt(8 * (d2 - 4))) / ((d2 - 6) * sqrt(d1 * (d1 + d2
254    /// - 2)))
255    /// ```
256    ///
257    /// where `d1` is the first degree of freedom and `d2` is
258    /// the second degree of freedom
259    fn skewness(&self) -> Option<f64> {
260        if self.freedom_2 <= 6.0 {
261            None
262        } else {
263            let val = ((2.0 * self.freedom_1 + self.freedom_2 - 2.0)
264                * (8.0 * (self.freedom_2 - 4.0)).sqrt())
265                / ((self.freedom_2 - 6.0)
266                    * (self.freedom_1 * (self.freedom_1 + self.freedom_2 - 2.0)).sqrt());
267            Some(val)
268        }
269    }
270}
271
272impl Mode<Option<f64>> for FisherSnedecor {
273    /// Returns the mode for the fisher-snedecor distribution
274    ///
275    /// # Panics
276    ///
277    /// If `freedom_1 <= 2.0`
278    ///
279    /// # Remarks
280    ///
281    /// Returns `NaN` if `freedom_1` or `freedom_2` is `INF`
282    ///
283    /// # Formula
284    ///
285    /// ```ignore
286    /// ((d1 - 2) / d1) * (d2 / (d2 + 2))
287    /// ```
288    ///
289    /// where `d1` is the first degree of freedom and `d2` is
290    /// the second degree of freedom
291    fn mode(&self) -> Option<f64> {
292        if self.freedom_1 <= 2.0 {
293            None
294        } else {
295            let val = (self.freedom_2 * (self.freedom_1 - 2.0))
296                / (self.freedom_1 * (self.freedom_2 + 2.0));
297            Some(val)
298        }
299    }
300}
301
302impl Continuous<f64, f64> for FisherSnedecor {
303    /// Calculates the probability density function for the fisher-snedecor
304    /// distribution
305    /// at `x`
306    ///
307    /// # Remarks
308    ///
309    /// Returns `NaN` if `freedom_1`, `freedom_2` is `INF`, or `x` is `+INF` or
310    /// `-INF`
311    ///
312    /// # Formula
313    ///
314    /// ```ignore
315    /// sqrt(((d1 * x) ^ d1 * d2 ^ d2) / (d1 * x + d2) ^ (d1 + d2)) / (x * β(d1
316    /// / 2, d2 / 2))
317    /// ```
318    ///
319    /// where `d1` is the first degree of freedom, `d2` is
320    /// the second degree of freedom, and `β` is the beta function
321    fn pdf(&self, x: f64) -> f64 {
322        if x.is_infinite() || x <= 0.0 {
323            0.0
324        } else {
325            ((self.freedom_1 * x).powf(self.freedom_1) * self.freedom_2.powf(self.freedom_2)
326                / (self.freedom_1 * x + self.freedom_2).powf(self.freedom_1 + self.freedom_2))
327            .sqrt()
328                / (x * beta::beta(self.freedom_1 / 2.0, self.freedom_2 / 2.0))
329        }
330    }
331
332    /// Calculates the log probability density function for the fisher-snedecor
333    /// distribution
334    /// at `x`
335    ///
336    /// # Remarks
337    ///
338    /// Returns `NaN` if `freedom_1`, `freedom_2` is `INF`, or `x` is `+INF` or
339    /// `-INF`
340    ///
341    /// # Formula
342    ///
343    /// ```ignore
344    /// ln(sqrt(((d1 * x) ^ d1 * d2 ^ d2) / (d1 * x + d2) ^ (d1 + d2)) / (x *
345    /// β(d1 / 2, d2 / 2)))
346    /// ```
347    ///
348    /// where `d1` is the first degree of freedom, `d2` is
349    /// the second degree of freedom, and `β` is the beta function
350    fn ln_pdf(&self, x: f64) -> f64 {
351        self.pdf(x).ln()
352    }
353}
354
355#[rustfmt::skip]
356#[cfg(all(test, feature = "nightly"))]
357mod tests {
358    use crate::statistics::*;
359    use crate::distribution::{ContinuousCDF, Continuous, FisherSnedecor};
360    use crate::distribution::internal::*;
361    use crate::consts::ACC;
362
363    fn try_create(freedom_1: f64, freedom_2: f64) -> FisherSnedecor {
364        let n = FisherSnedecor::new(freedom_1, freedom_2);
365        assert!(n.is_ok());
366        n.unwrap()
367    }
368
369    fn create_case(freedom_1: f64, freedom_2: f64) {
370        let n = try_create(freedom_1, freedom_2);
371        assert_eq!(freedom_1, n.freedom_1());
372        assert_eq!(freedom_2, n.freedom_2());
373    }
374
375    fn bad_create_case(freedom_1: f64, freedom_2: f64) {
376        let n = FisherSnedecor::new(freedom_1, freedom_2);
377        assert!(n.is_err());
378    }
379
380     fn get_value<F>(freedom_1: f64, freedom_2: f64, eval: F) -> f64
381        where F: Fn(FisherSnedecor) -> f64
382    {
383        let n = try_create(freedom_1, freedom_2);
384        eval(n)
385    }
386
387    fn test_case<F>(freedom_1: f64, freedom_2: f64, expected: f64, eval: F)
388        where F: Fn(FisherSnedecor) -> f64
389    {
390        let x = get_value(freedom_1, freedom_2, eval);
391        assert_eq!(expected, x);
392    }
393
394    fn test_almost<F>(freedom_1: f64, freedom_2: f64, expected: f64, acc: f64, eval: F)
395        where F: Fn(FisherSnedecor) -> f64
396    {
397        let x = get_value(freedom_1, freedom_2, eval);
398        assert_almost_eq!(expected, x, acc);
399    }
400
401    #[test]
402    fn test_create() {
403        create_case(0.1, 0.1);
404        create_case(1.0, 0.1);
405        create_case(10.0, 0.1);
406        create_case(0.1, 1.0);
407        create_case(1.0, 1.0);
408        create_case(10.0, 1.0);
409    }
410
411    #[test]
412    fn test_bad_create() {
413        bad_create_case(f64::NAN, f64::NAN);
414        bad_create_case(0.0, f64::NAN);
415        bad_create_case(-1.0, f64::NAN);
416        bad_create_case(-10.0, f64::NAN);
417        bad_create_case(f64::NAN, 0.0);
418        bad_create_case(0.0, 0.0);
419        bad_create_case(-1.0, 0.0);
420        bad_create_case(-10.0, 0.0);
421        bad_create_case(f64::NAN, -1.0);
422        bad_create_case(0.0, -1.0);
423        bad_create_case(-1.0, -1.0);
424        bad_create_case(-10.0, -1.0);
425        bad_create_case(f64::NAN, -10.0);
426        bad_create_case(0.0, -10.0);
427        bad_create_case(-1.0, -10.0);
428        bad_create_case(-10.0, -10.0);
429        bad_create_case(f64::INFINITY, 0.1);
430        bad_create_case(0.1, f64::INFINITY);
431        bad_create_case(f64::INFINITY, f64::INFINITY);
432    }
433
434    #[test]
435    fn test_mean() {
436        let mean = |x: FisherSnedecor| x.mean().unwrap();
437        test_case(0.1, 10.0, 1.25, mean);
438        test_case(1.0, 10.0, 1.25, mean);
439        test_case(10.0, 10.0, 1.25, mean);
440    }
441
442    #[test]
443    #[should_panic]
444    fn test_mean_with_low_d2() {
445        let mean = |x: FisherSnedecor| x.mean().unwrap();
446        get_value(0.1, 0.1, mean);
447    }
448
449    #[test]
450    fn test_variance() {
451        let variance = |x: FisherSnedecor| x.variance().unwrap();
452        test_almost(0.1, 10.0, 42.1875, 1e-14, variance);
453        test_case(1.0, 10.0, 4.6875, variance);
454        test_case(10.0, 10.0, 0.9375, variance);
455    }
456
457    #[test]
458    #[should_panic]
459    fn test_variance_with_low_d2() {
460        let variance = |x: FisherSnedecor| x.variance().unwrap();
461        get_value(0.1, 0.1, variance);
462    }
463
464    #[test]
465    fn test_skewness() {
466        let skewness = |x: FisherSnedecor| x.skewness().unwrap();
467        test_almost(0.1, 10.0, 15.78090735784977089658, 1e-14, skewness);
468        test_case(1.0, 10.0, 5.773502691896257645091, skewness);
469        test_case(10.0, 10.0, 3.614784456460255759501, skewness);
470    }
471
472    #[test]
473    #[should_panic]
474    fn test_skewness_with_low_d2() {
475        let skewness = |x: FisherSnedecor| x.skewness().unwrap();
476        get_value(0.1, 0.1, skewness);
477    }
478
479    #[test]
480    fn test_mode() {
481        let mode = |x: FisherSnedecor| x.mode().unwrap();
482        test_case(10.0, 0.1, 0.0380952380952380952381, mode);
483        test_case(10.0, 1.0, 4.0 / 15.0, mode);
484        test_case(10.0, 10.0, 2.0 / 3.0, mode);
485    }
486
487    #[test]
488    #[should_panic]
489    fn test_mode_with_low_d1() {
490        let mode = |x: FisherSnedecor| x.mode().unwrap();
491        get_value(0.1, 0.1, mode);
492    }
493
494    #[test]
495    fn test_min_max() {
496        let min = |x: FisherSnedecor| x.min();
497        let max = |x: FisherSnedecor| x.max();
498        test_case(1.0, 1.0, 0.0, min);
499        test_case(1.0, 1.0, f64::INFINITY, max);
500    }
501
502    #[test]
503    fn test_pdf() {
504        let pdf = |arg: f64| move |x: FisherSnedecor| x.pdf(arg);
505        test_almost(0.1, 0.1, 0.0234154207226588982471, 1e-16, pdf(1.0));
506        test_almost(1.0, 0.1, 0.0396064560910663979961, 1e-16, pdf(1.0));
507        test_almost(10.0, 0.1, 0.0418440630400545297349, 1e-14, pdf(1.0));
508        test_almost(0.1, 1.0, 0.0396064560910663979961, 1e-16, pdf(1.0));
509        test_almost(1.0, 1.0, 0.1591549430918953357689, 1e-16, pdf(1.0));
510        test_almost(10.0, 1.0, 0.230361989229138647108, 1e-16, pdf(1.0));
511        test_almost(0.1, 0.1, 0.00221546909694001013517, 1e-18, pdf(10.0));
512        test_almost(1.0, 0.1, 0.00369960370387922619592, 1e-17, pdf(10.0));
513        test_almost(10.0, 0.1, 0.00390179721174142927402, 1e-15, pdf(10.0));
514        test_almost(0.1, 1.0, 0.00319864073359931548273, 1e-17, pdf(10.0));
515        test_almost(1.0, 1.0, 0.009150765837179460915678, 1e-17, pdf(10.0));
516        test_almost(10.0, 1.0, 0.0116493859171442148446, 1e-17, pdf(10.0));
517        test_almost(0.1, 10.0, 0.00305087016058573989694, 1e-15, pdf(10.0));
518        test_almost(1.0, 10.0, 0.00271897749113479577864, 1e-17, pdf(10.0));
519        test_almost(10.0, 10.0, 2.4289227234060500084E-4, 1e-18, pdf(10.0));
520    }
521
522    #[test]
523    fn test_ln_pdf() {
524        let ln_pdf = |arg: f64| move |x: FisherSnedecor| x.ln_pdf(arg);
525        test_almost(0.1, 0.1, 0.0234154207226588982471f64.ln(), 1e-15, ln_pdf(1.0));
526        test_almost(1.0, 0.1, 0.0396064560910663979961f64.ln(), 1e-15, ln_pdf(1.0));
527        test_almost(10.0, 0.1, 0.0418440630400545297349f64.ln(), 1e-13, ln_pdf(1.0));
528        test_almost(0.1, 1.0, 0.0396064560910663979961f64.ln(), 1e-15, ln_pdf(1.0));
529        test_almost(1.0, 1.0, 0.1591549430918953357689f64.ln(), 1e-15, ln_pdf(1.0));
530        test_almost(10.0, 1.0, 0.230361989229138647108f64.ln(), 1e-15, ln_pdf(1.0));
531        test_case(0.1, 0.1, 0.00221546909694001013517f64.ln(), ln_pdf(10.0));
532        test_almost(1.0, 0.1, 0.00369960370387922619592f64.ln(), 1e-15, ln_pdf(10.0));
533        test_almost(10.0, 0.1, 0.00390179721174142927402f64.ln(), 1e-13, ln_pdf(10.0));
534        test_almost(0.1, 1.0, 0.00319864073359931548273f64.ln(), 1e-15, ln_pdf(10.0));
535        test_almost(1.0, 1.0, 0.009150765837179460915678f64.ln(), 1e-15, ln_pdf(10.0));
536        test_case(10.0, 1.0, 0.0116493859171442148446f64.ln(), ln_pdf(10.0));
537        test_almost(0.1, 10.0, 0.00305087016058573989694f64.ln(), 1e-13, ln_pdf(10.0));
538        test_case(1.0, 10.0, 0.00271897749113479577864f64.ln(), ln_pdf(10.0));
539        test_almost(10.0, 10.0, 2.4289227234060500084E-4f64.ln(), 1e-14, ln_pdf(10.0));
540    }
541
542    #[test]
543    fn test_cdf() {
544        let cdf = |arg: f64| move |x: FisherSnedecor| x.cdf(arg);
545        test_almost(0.1, 0.1, 0.44712986033425140335, 1e-15, cdf(0.1));
546        test_almost(1.0, 0.1, 0.08156522095104674015, 1e-15, cdf(0.1));
547        test_almost(10.0, 0.1, 0.033184005716276536322, 1e-13, cdf(0.1));
548        test_almost(0.1, 1.0, 0.74378710917986379989, 1e-15, cdf(0.1));
549        test_almost(1.0, 1.0, 0.1949822290421366451595, 1e-16, cdf(0.1));
550        test_almost(10.0, 1.0, 0.0101195597354337146205, 1e-17, cdf(0.1));
551        test_almost(0.1, 0.1, 0.5, 1e-15, cdf(1.0));
552        test_almost(1.0, 0.1, 0.16734351500944271141, 1e-14, cdf(1.0));
553        test_almost(10.0, 0.1, 0.12207560664741704938, 1e-13, cdf(1.0));
554        test_almost(0.1, 1.0, 0.83265648499055728859, 1e-15, cdf(1.0));
555        test_almost(1.0, 1.0, 0.5, 1e-15, cdf(1.0));
556        test_almost(10.0, 1.0, 0.340893132302059872675, 1e-15, cdf(1.0));
557    }
558
559    #[test]
560    fn test_cdf_lower_bound() {
561        let cdf = |arg: f64| move |x: FisherSnedecor| x.cdf(arg);
562        test_case(0.1, 0.1, 0.0, cdf(-1.0));
563    }
564
565    #[test]
566    fn test_sf() {
567        let sf = |arg: f64| move |x: FisherSnedecor| x.sf(arg);
568        test_almost(0.1, 0.1, 0.5528701396657489, 1e-12, sf(0.1));
569        test_almost(1.0, 0.1, 0.9184347790489533, 1e-12, sf(0.1));
570        test_almost(10.0, 0.1, 0.9668159942836896, 1e-12, sf(0.1));
571        test_almost(0.1, 1.0, 0.25621289082013654, 1e-12, sf(0.1));
572        test_almost(1.0, 1.0, 0.8050177709578634, 1e-12, sf(0.1));
573        test_almost(10.0, 1.0, 0.9898804402645662, 1e-12, sf(0.1));
574        test_almost(0.1, 0.1, 0.5, 1e-15, sf(1.0));
575        test_almost(1.0, 0.1, 0.8326564849905562, 1e-12, sf(1.0));
576        test_almost(10.0, 0.1, 0.8779243933525519, 1e-12, sf(1.0));
577        test_almost(0.1, 1.0, 0.16734351500944344, 1e-12, sf(1.0));
578        test_almost(1.0, 1.0, 0.5, 1e-12, sf(1.0));
579        test_almost(10.0, 1.0, 0.65910686769794, 1e-12, sf(1.0));
580    }
581
582    #[test]
583    fn test_sf_lower_bound() {
584        let sf = |arg: f64| move |x: FisherSnedecor| x.sf(arg);
585        test_case(0.1, 0.1, 1.0, sf(-1.0));
586    }
587
588    #[test]
589    fn test_continuous() {
590        test::check_continuous_distribution(&try_create(10.0, 10.0), 0.0, 10.0);
591    }
592}
statrs/distribution/fisher_snedecor.rs

statrs/distribution/
fisher_snedecor.rs
