rand_distr/
chi_squared.rs

1// Copyright 2018 Developers of the Rand project.
2// Copyright 2013 The Rust Project Developers.
3//
4// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
5// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
6// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
7// option. This file may not be copied, modified, or distributed
8// except according to those terms.
9
10//! The Chi-squared distribution.
11
12use self::ChiSquaredRepr::*;
13
14use crate::{Distribution, Exp1, Gamma, Open01, StandardNormal};
15use core::fmt;
16use num_traits::Float;
17use rand::Rng;
18#[cfg(feature = "serde")]
19use serde::{Deserialize, Serialize};
20
21/// The [chi-squared distribution](https://en.wikipedia.org/wiki/Chi-squared_distribution) `χ²(k)`.
22///
23/// The chi-squared distribution is a continuous probability
24/// distribution with parameter `k > 0` degrees of freedom.
25///
26/// For `k > 0` integral, this distribution is the sum of the squares
27/// of `k` independent standard normal random variables. For other
28/// `k`, this uses the equivalent characterisation
29/// `χ²(k) = Gamma(k/2, 2)`.
30///
31/// # Plot
32///
33/// The plot shows the chi-squared distribution with various degrees
34/// of freedom.
35///
36/// ![Chi-squared distribution](https://raw.githubusercontent.com/rust-random/charts/main/charts/chi_squared.svg)
37///
38/// # Example
39///
40/// ```
41/// use rand_distr::{ChiSquared, Distribution};
42///
43/// let chi = ChiSquared::new(11.0).unwrap();
44/// let v = chi.sample(&mut rand::rng());
45/// println!("{} is from a χ²(11) distribution", v)
46/// ```
47#[derive(Clone, Copy, Debug, PartialEq)]
48#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
49pub struct ChiSquared<F>
50where
51    F: Float,
52    StandardNormal: Distribution<F>,
53    Exp1: Distribution<F>,
54    Open01: Distribution<F>,
55{
56    repr: ChiSquaredRepr<F>,
57}
58
59/// Error type returned from [`ChiSquared::new`] and [`StudentT::new`](crate::StudentT::new).
60#[derive(Clone, Copy, Debug, PartialEq, Eq)]
61#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
62pub enum Error {
63    /// `0.5 * k <= 0` or `nan`.
64    DoFTooSmall,
65}
66
67impl fmt::Display for Error {
68    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
69        f.write_str(match self {
70            Error::DoFTooSmall => {
71                "degrees-of-freedom k is not positive in chi-squared distribution"
72            }
73        })
74    }
75}
76
77#[cfg(feature = "std")]
78impl std::error::Error for Error {}
79
80#[derive(Clone, Copy, Debug, PartialEq)]
81#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
82enum ChiSquaredRepr<F>
83where
84    F: Float,
85    StandardNormal: Distribution<F>,
86    Exp1: Distribution<F>,
87    Open01: Distribution<F>,
88{
89    // k == 1, Gamma(alpha, ..) is particularly slow for alpha < 1,
90    // e.g. when alpha = 1/2 as it would be for this case, so special-
91    // casing and using the definition of N(0,1)^2 is faster.
92    DoFExactlyOne,
93    DoFAnythingElse(Gamma<F>),
94}
95
96impl<F> ChiSquared<F>
97where
98    F: Float,
99    StandardNormal: Distribution<F>,
100    Exp1: Distribution<F>,
101    Open01: Distribution<F>,
102{
103    /// Create a new chi-squared distribution with degrees-of-freedom
104    /// `k`.
105    pub fn new(k: F) -> Result<ChiSquared<F>, Error> {
106        let repr = if k == F::one() {
107            DoFExactlyOne
108        } else {
109            if !(F::from(0.5).unwrap() * k > F::zero()) {
110                return Err(Error::DoFTooSmall);
111            }
112            DoFAnythingElse(Gamma::new(F::from(0.5).unwrap() * k, F::from(2.0).unwrap()).unwrap())
113        };
114        Ok(ChiSquared { repr })
115    }
116}
117impl<F> Distribution<F> for ChiSquared<F>
118where
119    F: Float,
120    StandardNormal: Distribution<F>,
121    Exp1: Distribution<F>,
122    Open01: Distribution<F>,
123{
124    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> F {
125        match self.repr {
126            DoFExactlyOne => {
127                // k == 1 => N(0,1)^2
128                let norm: F = rng.sample(StandardNormal);
129                norm * norm
130            }
131            DoFAnythingElse(ref g) => g.sample(rng),
132        }
133    }
134}
135
136#[cfg(test)]
137mod test {
138    use super::*;
139
140    #[test]
141    fn test_chi_squared_one() {
142        let chi = ChiSquared::new(1.0).unwrap();
143        let mut rng = crate::test::rng(201);
144        for _ in 0..1000 {
145            chi.sample(&mut rng);
146        }
147    }
148    #[test]
149    fn test_chi_squared_small() {
150        let chi = ChiSquared::new(0.5).unwrap();
151        let mut rng = crate::test::rng(202);
152        for _ in 0..1000 {
153            chi.sample(&mut rng);
154        }
155    }
156    #[test]
157    fn test_chi_squared_large() {
158        let chi = ChiSquared::new(30.0).unwrap();
159        let mut rng = crate::test::rng(203);
160        for _ in 0..1000 {
161            chi.sample(&mut rng);
162        }
163    }
164    #[test]
165    #[should_panic]
166    fn test_chi_squared_invalid_dof() {
167        ChiSquared::new(-1.0).unwrap();
168    }
169
170    #[test]
171    fn gamma_distributions_can_be_compared() {
172        assert_eq!(Gamma::new(1.0, 2.0), Gamma::new(1.0, 2.0));
173    }
174
175    #[test]
176    fn chi_squared_distributions_can_be_compared() {
177        assert_eq!(ChiSquared::new(1.0), ChiSquared::new(1.0));
178    }
179}
rand_distr/chi_squared.rs

rand_distr/
chi_squared.rs
