statrs/distribution/
chi.rs

1use crate::distribution::{Continuous, ContinuousCDF};
2use crate::function::gamma;
3use crate::statistics::*;
4use crate::{Result, StatsError};
5use rand::Rng;
6use std::f64;
7
8/// Implements the [Chi](https://en.wikipedia.org/wiki/Chi_distribution)
9/// distribution
10///
11/// # Examples
12///
13/// ```
14/// use statrs::distribution::{Chi, Continuous};
15/// use statrs::statistics::Distribution;
16/// use statrs::prec;
17///
18/// let n = Chi::new(2.0).unwrap();
19/// assert!(prec::almost_eq(n.mean().unwrap(), 1.25331413731550025121, 1e-14));
20/// assert!(prec::almost_eq(n.pdf(1.0), 0.60653065971263342360, 1e-15));
21/// ```
22#[derive(Debug, Copy, Clone, PartialEq)]
23pub struct Chi {
24    freedom: f64,
25}
26
27impl Chi {
28    /// Constructs a new chi distribution
29    /// with `freedom` degrees of freedom
30    ///
31    /// # Errors
32    ///
33    /// Returns an error if `freedom` is `NaN` or
34    /// less than or equal to `0.0`
35    ///
36    /// # Examples
37    ///
38    /// ```
39    /// use statrs::distribution::Chi;
40    ///
41    /// let mut result = Chi::new(2.0);
42    /// assert!(result.is_ok());
43    ///
44    /// result = Chi::new(0.0);
45    /// assert!(result.is_err());
46    /// ```
47    pub fn new(freedom: f64) -> Result<Chi> {
48        if freedom.is_nan() || freedom <= 0.0 {
49            Err(StatsError::BadParams)
50        } else {
51            Ok(Chi { freedom })
52        }
53    }
54
55    /// Returns the degrees of freedom of
56    /// the chi distribution.
57    ///
58    /// # Examples
59    ///
60    /// ```
61    /// use statrs::distribution::Chi;
62    ///
63    /// let n = Chi::new(2.0).unwrap();
64    /// assert_eq!(n.freedom(), 2.0);
65    /// ```
66    pub fn freedom(&self) -> f64 {
67        self.freedom
68    }
69}
70
71impl ::rand::distributions::Distribution<f64> for Chi {
72    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
73        (0..self.freedom as i64)
74            .fold(0.0, |acc, _| {
75                acc + super::normal::sample_unchecked(rng, 0.0, 1.0).powf(2.0)
76            })
77            .sqrt()
78    }
79}
80
81impl ContinuousCDF<f64, f64> for Chi {
82    /// Calculates the cumulative distribution function for the chi
83    /// distribution at `x`.
84    ///
85    /// # Formula
86    ///
87    /// ```ignore
88    /// P(k / 2, x^2 / 2)
89    /// ```
90    ///
91    /// where `k` is the degrees of freedom and `P` is
92    /// the regularized lower incomplete Gamma function
93    fn cdf(&self, x: f64) -> f64 {
94        if self.freedom == f64::INFINITY || x == f64::INFINITY {
95            1.0
96        } else if x <= 0.0 {
97            0.0
98        } else {
99            gamma::gamma_lr(self.freedom / 2.0, x * x / 2.0)
100        }
101    }
102
103    /// Calculates the survival function for the chi
104    /// distribution at `x`.
105    ///
106    /// # Formula
107    ///
108    /// ```ignore
109    /// P(k / 2, x^2 / 2)
110    /// ```
111    ///
112    /// where `k` is the degrees of freedom and `P` is
113    /// the regularized upper incomplete Gamma function
114    fn sf(&self, x: f64) -> f64 {
115        if self.freedom == f64::INFINITY || x == f64::INFINITY {
116            0.0
117        } else if x <= 0.0 {
118            1.0
119        } else {
120            gamma::gamma_ur(self.freedom / 2.0, x * x / 2.0)
121        }
122    }
123}
124
125impl Min<f64> for Chi {
126    /// Returns the minimum value in the domain of the chi distribution
127    /// representable by a double precision float
128    ///
129    /// # Formula
130    ///
131    /// ```ignore
132    /// 0
133    /// ```
134    fn min(&self) -> f64 {
135        0.0
136    }
137}
138
139impl Max<f64> for Chi {
140    /// Returns the maximum value in the domain of the chi distribution
141    /// representable by a double precision float
142    ///
143    /// # Formula
144    ///
145    /// ```ignore
146    /// INF
147    /// ```
148    fn max(&self) -> f64 {
149        f64::INFINITY
150    }
151}
152
153impl Distribution<f64> for Chi {
154    /// Returns the mean of the chi distribution
155    ///
156    /// # Remarks
157    ///
158    /// Returns `NaN` if `freedom` is `INF`
159    ///
160    /// # Formula
161    ///
162    /// ```ignore
163    /// sqrt2 * Γ((k + 1) / 2) / Γ(k / 2)
164    /// ```
165    ///
166    /// where `k` is degrees of freedom and `Γ` is the gamma function
167    fn mean(&self) -> Option<f64> {
168        if self.freedom.is_infinite() {
169            None
170        } else if self.freedom > 300.0 {
171            // Large n approximation based on the Stirling series approximation to the Gamma function
172            // This avoids call the Gamma function with large arguments and returning NaN
173            //
174            // Relative accuracy follows O(1/n^4) and at 300 d.o.f. is better than 1e-12
175            // For a f32 impl the threshold should be changed to 150
176            Some(
177                self.freedom.sqrt()
178                    / ((1.0 + 0.25 / self.freedom)
179                        * (1.0 + 0.03125 / (self.freedom * self.freedom))
180                        * (1.0 - 0.046875 / (self.freedom * self.freedom * self.freedom))),
181            )
182        } else {
183            let mean = f64::consts::SQRT_2 * gamma::gamma((self.freedom + 1.0) / 2.0)
184                / gamma::gamma(self.freedom / 2.0);
185            Some(mean)
186        }
187    }
188    /// Returns the variance of the chi distribution
189    ///
190    /// # Remarks
191    ///
192    /// Returns `NaN` if `freedom` is `INF`
193    ///
194    /// # Formula
195    ///
196    /// ```ignore
197    /// k - μ^2
198    /// ```
199    ///
200    /// where `k` is degrees of freedom and `μ` is the mean
201    /// of the distribution
202    fn variance(&self) -> Option<f64> {
203        let mean = self.mean()?;
204        Some(self.freedom - mean * mean)
205    }
206    /// Returns the entropy of the chi distribution
207    ///
208    /// # Remarks
209    ///
210    /// Returns `None` if `freedom` is `INF`
211    ///
212    /// # Formula
213    ///
214    /// ```ignore
215    /// ln(Γ(k / 2)) + 0.5 * (k - ln2 - (k - 1) * ψ(k / 2))
216    /// ```
217    ///
218    /// where `k` is degrees of freedom, `Γ` is the gamma function,
219    /// and `ψ` is the digamma function
220    fn entropy(&self) -> Option<f64> {
221        if self.freedom.is_infinite() {
222            return None;
223        }
224        let entr = gamma::ln_gamma(self.freedom / 2.0)
225            + (self.freedom
226                - (2.0f64).ln()
227                - (self.freedom - 1.0) * gamma::digamma(self.freedom / 2.0))
228                / 2.0;
229        Some(entr)
230    }
231    /// Returns the skewness of the chi distribution
232    ///
233    /// # Remarks
234    ///
235    /// Returns `NaN` if `freedom` is `INF`
236    ///
237    /// # Formula
238    ///
239    /// ```ignore
240    /// (μ / σ^3) * (1 - 2σ^2)
241    /// ```
242    /// where `μ` is the mean and `σ` the standard deviation
243    /// of the distribution
244    fn skewness(&self) -> Option<f64> {
245        let sigma = self.std_dev()?;
246        let skew = self.mean()? * (1.0 - 2.0 * sigma * sigma) / (sigma * sigma * sigma);
247        Some(skew)
248    }
249}
250
251impl Mode<Option<f64>> for Chi {
252    /// Returns the mode for the chi distribution
253    ///
254    /// # Panics
255    ///
256    /// If `freedom < 1.0`
257    ///
258    /// # Formula
259    ///
260    /// ```ignore
261    /// sqrt(k - 1)
262    /// ```
263    ///
264    /// where `k` is the degrees of freedom
265    fn mode(&self) -> Option<f64> {
266        if self.freedom - 1.0 < 0.0 {
267            return None;
268        }
269        Some((self.freedom - 1.0).sqrt())
270    }
271}
272
273impl Continuous<f64, f64> for Chi {
274    /// Calculates the probability density function for the chi
275    /// distribution at `x`
276    ///
277    /// # Formula
278    ///
279    /// ```ignore
280    /// (2^(1 - (k / 2)) * x^(k - 1) * e^(-x^2 / 2)) / Γ(k / 2)
281    /// ```
282    ///
283    /// where `k` is the degrees of freedom and `Γ` is the gamma function
284    fn pdf(&self, x: f64) -> f64 {
285        if self.freedom == f64::INFINITY || x == f64::INFINITY || x <= 0.0 {
286            0.0
287        } else if self.freedom > 160.0 {
288            self.ln_pdf(x).exp()
289        } else {
290            (2.0f64).powf(1.0 - self.freedom / 2.0)
291                * x.powf(self.freedom - 1.0)
292                * (-x * x / 2.0).exp()
293                / gamma::gamma(self.freedom / 2.0)
294        }
295    }
296
297    /// Calculates the log probability density function for the chi distribution
298    /// at `x`
299    ///
300    /// # Formula
301    ///
302    /// ```ignore
303    /// ln((2^(1 - (k / 2)) * x^(k - 1) * e^(-x^2 / 2)) / Γ(k / 2))
304    /// ```
305    fn ln_pdf(&self, x: f64) -> f64 {
306        if self.freedom == f64::INFINITY || x == f64::INFINITY || x <= 0.0 {
307            f64::NEG_INFINITY
308        } else {
309            (1.0 - self.freedom / 2.0) * (2.0f64).ln() + ((self.freedom - 1.0) * x.ln())
310                - x * x / 2.0
311                - gamma::ln_gamma(self.freedom / 2.0)
312        }
313    }
314}
315
316#[rustfmt::skip]
317#[cfg(all(test, feature = "nightly"))]
318mod tests {
319    use std::f64;
320    use crate::distribution::internal::*;
321    use crate::distribution::{Chi, Continuous, ContinuousCDF};
322    use crate::statistics::*;
323    use crate::consts::ACC;
324
325    fn try_create(freedom: f64) -> Chi {
326        let n = Chi::new(freedom);
327        assert!(n.is_ok());
328        n.unwrap()
329    }
330
331    fn create_case(freedom: f64) {
332        let n = try_create(freedom);
333        assert_eq!(freedom, n.freedom());
334    }
335
336    fn bad_create_case(freedom: f64) {
337        let n = Chi::new(freedom);
338        assert!(n.is_err());
339    }
340
341    fn get_value<F>(freedom: f64, eval: F) -> f64
342    where
343        F: Fn(Chi) -> f64,
344    {
345        let n = try_create(freedom);
346        eval(n)
347    }
348
349    fn test_case<F>(freedom: f64, expected: f64, eval: F)
350    where
351        F: Fn(Chi) -> f64,
352    {
353        let x = get_value(freedom, eval);
354        assert_eq!(expected, x);
355    }
356
357    fn test_almost<F>(freedom: f64, expected: f64, acc: f64, eval: F)
358    where
359        F: Fn(Chi) -> f64,
360    {
361        let x = get_value(freedom, eval);
362        assert_almost_eq!(expected, x, acc);
363    }
364
365    #[test]
366    fn test_create() {
367        create_case(1.0);
368        create_case(3.0);
369        create_case(f64::INFINITY);
370    }
371
372    #[test]
373    fn test_bad_create() {
374        bad_create_case(0.0);
375        bad_create_case(-1.0);
376        bad_create_case(-100.0);
377        bad_create_case(f64::NEG_INFINITY);
378        bad_create_case(f64::NAN);
379    }
380
381    #[test]
382    fn test_mean() {
383        let mean = |x: Chi| x.mean().unwrap();
384        test_almost(1.0, 0.7978845608028653558799, 1e-15, mean);
385        test_almost(2.0, 1.25331413731550025121, 1e-14, mean);
386        test_almost(2.5, 1.43396639245837498609, 1e-14, mean);
387        test_almost(5.0, 2.12769216214097428235, 1e-14, mean);
388        test_almost(336.0, 18.31666925443713, 1e-12, mean);
389    }
390
391    #[test]
392    fn test_large_dof_mean_not_nan() {
393        for i in 1..1000 {
394            let mean = Chi::new(i as f64).unwrap().mean().unwrap();
395            assert!(!mean.is_nan(), "Chi mean for {} dof was {}", i, mean);
396        }
397    }
398
399    #[test]
400    #[should_panic]
401    fn test_mean_degen() {
402        let mean = |x: Chi| x.mean().unwrap();
403        get_value(f64::INFINITY, mean);
404    }
405
406    #[test]
407    fn test_variance() {
408        let variance = |x: Chi| x.variance().unwrap();
409        test_almost(1.0, 0.3633802276324186569245, 1e-15, variance);
410        test_almost(2.0, 0.42920367320510338077, 1e-14, variance);
411        test_almost(2.5, 0.44374038529991368581, 1e-13, variance);
412        test_almost(3.0, 0.4535209105296746277, 1e-14, variance);
413    }
414
415    #[test]
416    #[should_panic]
417    fn test_variance_degen() {
418        let variance = |x: Chi| x.variance().unwrap();
419        get_value(f64::INFINITY, variance);
420    }
421
422    #[test]
423    fn test_entropy() {
424        let entropy = |x: Chi| x.entropy().unwrap();
425        test_almost(1.0, 0.7257913526447274323631, 1e-15, entropy);
426        test_almost(2.0, 0.9420342421707937755946, 1e-15, entropy);
427        test_almost(2.5, 0.97574472333041323989, 1e-14, entropy);
428        test_almost(3.0, 0.99615419810620560239, 1e-14, entropy);
429    }
430
431    #[test]
432    #[should_panic]
433    fn test_entropy_degen() {
434        let entropy = |x: Chi| x.entropy().unwrap();
435        get_value(f64::INFINITY, entropy);
436    }
437
438    #[test]
439    fn test_skewness() {
440        let skewness = |x: Chi| x.skewness().unwrap();
441        test_almost(1.0, 0.995271746431156042444, 1e-14, skewness);
442        test_almost(2.0, 0.6311106578189371382, 1e-13, skewness);
443        test_almost(2.5, 0.5458487096285153216, 1e-12, skewness);
444        test_almost(3.0, 0.485692828049590809, 1e-12, skewness);
445    }
446
447    #[test]
448    #[should_panic]
449    fn test_skewness_degen() {
450        let skewness = |x: Chi| x.skewness().unwrap();
451        get_value(f64::INFINITY, skewness);
452    }
453
454    #[test]
455    fn test_mode() {
456        let mode = |x: Chi| x.mode().unwrap();
457        test_case(1.0, 0.0, mode);
458        test_case(2.0, 1.0, mode);
459        test_case(2.5, 1.224744871391589049099, mode);
460        test_case(3.0, f64::consts::SQRT_2, mode);
461        test_case(f64::INFINITY, f64::INFINITY, mode);
462    }
463
464    #[test]
465    #[should_panic]
466    fn test_mode_freedom_lt_1() {
467        let mode = |x: Chi| x.mode().unwrap();
468        get_value(0.5, mode);
469    }
470
471    #[test]
472    fn test_min_max() {
473        let min = |x: Chi| x.min();
474        let max = |x: Chi| x.max();
475        test_case(1.0, 0.0, min);
476        test_case(2.0, 0.0, min);
477        test_case(2.5, 0.0, min);
478        test_case(3.0, 0.0, min);
479        test_case(f64::INFINITY, 0.0, min);
480        test_case(1.0, f64::INFINITY, max);
481        test_case(2.0, f64::INFINITY, max);
482        test_case(2.5, f64::INFINITY, max);
483        test_case(3.0, f64::INFINITY, max);
484        test_case(f64::INFINITY, f64::INFINITY, max);
485    }
486
487    #[test]
488    fn test_pdf() {
489        let pdf = |arg: f64| move |x: Chi| x.pdf(arg);
490        test_case(1.0, 0.0, pdf(0.0));
491        test_almost(1.0, 0.79390509495402353102, 1e-15, pdf(0.1));
492        test_almost(1.0, 0.48394144903828669960, 1e-15, pdf(1.0));
493        test_almost(1.0, 2.1539520085086552718e-7, 1e-22, pdf(5.5));
494        test_case(1.0, 0.0, pdf(f64::INFINITY));
495        test_case(2.0, 0.0, pdf(0.0));
496        test_almost(2.0, 0.099501247919268231335, 1e-16, pdf(0.1));
497        test_almost(2.0, 0.60653065971263342360, 1e-15, pdf(1.0));
498        test_almost(2.0, 1.4847681768496578863e-6, 1e-21, pdf(5.5));
499        test_case(2.0, 0.0, pdf(f64::INFINITY));
500        test_case(2.5, 0.0, pdf(0.0));
501        test_almost(2.5, 0.029191065334961657461, 1e-16, pdf(0.1));
502        test_almost(2.5, 0.56269645152636456261, 1e-15, pdf(1.0));
503        test_almost(2.5, 3.2304380188895211768e-6, 1e-20, pdf(5.5));
504        test_case(2.5, 0.0, pdf(f64::INFINITY));
505        test_case(f64::INFINITY, 0.0, pdf(0.0));
506        test_case(f64::INFINITY, 0.0, pdf(0.1));
507        test_case(f64::INFINITY, 0.0, pdf(1.0));
508        test_case(f64::INFINITY, 0.0, pdf(5.5));
509        test_case(f64::INFINITY, 0.0, pdf(f64::INFINITY));
510        test_almost(170.0, 0.5644678498668440878, 1e-13, pdf(13.0));
511    }
512
513    #[test]
514    fn test_neg_pdf() {
515        let pdf = |arg: f64| move |x: Chi| x.pdf(arg);
516        test_case(1.0, 0.0, pdf(-1.0));
517    }
518
519    #[test]
520    fn test_ln_pdf() {
521        let ln_pdf = |arg: f64| move |x: Chi| x.ln_pdf(arg);
522        test_case(1.0, f64::NEG_INFINITY, ln_pdf(0.0));
523        test_almost(1.0, -0.23079135264472743236, 1e-15, ln_pdf(0.1));
524        test_almost(1.0, -0.72579135264472743236, 1e-15, ln_pdf(1.0));
525        test_almost(1.0, -15.350791352644727432, 1e-14, ln_pdf(5.5));
526        test_case(1.0, f64::NEG_INFINITY, ln_pdf(f64::INFINITY));
527        test_case(2.0, f64::NEG_INFINITY, ln_pdf(0.0));
528        test_almost(2.0, -2.3075850929940456840, 1e-15, ln_pdf(0.1));
529        test_almost(2.0, -0.5, 1e-15, ln_pdf(1.0));
530        test_almost(2.0, -13.420251907761574765, 1e-15, ln_pdf(5.5));
531        test_case(2.0, f64::NEG_INFINITY, ln_pdf(f64::INFINITY));
532        test_case(2.5, f64::NEG_INFINITY, ln_pdf(0.0));
533        test_almost(2.5, -3.5338925982092416919, 1e-15, ln_pdf(0.1));
534        test_almost(2.5, -0.57501495871817316589, 1e-15, ln_pdf(1.0));
535        test_almost(2.5, -12.642892820360535314, 1e-16, ln_pdf(5.5));
536        test_case(2.5, f64::NEG_INFINITY, ln_pdf(f64::INFINITY));
537        test_case(f64::INFINITY, f64::NEG_INFINITY, ln_pdf(0.0));
538        test_case(f64::INFINITY, f64::NEG_INFINITY, ln_pdf(0.1));
539        test_case(f64::INFINITY, f64::NEG_INFINITY, ln_pdf(1.0));
540        test_case(f64::INFINITY, f64::NEG_INFINITY, ln_pdf(5.5));
541        test_case(f64::INFINITY, f64::NEG_INFINITY, ln_pdf(f64::INFINITY));
542        test_almost(170.0, -0.57187185030600516424237, 1e-13, ln_pdf(13.0));
543    }
544
545    #[test]
546    fn test_neg_ln_pdf() {
547        let ln_pdf = |arg: f64| move |x: Chi| x.ln_pdf(arg);
548        test_case(1.0, f64::NEG_INFINITY, ln_pdf(-1.0));
549    }
550
551    #[test]
552    fn test_cdf() {
553        let cdf = |arg: f64| move |x: Chi| x.cdf(arg);
554        test_case(1.0, 0.0, cdf(0.0));
555        test_almost(1.0, 0.079655674554057962931, 1e-16, cdf(0.1));
556        test_almost(1.0, 0.68268949213708589717, 1e-15, cdf(1.0));
557        test_case(1.0, 0.99999996202087506822, cdf(5.5));
558        test_case(1.0, 1.0, cdf(f64::INFINITY));
559        test_case(2.0, 0.0, cdf(0.0));
560        test_almost(2.0, 0.0049875208073176866474, 1e-17, cdf(0.1));
561        test_almost(2.0, 0.39346934028736657640, 1e-15, cdf(1.0));
562        test_case(2.0, 0.99999973004214966370, cdf(5.5));
563        test_case(2.0, 1.0, cdf(f64::INFINITY));
564        test_case(2.5, 0.0, cdf(0.0));
565        test_almost(2.5, 0.0011702413714030096290, 1e-18, cdf(0.1));
566        test_almost(2.5, 0.28378995266531297417, 1e-16, cdf(1.0));
567        test_case(2.5, 0.99999940337322804750, cdf(5.5));
568        test_case(2.5, 1.0, cdf(f64::INFINITY));
569        test_case(f64::INFINITY, 1.0, cdf(0.0));
570        test_case(f64::INFINITY, 1.0, cdf(0.1));
571        test_case(f64::INFINITY, 1.0, cdf(1.0));
572        test_case(f64::INFINITY, 1.0, cdf(5.5));
573        test_case(f64::INFINITY, 1.0, cdf(f64::INFINITY));
574    }
575
576    #[test]
577    fn test_sf() {
578        let sf = |arg: f64| move |x: Chi| x.sf(arg);
579        test_case(1.0, 1.0, sf(0.0));
580        test_almost(1.0, 0.920344325445942, 1e-16, sf(0.1));
581        test_almost(1.0, 0.31731050786291404, 1e-15, sf(1.0));
582        test_almost(1.0, 3.797912493177544e-8, 1e-15, sf(5.5));
583        test_case(1.0, 0.0, sf(f64::INFINITY));
584        test_case(2.0, 1.0, sf(0.0));
585        test_almost(2.0, 0.9950124791926823, 1e-17, sf(0.1));
586        test_almost(2.0, 0.6065306597126333, 1e-15, sf(1.0));
587        test_almost(2.0, 2.699578503363014e-7, 1e-15, sf(5.5));
588        test_case(2.0, 0.0, sf(f64::INFINITY));
589        test_case(2.5, 1.0, sf(0.0));
590        test_almost(2.5, 0.998829758628597, 1e-18, sf(0.1));
591        test_almost(2.5, 0.716210047334687, 1e-16, sf(1.0));
592        test_almost(2.5, 5.966267719870189e-7, 1e-15, sf(5.5));
593        test_case(2.5, 0.0, sf(f64::INFINITY));
594        test_case(f64::INFINITY, 0.0, sf(0.0));
595        test_case(f64::INFINITY, 0.0, sf(0.1));
596        test_case(f64::INFINITY, 0.0, sf(1.0));
597        test_case(f64::INFINITY, 0.0, sf(5.5));
598        test_case(f64::INFINITY, 0.0, sf(f64::INFINITY));
599    }
600
601    #[test]
602    fn test_neg_cdf() {
603        let cdf = |arg: f64| move |x: Chi| x.cdf(arg);
604        test_case(1.0, 0.0, cdf(-1.0));
605    }
606
607    #[test]
608    fn test_neg_sf() {
609        let sf = |arg: f64| move |x: Chi| x.sf(arg);
610        test_case(1.0, 1.0, sf(-1.0));
611    }
612
613    #[test]
614    fn test_continuous() {
615        test::check_continuous_distribution(&try_create(1.0), 0.0, 10.0);
616        test::check_continuous_distribution(&try_create(2.0), 0.0, 10.0);
617        test::check_continuous_distribution(&try_create(5.0), 0.0, 10.0);
618    }
619}
statrs/distribution/chi.rs

statrs/distribution/
chi.rs
