statrs/distribution/
chi_squared.rs

1use crate::distribution::{Continuous, ContinuousCDF, Gamma, GammaError};
2use crate::statistics::*;
3use std::f64;
4
5/// Implements the
6/// [Chi-squared](https://en.wikipedia.org/wiki/Chi-squared_distribution)
7/// distribution which is a special case of the
8/// [Gamma](https://en.wikipedia.org/wiki/Gamma_distribution) distribution
9/// (referenced [Here](./struct.Gamma.html))
10///
11/// # Examples
12///
13/// ```
14/// use statrs::distribution::{ChiSquared, Continuous};
15/// use statrs::statistics::Distribution;
16/// use statrs::prec;
17///
18/// let n = ChiSquared::new(3.0).unwrap();
19/// assert_eq!(n.mean().unwrap(), 3.0);
20/// assert!(prec::almost_eq(n.pdf(4.0), 0.107981933026376103901, 1e-15));
21/// ```
22#[derive(Copy, Clone, PartialEq, Debug)]
23pub struct ChiSquared {
24    freedom: f64,
25    g: Gamma,
26}
27
28impl ChiSquared {
29    /// Constructs a new chi-squared distribution with `freedom`
30    /// degrees of freedom. This is equivalent to a Gamma distribution
31    /// with a shape of `freedom / 2.0` and a rate of `0.5`.
32    ///
33    /// # Errors
34    ///
35    /// Returns an error if `freedom` is `NaN` or less than
36    /// or equal to `0.0`
37    ///
38    /// # Examples
39    ///
40    /// ```
41    /// use statrs::distribution::ChiSquared;
42    ///
43    /// let mut result = ChiSquared::new(3.0);
44    /// assert!(result.is_ok());
45    ///
46    /// result = ChiSquared::new(0.0);
47    /// assert!(result.is_err());
48    /// ```
49    pub fn new(freedom: f64) -> Result<ChiSquared, GammaError> {
50        Gamma::new(freedom / 2.0, 0.5).map(|g| ChiSquared { freedom, g })
51    }
52
53    /// Returns the degrees of freedom of the chi-squared
54    /// distribution
55    ///
56    /// # Examples
57    ///
58    /// ```
59    /// use statrs::distribution::ChiSquared;
60    ///
61    /// let n = ChiSquared::new(3.0).unwrap();
62    /// assert_eq!(n.freedom(), 3.0);
63    /// ```
64    pub fn freedom(&self) -> f64 {
65        self.freedom
66    }
67
68    /// Returns the shape of the underlying Gamma distribution
69    ///
70    /// # Examples
71    ///
72    /// ```
73    /// use statrs::distribution::ChiSquared;
74    ///
75    /// let n = ChiSquared::new(3.0).unwrap();
76    /// assert_eq!(n.shape(), 3.0 / 2.0);
77    /// ```
78    pub fn shape(&self) -> f64 {
79        self.g.shape()
80    }
81
82    /// Returns the rate of the underlying Gamma distribution
83    ///
84    /// # Examples
85    ///
86    /// ```
87    /// use statrs::distribution::ChiSquared;
88    ///
89    /// let n = ChiSquared::new(3.0).unwrap();
90    /// assert_eq!(n.rate(), 0.5);
91    /// ```
92    pub fn rate(&self) -> f64 {
93        self.g.rate()
94    }
95}
96
97impl std::fmt::Display for ChiSquared {
98    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
99        write!(f, "χ^2_{}", self.freedom)
100    }
101}
102
103#[cfg(feature = "rand")]
104#[cfg_attr(docsrs, doc(cfg(feature = "rand")))]
105impl ::rand::distributions::Distribution<f64> for ChiSquared {
106    fn sample<R: ::rand::Rng + ?Sized>(&self, r: &mut R) -> f64 {
107        ::rand::distributions::Distribution::sample(&self.g, r)
108    }
109}
110
111impl ContinuousCDF<f64, f64> for ChiSquared {
112    /// Calculates the cumulative distribution function for the
113    /// chi-squared distribution at `x`
114    ///
115    /// # Formula
116    ///
117    /// ```text
118    /// (1 / Γ(k / 2)) * γ(k / 2, x / 2)
119    /// ```
120    ///
121    /// where `k` is the degrees of freedom, `Γ` is the gamma function,
122    /// and `γ` is the lower incomplete gamma function
123    fn cdf(&self, x: f64) -> f64 {
124        self.g.cdf(x)
125    }
126
127    /// Calculates the cumulative distribution function for the
128    /// chi-squared distribution at `x`
129    ///
130    /// # Formula
131    ///
132    /// ```text
133    /// (1 / Γ(k / 2)) * γ(k / 2, x / 2)
134    /// ```
135    ///
136    /// where `k` is the degrees of freedom, `Γ` is the gamma function,
137    /// and `γ` is the upper incomplete gamma function
138    fn sf(&self, x: f64) -> f64 {
139        self.g.sf(x)
140    }
141
142    /// Calculates the inverse cumulative distribution function for the
143    /// chi-squared distribution at `x`
144    ///
145    /// # Formula
146    ///
147    /// ```text
148    /// γ^{-1}(k / 2, x * Γ(k / 2) / 2)
149    /// ```
150    ///
151    /// where `k` is the degrees of freedom, `Γ` is the gamma function,
152    /// and `γ` is the lower incomplete gamma function
153    fn inverse_cdf(&self, p: f64) -> f64 {
154        self.g.inverse_cdf(p)
155    }
156}
157
158impl Min<f64> for ChiSquared {
159    /// Returns the minimum value in the domain of the
160    /// chi-squared distribution representable by a double precision
161    /// float
162    ///
163    /// # Formula
164    ///
165    /// ```text
166    /// 0
167    /// ```
168    fn min(&self) -> f64 {
169        0.0
170    }
171}
172
173impl Max<f64> for ChiSquared {
174    /// Returns the maximum value in the domain of the
175    /// chi-squared distribution representable by a double precision
176    /// float
177    ///
178    /// # Formula
179    ///
180    /// ```text
181    /// f64::INFINITY
182    /// ```
183    fn max(&self) -> f64 {
184        f64::INFINITY
185    }
186}
187
188impl Distribution<f64> for ChiSquared {
189    /// Returns the mean of the chi-squared distribution
190    ///
191    /// # Formula
192    ///
193    /// ```text
194    /// k
195    /// ```
196    ///
197    /// where `k` is the degrees of freedom
198    fn mean(&self) -> Option<f64> {
199        self.g.mean()
200    }
201
202    /// Returns the variance of the chi-squared distribution
203    ///
204    /// # Formula
205    ///
206    /// ```text
207    /// 2k
208    /// ```
209    ///
210    /// where `k` is the degrees of freedom
211    fn variance(&self) -> Option<f64> {
212        self.g.variance()
213    }
214
215    /// Returns the entropy of the chi-squared distribution
216    ///
217    /// # Formula
218    ///
219    /// ```text
220    /// (k / 2) + ln(2 * Γ(k / 2)) + (1 - (k / 2)) * ψ(k / 2)
221    /// ```
222    ///
223    /// where `k` is the degrees of freedom, `Γ` is the gamma function,
224    /// and `ψ` is the digamma function
225    fn entropy(&self) -> Option<f64> {
226        self.g.entropy()
227    }
228
229    /// Returns the skewness of the chi-squared distribution
230    ///
231    /// # Formula
232    ///
233    /// ```text
234    /// sqrt(8 / k)
235    /// ```
236    ///
237    /// where `k` is the degrees of freedom
238    fn skewness(&self) -> Option<f64> {
239        self.g.skewness()
240    }
241}
242
243impl Median<f64> for ChiSquared {
244    /// Returns the median  of the chi-squared distribution
245    ///
246    /// # Formula
247    ///
248    /// ```text
249    /// k * (1 - (2 / 9k))^3
250    /// ```
251    fn median(&self) -> f64 {
252        if self.freedom < 1.0 {
253            // if k is small, calculate using expansion of formula
254            self.freedom - 2.0 / 3.0 + 12.0 / (81.0 * self.freedom)
255                - 8.0 / (729.0 * self.freedom * self.freedom)
256        } else {
257            // if k is large enough, median heads toward k - 2/3
258            self.freedom - 2.0 / 3.0
259        }
260    }
261}
262
263impl Mode<Option<f64>> for ChiSquared {
264    /// Returns the mode of the chi-squared distribution
265    ///
266    /// # Formula
267    ///
268    /// ```text
269    /// k - 2
270    /// ```
271    ///
272    /// where `k` is the degrees of freedom
273    fn mode(&self) -> Option<f64> {
274        self.g.mode()
275    }
276}
277
278impl Continuous<f64, f64> for ChiSquared {
279    /// Calculates the probability density function for the chi-squared
280    /// distribution at `x`
281    ///
282    /// # Formula
283    ///
284    /// ```text
285    /// 1 / (2^(k / 2) * Γ(k / 2)) * x^((k / 2) - 1) * e^(-x / 2)
286    /// ```
287    ///
288    /// where `k` is the degrees of freedom and `Γ` is the gamma function
289    fn pdf(&self, x: f64) -> f64 {
290        self.g.pdf(x)
291    }
292
293    /// Calculates the log probability density function for the chi-squared
294    /// distribution at `x`
295    ///
296    /// # Formula
297    ///
298    /// ```text
299    /// ln(1 / (2^(k / 2) * Γ(k / 2)) * x^((k / 2) - 1) * e^(-x / 2))
300    /// ```
301    fn ln_pdf(&self, x: f64) -> f64 {
302        self.g.ln_pdf(x)
303    }
304}
305
306#[rustfmt::skip]
307#[cfg(test)]
308mod tests {
309    use super::*;
310    use crate::distribution::internal::*;
311    use crate::testing_boiler;
312
313    testing_boiler!(freedom: f64; ChiSquared; GammaError);
314
315    #[test]
316    fn test_median() {
317        let median = |x: ChiSquared| x.median();
318        test_absolute(0.5, 0.0857338820301783264746, 1e-16, median);
319        test_exact(1.0, 1.0 - 2.0 / 3.0, median);
320        test_exact(2.0, 2.0 - 2.0 / 3.0, median);
321        test_exact(2.5, 2.5 - 2.0 / 3.0, median);
322        test_exact(3.0, 3.0 - 2.0 / 3.0, median);
323    }
324
325    #[test]
326    fn test_continuous() {
327        // TODO: figure out why this test fails:
328        //test::check_continuous_distribution(&create_ok(1.0), 0.0, 10.0);
329        test::check_continuous_distribution(&create_ok(2.0), 0.0, 10.0);
330        test::check_continuous_distribution(&create_ok(5.0), 0.0, 50.0);
331    }
332}
statrs/distribution/chi_squared.rs

statrs/distribution/
chi_squared.rs
