statrs/distribution/
chi_squared.rs

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

statrs/distribution/
chi_squared.rs
