statrs/distribution/
chi.rs

1use crate::distribution::{Continuous, ContinuousCDF};
2use crate::function::gamma;
3use crate::statistics::*;
4use std::f64;
5use std::num::NonZeroU64;
6
7/// Implements the [Chi](https://en.wikipedia.org/wiki/Chi_distribution)
8/// distribution
9///
10/// # Examples
11///
12/// ```
13/// use statrs::distribution::{Chi, Continuous};
14/// use statrs::statistics::Distribution;
15/// use statrs::prec;
16///
17/// let n = Chi::new(2).unwrap();
18/// assert!(prec::almost_eq(n.mean().unwrap(), 1.25331413731550025121, 1e-14));
19/// assert!(prec::almost_eq(n.pdf(1.0), 0.60653065971263342360, 1e-15));
20/// ```
21#[derive(Copy, Clone, PartialEq, Debug)]
22pub struct Chi {
23    freedom: NonZeroU64,
24}
25
26/// Represents the errors that can occur when creating a [`Chi`].
27#[derive(Copy, Clone, PartialEq, Eq, Debug, Hash)]
28#[non_exhaustive]
29pub enum ChiError {
30    /// The degrees of freedom are zero.
31    FreedomInvalid,
32}
33
34impl std::fmt::Display for ChiError {
35    #[cfg_attr(coverage_nightly, coverage(off))]
36    fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
37        match self {
38            ChiError::FreedomInvalid => {
39                write!(f, "Degrees of freedom are zero")
40            }
41        }
42    }
43}
44
45impl std::error::Error for ChiError {}
46
47impl Chi {
48    /// Constructs a new chi distribution
49    /// with `freedom` degrees of freedom
50    ///
51    /// # Errors
52    ///
53    /// Returns an error if `freedom` is equal to `0`.
54    ///
55    /// # Examples
56    ///
57    /// ```
58    /// use statrs::distribution::Chi;
59    ///
60    /// let mut result = Chi::new(2);
61    /// assert!(result.is_ok());
62    ///
63    /// result = Chi::new(0);
64    /// assert!(result.is_err());
65    /// ```
66    pub fn new(freedom: u64) -> Result<Chi, ChiError> {
67        match NonZeroU64::new(freedom) {
68            Some(freedom) => Ok(Self { freedom }),
69            None => Err(ChiError::FreedomInvalid),
70        }
71    }
72
73    /// Returns the degrees of freedom of the chi distribution.
74    /// Guaranteed to be non-zero.
75    ///
76    /// # Examples
77    ///
78    /// ```
79    /// use statrs::distribution::Chi;
80    ///
81    /// let n = Chi::new(2).unwrap();
82    /// assert_eq!(n.freedom(), 2);
83    /// ```
84    pub fn freedom(&self) -> u64 {
85        self.freedom.get()
86    }
87}
88
89impl std::fmt::Display for Chi {
90    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
91        write!(f, "χ_{}", self.freedom)
92    }
93}
94
95#[cfg(feature = "rand")]
96#[cfg_attr(docsrs, doc(cfg(feature = "rand")))]
97impl ::rand::distributions::Distribution<f64> for Chi {
98    fn sample<R: ::rand::Rng + ?Sized>(&self, rng: &mut R) -> f64 {
99        (0..self.freedom())
100            .fold(0.0, |acc, _| {
101                acc + super::normal::sample_unchecked(rng, 0.0, 1.0).powf(2.0)
102            })
103            .sqrt()
104    }
105}
106
107impl ContinuousCDF<f64, f64> for Chi {
108    /// Calculates the cumulative distribution function for the chi
109    /// distribution at `x`.
110    ///
111    /// # Formula
112    ///
113    /// ```text
114    /// P(k / 2, x^2 / 2)
115    /// ```
116    ///
117    /// where `k` is the degrees of freedom and `P` is
118    /// the regularized lower incomplete Gamma function
119    fn cdf(&self, x: f64) -> f64 {
120        if x == f64::INFINITY {
121            1.0
122        } else if x <= 0.0 {
123            0.0
124        } else {
125            gamma::gamma_lr(self.freedom() as f64 / 2.0, x * x / 2.0)
126        }
127    }
128
129    /// Calculates the survival function for the chi
130    /// distribution at `x`.
131    ///
132    /// # Formula
133    ///
134    /// ```text
135    /// P(k / 2, x^2 / 2)
136    /// ```
137    ///
138    /// where `k` is the degrees of freedom and `P` is
139    /// the regularized upper incomplete Gamma function
140    fn sf(&self, x: f64) -> f64 {
141        if x == f64::INFINITY {
142            0.0
143        } else if x <= 0.0 {
144            1.0
145        } else {
146            gamma::gamma_ur(self.freedom() as f64 / 2.0, x * x / 2.0)
147        }
148    }
149}
150
151impl Min<f64> for Chi {
152    /// Returns the minimum value in the domain of the chi distribution
153    /// representable by a double precision float
154    ///
155    /// # Formula
156    ///
157    /// ```text
158    /// 0
159    /// ```
160    fn min(&self) -> f64 {
161        0.0
162    }
163}
164
165impl Max<f64> for Chi {
166    /// Returns the maximum value in the domain of the chi distribution
167    /// representable by a double precision float
168    ///
169    /// # Formula
170    ///
171    /// ```text
172    /// f64::INFINITY
173    /// ```
174    fn max(&self) -> f64 {
175        f64::INFINITY
176    }
177}
178
179impl Distribution<f64> for Chi {
180    /// Returns the mean of the chi distribution
181    ///
182    /// # Remarks
183    ///
184    /// Returns `NaN` if `freedom` is `INF`
185    ///
186    /// # Formula
187    ///
188    /// ```text
189    /// sqrt2 * Γ((k + 1) / 2) / Γ(k / 2)
190    /// ```
191    ///
192    /// where `k` is degrees of freedom and `Γ` is the gamma function
193    fn mean(&self) -> Option<f64> {
194        let freedom = self.freedom() as f64;
195
196        if self.freedom() > 300 {
197            // Large n approximation based on the Stirling series approximation to the Gamma function
198            // This avoids call the Gamma function with large arguments and returning NaN
199            //
200            // Relative accuracy follows O(1/n^4) and at 300 d.o.f. is better than 1e-12
201            // For a f32 impl the threshold should be changed to 150
202            Some(
203                (freedom.sqrt())
204                    / ((1.0 + 0.25 / freedom)
205                        * (1.0 + 0.03125 / (freedom * freedom))
206                        * (1.0 - 0.046875 / (freedom * freedom * freedom))),
207            )
208        } else {
209            let mean = f64::consts::SQRT_2 * gamma::gamma((freedom + 1.0) / 2.0)
210                / gamma::gamma(freedom / 2.0);
211            Some(mean)
212        }
213    }
214
215    /// Returns the variance of the chi distribution
216    ///
217    /// # Remarks
218    ///
219    /// Returns `NaN` if `freedom` is `INF`
220    ///
221    /// # Formula
222    ///
223    /// ```text
224    /// k - μ^2
225    /// ```
226    ///
227    /// where `k` is degrees of freedom and `μ` is the mean
228    /// of the distribution
229    fn variance(&self) -> Option<f64> {
230        let mean = self.mean()?;
231        Some(self.freedom() as f64 - mean * mean)
232    }
233
234    /// Returns the entropy of the chi distribution
235    ///
236    /// # Remarks
237    ///
238    /// Returns `None` if `freedom` is `INF`
239    ///
240    /// # Formula
241    ///
242    /// ```text
243    /// ln(Γ(k / 2)) + 0.5 * (k - ln2 - (k - 1) * ψ(k / 2))
244    /// ```
245    ///
246    /// where `k` is degrees of freedom, `Γ` is the gamma function,
247    /// and `ψ` is the digamma function
248    fn entropy(&self) -> Option<f64> {
249        let freedom = self.freedom() as f64;
250        let entr = gamma::ln_gamma(freedom / 2.0)
251            + (freedom - (2.0f64).ln() - (freedom - 1.0) * gamma::digamma(freedom / 2.0)) / 2.0;
252        Some(entr)
253    }
254
255    /// Returns the skewness of the chi distribution
256    ///
257    /// # Remarks
258    ///
259    /// Returns `NaN` if `freedom` is `INF`
260    ///
261    /// # Formula
262    ///
263    /// ```text
264    /// (μ / σ^3) * (1 - 2σ^2)
265    /// ```
266    /// where `μ` is the mean and `σ` the standard deviation
267    /// of the distribution
268    fn skewness(&self) -> Option<f64> {
269        let sigma = self.std_dev()?;
270        let skew = self.mean()? * (1.0 - 2.0 * sigma * sigma) / (sigma * sigma * sigma);
271        Some(skew)
272    }
273}
274
275impl Mode<Option<f64>> for Chi {
276    /// Returns the mode for the chi distribution
277    ///
278    /// # Panics
279    ///
280    /// If `freedom < 1.0`
281    ///
282    /// # Formula
283    ///
284    /// ```text
285    /// sqrt(k - 1)
286    /// ```
287    ///
288    /// where `k` is the degrees of freedom
289    fn mode(&self) -> Option<f64> {
290        Some(((self.freedom() - 1) as f64).sqrt())
291    }
292}
293
294impl Continuous<f64, f64> for Chi {
295    /// Calculates the probability density function for the chi
296    /// distribution at `x`
297    ///
298    /// # Formula
299    ///
300    /// ```text
301    /// (2^(1 - (k / 2)) * x^(k - 1) * e^(-x^2 / 2)) / Γ(k / 2)
302    /// ```
303    ///
304    /// where `k` is the degrees of freedom and `Γ` is the gamma function
305    fn pdf(&self, x: f64) -> f64 {
306        if x == f64::INFINITY || x <= 0.0 {
307            0.0
308        } else if self.freedom() > 160 {
309            self.ln_pdf(x).exp()
310        } else {
311            let freedom = self.freedom() as f64;
312            (2.0f64).powf(1.0 - freedom / 2.0) * x.powf(freedom - 1.0) * (-x * x / 2.0).exp()
313                / gamma::gamma(freedom / 2.0)
314        }
315    }
316
317    /// Calculates the log probability density function for the chi distribution
318    /// at `x`
319    ///
320    /// # Formula
321    ///
322    /// ```text
323    /// ln((2^(1 - (k / 2)) * x^(k - 1) * e^(-x^2 / 2)) / Γ(k / 2))
324    /// ```
325    fn ln_pdf(&self, x: f64) -> f64 {
326        if x == f64::INFINITY || x <= 0.0 {
327            f64::NEG_INFINITY
328        } else {
329            let freedom = self.freedom() as f64;
330            (1.0 - freedom / 2.0) * (2.0f64).ln() + ((freedom - 1.0) * x.ln())
331                - x * x / 2.0
332                - gamma::ln_gamma(freedom / 2.0)
333        }
334    }
335}
336
337#[rustfmt::skip]
338#[cfg(test)]
339mod tests {
340    use super::*;
341    use crate::distribution::internal::*;
342    use crate::testing_boiler;
343
344    testing_boiler!(freedom: u64; Chi; ChiError);
345
346    #[test]
347    fn test_create() {
348        create_ok(1);
349        create_ok(3);
350    }
351
352    #[test]
353    fn test_bad_create() {
354        create_err(0);
355    }
356
357    #[test]
358    fn test_mean() {
359        let mean = |x: Chi| x.mean().unwrap();
360        test_absolute(1, 0.7978845608028653558799, 1e-15, mean);
361        test_absolute(2, 1.25331413731550025121, 1e-14, mean);
362        test_absolute(5, 2.12769216214097428235, 1e-14, mean);
363        test_absolute(336, 18.31666925443713, 1e-12, mean);
364    }
365
366    #[test]
367    fn test_large_dof_mean_not_nan() {
368        for i in 1..1000 {
369            let mean = Chi::new(i).unwrap().mean().unwrap();
370            assert!(!mean.is_nan(), "Chi mean for {i} dof was {mean}");
371        }
372    }
373
374    #[test]
375    fn test_variance() {
376        let variance = |x: Chi| x.variance().unwrap();
377        test_absolute(1, 0.3633802276324186569245, 1e-15, variance);
378        test_absolute(2, 0.42920367320510338077, 1e-14, variance);
379        test_absolute(3, 0.4535209105296746277, 1e-14, variance);
380    }
381
382    #[test]
383    fn test_entropy() {
384        let entropy = |x: Chi| x.entropy().unwrap();
385        test_absolute(1, 0.7257913526447274323631, 1e-15, entropy);
386        test_absolute(2, 0.9420342421707937755946, 1e-15, entropy);
387        test_absolute(3, 0.99615419810620560239, 1e-14, entropy);
388    }
389
390    #[test]
391    fn test_skewness() {
392        let skewness = |x: Chi| x.skewness().unwrap();
393        test_absolute(1, 0.995271746431156042444, 1e-14, skewness);
394        test_absolute(3, 0.485692828049590809, 1e-12, skewness);
395    }
396
397    #[test]
398    fn test_mode() {
399        let mode = |x: Chi| x.mode().unwrap();
400        test_exact(1, 0.0, mode);
401        test_exact(2, 1.0, mode);
402        test_exact(3, f64::consts::SQRT_2, mode);
403    }
404
405    #[test]
406    fn test_min_max() {
407        let min = |x: Chi| x.min();
408        let max = |x: Chi| x.max();
409        test_exact(1, 0.0, min);
410        test_exact(2, 0.0, min);
411        test_exact(2, 0.0, min);
412        test_exact(3, 0.0, min);
413        test_exact(1, f64::INFINITY, max);
414        test_exact(2, f64::INFINITY, max);
415        test_exact(2, f64::INFINITY, max);
416        test_exact(3, f64::INFINITY, max);
417    }
418
419    #[test]
420    fn test_pdf() {
421        let pdf = |arg: f64| move |x: Chi| x.pdf(arg);
422        test_exact(1, 0.0, pdf(0.0));
423        test_absolute(1, 0.79390509495402353102, 1e-15, pdf(0.1));
424        test_absolute(1, 0.48394144903828669960, 1e-15, pdf(1.0));
425        test_absolute(1, 2.1539520085086552718e-7, 1e-22, pdf(5.5));
426        test_exact(1, 0.0, pdf(f64::INFINITY));
427        test_exact(2, 0.0, pdf(0.0));
428        test_absolute(2, 0.099501247919268231335, 1e-16, pdf(0.1));
429        test_absolute(2, 0.60653065971263342360, 1e-15, pdf(1.0));
430        test_absolute(2, 1.4847681768496578863e-6, 1e-21, pdf(5.5));
431        test_exact(2, 0.0, pdf(f64::INFINITY));
432        test_exact(2, 0.0, pdf(0.0));
433        test_exact(2, 0.0, pdf(f64::INFINITY));
434        test_absolute(170, 0.5644678498668440878, 1e-13, pdf(13.0));
435    }
436
437    #[test]
438    fn test_neg_pdf() {
439        let pdf = |arg: f64| move |x: Chi| x.pdf(arg);
440        test_exact(1, 0.0, pdf(-1.0));
441    }
442
443    #[test]
444    fn test_ln_pdf() {
445        let ln_pdf = |arg: f64| move |x: Chi| x.ln_pdf(arg);
446        test_exact(1, f64::NEG_INFINITY, ln_pdf(0.0));
447        test_absolute(1, -0.23079135264472743236, 1e-15, ln_pdf(0.1));
448        test_absolute(1, -0.72579135264472743236, 1e-15, ln_pdf(1.0));
449        test_absolute(1, -15.350791352644727432, 1e-14, ln_pdf(5.5));
450        test_exact(1, f64::NEG_INFINITY, ln_pdf(f64::INFINITY));
451        test_exact(2, f64::NEG_INFINITY, ln_pdf(0.0));
452        test_absolute(2, -2.3075850929940456840, 1e-15, ln_pdf(0.1));
453        test_absolute(2, -0.5, 1e-15, ln_pdf(1.0));
454        test_absolute(2, -13.420251907761574765, 1e-15, ln_pdf(5.5));
455        test_exact(2, f64::NEG_INFINITY, ln_pdf(f64::INFINITY));
456        test_exact(2, f64::NEG_INFINITY, ln_pdf(0.0));
457        test_exact(2, f64::NEG_INFINITY, ln_pdf(f64::INFINITY));
458        test_absolute(170, -0.57187185030600516424237, 1e-13, ln_pdf(13.0));
459    }
460
461    #[test]
462    fn test_neg_ln_pdf() {
463        let ln_pdf = |arg: f64| move |x: Chi| x.ln_pdf(arg);
464        test_exact(1, f64::NEG_INFINITY, ln_pdf(-1.0));
465    }
466
467    #[test]
468    fn test_cdf() {
469        let cdf = |arg: f64| move |x: Chi| x.cdf(arg);
470        test_exact(1, 0.0, cdf(0.0));
471        test_absolute(1, 0.079655674554057962931, 1e-16, cdf(0.1));
472        test_absolute(1, 0.68268949213708589717, 1e-15, cdf(1.0));
473        test_exact(1, 0.99999996202087506822, cdf(5.5));
474        test_exact(1, 1.0, cdf(f64::INFINITY));
475        test_exact(2, 0.0, cdf(0.0));
476        test_absolute(2, 0.0049875208073176866474, 1e-17, cdf(0.1));
477        test_exact(2, 1.0, cdf(f64::INFINITY));
478        test_exact(2, 0.0, cdf(0.0));
479        test_exact(2, 1.0, cdf(f64::INFINITY));
480    }
481
482    #[test]
483    fn test_sf() {
484        let sf = |arg: f64| move |x: Chi| x.sf(arg);
485        test_exact(1, 1.0, sf(0.0));
486        test_absolute(1, 0.920344325445942, 1e-16, sf(0.1));
487        test_absolute(1, 0.31731050786291404, 1e-15, sf(1.0));
488        test_absolute(1, 3.797912493177544e-8, 1e-15, sf(5.5));
489        test_exact(1, 0.0, sf(f64::INFINITY));
490        test_exact(2, 1.0, sf(0.0));
491        test_absolute(2, 0.9950124791926823, 1e-17, sf(0.1));
492        test_absolute(2, 0.6065306597126333, 1e-15, sf(1.0));
493        test_absolute(2, 2.699578503363014e-7, 1e-15, sf(5.5));
494        test_exact(2, 0.0, sf(f64::INFINITY));
495        test_exact(2, 1.0, sf(0.0));
496        test_exact(2, 0.0, sf(f64::INFINITY));
497    }
498
499    #[test]
500    fn test_neg_cdf() {
501        let cdf = |arg: f64| move |x: Chi| x.cdf(arg);
502        test_exact(1, 0.0, cdf(-1.0));
503    }
504
505    #[test]
506    fn test_neg_sf() {
507        let sf = |arg: f64| move |x: Chi| x.sf(arg);
508        test_exact(1, 1.0, sf(-1.0));
509    }
510
511    #[test]
512    fn test_continuous() {
513        test::check_continuous_distribution(&create_ok(1), 0.0, 10.0);
514        test::check_continuous_distribution(&create_ok(2), 0.0, 10.0);
515        test::check_continuous_distribution(&create_ok(5), 0.0, 10.0);
516    }
517}
statrs/distribution/chi.rs

statrs/distribution/
chi.rs
