statrs/distribution/
students_t.rs

1use crate::distribution::{Continuous, ContinuousCDF};
2use crate::function::{beta, gamma};
3use crate::is_zero;
4use crate::statistics::*;
5use crate::{Result, StatsError};
6use rand::Rng;
7use std::f64;
8
9/// Implements the [Student's
10/// T](https://en.wikipedia.org/wiki/Student%27s_t-distribution) distribution
11///
12/// # Examples
13///
14/// ```
15/// use statrs::distribution::{StudentsT, Continuous};
16/// use statrs::statistics::Distribution;
17/// use statrs::prec;
18///
19/// let n = StudentsT::new(0.0, 1.0, 2.0).unwrap();
20/// assert_eq!(n.mean().unwrap(), 0.0);
21/// assert!(prec::almost_eq(n.pdf(0.0), 0.353553390593274, 1e-15));
22/// ```
23#[derive(Debug, Copy, Clone, PartialEq)]
24pub struct StudentsT {
25    location: f64,
26    scale: f64,
27    freedom: f64,
28}
29
30impl StudentsT {
31    /// Constructs a new student's t-distribution with location `location`,
32    /// scale `scale`,
33    /// and `freedom` freedom.
34    ///
35    /// # Errors
36    ///
37    /// Returns an error if any of `location`, `scale`, or `freedom` are `NaN`.
38    /// Returns an error if `scale <= 0.0` or `freedom <= 0.0`
39    ///
40    /// # Examples
41    ///
42    /// ```
43    /// use statrs::distribution::StudentsT;
44    ///
45    /// let mut result = StudentsT::new(0.0, 1.0, 2.0);
46    /// assert!(result.is_ok());
47    ///
48    /// result = StudentsT::new(0.0, 0.0, 0.0);
49    /// assert!(result.is_err());
50    /// ```
51    pub fn new(location: f64, scale: f64, freedom: f64) -> Result<StudentsT> {
52        let is_nan = location.is_nan() || scale.is_nan() || freedom.is_nan();
53        if is_nan || scale <= 0.0 || freedom <= 0.0 {
54            Err(StatsError::BadParams)
55        } else {
56            Ok(StudentsT {
57                location,
58                scale,
59                freedom,
60            })
61        }
62    }
63
64    /// Returns the location of the student's t-distribution
65    ///
66    /// # Examples
67    ///
68    /// ```
69    /// use statrs::distribution::StudentsT;
70    ///
71    /// let n = StudentsT::new(0.0, 1.0, 2.0).unwrap();
72    /// assert_eq!(n.location(), 0.0);
73    /// ```
74    pub fn location(&self) -> f64 {
75        self.location
76    }
77
78    /// Returns the scale of the student's t-distribution
79    ///
80    /// # Examples
81    ///
82    /// ```
83    /// use statrs::distribution::StudentsT;
84    ///
85    /// let n = StudentsT::new(0.0, 1.0, 2.0).unwrap();
86    /// assert_eq!(n.scale(), 1.0);
87    /// ```
88    pub fn scale(&self) -> f64 {
89        self.scale
90    }
91
92    /// Returns the freedom of the student's t-distribution
93    ///
94    /// # Examples
95    ///
96    /// ```
97    /// use statrs::distribution::StudentsT;
98    ///
99    /// let n = StudentsT::new(0.0, 1.0, 2.0).unwrap();
100    /// assert_eq!(n.freedom(), 2.0);
101    /// ```
102    pub fn freedom(&self) -> f64 {
103        self.freedom
104    }
105}
106
107impl ::rand::distributions::Distribution<f64> for StudentsT {
108    fn sample<R: Rng + ?Sized>(&self, r: &mut R) -> f64 {
109        // based on method 2, section 5 in chapter 9 of L. Devroye's
110        // "Non-Uniform Random Variate Generation"
111        let gamma = super::gamma::sample_unchecked(r, 0.5 * self.freedom, 0.5);
112        super::normal::sample_unchecked(
113            r,
114            self.location,
115            self.scale * (self.freedom / gamma).sqrt(),
116        )
117    }
118}
119
120impl ContinuousCDF<f64, f64> for StudentsT {
121    /// Calculates the cumulative distribution function for the student's
122    /// t-distribution
123    /// at `x`
124    ///
125    /// # Formula
126    ///
127    /// ```ignore
128    /// if x < μ {
129    ///     (1 / 2) * I(t, v / 2, 1 / 2)
130    /// } else {
131    ///     1 - (1 / 2) * I(t, v / 2, 1 / 2)
132    /// }
133    /// ```
134    ///
135    /// where `t = v / (v + k^2)`, `k = (x - μ) / σ`, `μ` is the location,
136    /// `σ` is the scale, `v` is the freedom, and `I` is the regularized
137    /// incomplete beta function
138    fn cdf(&self, x: f64) -> f64 {
139        if self.freedom.is_infinite() {
140            super::normal::cdf_unchecked(x, self.location, self.scale)
141        } else {
142            let k = (x - self.location) / self.scale;
143            let h = self.freedom / (self.freedom + k * k);
144            let ib = 0.5 * beta::beta_reg(self.freedom / 2.0, 0.5, h);
145            if x <= self.location {
146                ib
147            } else {
148                1.0 - ib
149            }
150        }
151    }
152
153    /// Calculates the cumulative distribution function for the student's
154    /// t-distribution
155    /// at `x`
156    ///
157    /// # Formula
158    ///
159    /// ```ignore
160    /// if x < μ {
161    ///     1 - (1 / 2) * I(t, v / 2, 1 / 2)
162    /// } else {
163    ///     (1 / 2) * I(t, v / 2, 1 / 2)
164    /// }
165    /// ```
166    ///
167    /// where `t = v / (v + k^2)`, `k = (x - μ) / σ`, `μ` is the location,
168    /// `σ` is the scale, `v` is the freedom, and `I` is the regularized
169    /// incomplete beta function
170    fn sf(&self, x: f64) -> f64 {
171        if self.freedom.is_infinite() {
172            super::normal::sf_unchecked(x, self.location, self.scale)
173        } else {
174            let k = (x - self.location) / self.scale;
175            let h = self.freedom / (self.freedom + k * k);
176            let ib = 0.5 * beta::beta_reg(self.freedom / 2.0, 0.5, h);
177            if x <= self.location {
178                1.0 - ib
179            } else {
180                ib
181            }
182        }
183    }
184
185    /// Calculates the inverse cumulative distribution function for the
186    /// Student's T-distribution at `x`
187    fn inverse_cdf(&self, x: f64) -> f64 {
188        // first calculate inverse_cdf for normal Student's T
189        assert!((0.0..=1.0).contains(&x));
190        let x1 = if x >= 0.5 { 1.0 - x } else { x };
191        let a = 0.5 * self.freedom;
192        let b = 0.5;
193        let mut y = beta::inv_beta_reg(a, b, 2.0 * x1);
194        y = (self.freedom * (1. - y) / y).sqrt();
195        y = if x >= 0.5 { y } else { -y };
196        // generalised Student's T is related to normal Student's T by `Y = μ + σ X`
197        // where `X` is distributed as Student's T, so this result has to be scaled and shifted back
198        // formally: F_Y(t) = P(Y <= t) = P(X <= (t - μ) / σ) = F_X((t - μ) / σ)
199        // F_Y^{-1}(p) = inf { t' | F_Y(t') >= p } = inf { t' = μ + σ t | F_X((t' - μ) / σ) >= p }
200        // because scale is positive: loc + scale * t is strictly monotonic function
201        // = μ + σ inf { t | F_X(t) >= p } = μ + σ F_X^{-1}(p)
202        self.location + self.scale * y
203    }
204}
205
206impl Min<f64> for StudentsT {
207    /// Returns the minimum value in the domain of the student's t-distribution
208    /// representable by a double precision float
209    ///
210    /// # Formula
211    ///
212    /// ```ignore
213    /// -INF
214    /// ```
215    fn min(&self) -> f64 {
216        f64::NEG_INFINITY
217    }
218}
219
220impl Max<f64> for StudentsT {
221    /// Returns the maximum value in the domain of the student's t-distribution
222    /// representable by a double precision float
223    ///
224    /// # Formula
225    ///
226    /// ```ignore
227    /// INF
228    /// ```
229    fn max(&self) -> f64 {
230        f64::INFINITY
231    }
232}
233
234impl Distribution<f64> for StudentsT {
235    /// Returns the mean of the student's t-distribution
236    ///
237    /// # None
238    ///
239    /// If `freedom <= 1.0`
240    ///
241    /// # Formula
242    ///
243    /// ```ignore
244    /// μ
245    /// ```
246    ///
247    /// where `μ` is the location
248    fn mean(&self) -> Option<f64> {
249        if self.freedom <= 1.0 {
250            None
251        } else {
252            Some(self.location)
253        }
254    }
255    /// Returns the variance of the student's t-distribution
256    ///
257    /// # None
258    ///
259    /// If `freedom <= 2.0`
260    ///
261    /// # Formula
262    ///
263    /// ```ignore
264    /// if v == INF {
265    ///     Some(σ^2)
266    /// } else if freedom > 2.0 {
267    ///     Some(v * σ^2 / (v - 2))
268    /// } else {
269    ///     None
270    /// }
271    /// ```
272    ///
273    /// where `σ` is the scale and `v` is the freedom
274    fn variance(&self) -> Option<f64> {
275        if self.freedom.is_infinite() {
276            Some(self.scale * self.scale)
277        } else if self.freedom > 2.0 {
278            Some(self.freedom * self.scale * self.scale / (self.freedom - 2.0))
279        } else {
280            None
281        }
282    }
283    /// Returns the entropy for the student's t-distribution
284    ///
285    /// # Formula
286    ///
287    /// ```ignore
288    /// - ln(σ) + (v + 1) / 2 * (ψ((v + 1) / 2) - ψ(v / 2)) + ln(sqrt(v) * B(v / 2, 1 /
289    /// 2))
290    /// ```
291    ///
292    /// where `σ` is the scale, `v` is the freedom, `ψ` is the digamma function, and `B` is the
293    /// beta function
294    fn entropy(&self) -> Option<f64> {
295        // generalised Student's T is related to normal Student's T by `Y = μ + σ X`
296        // where `X` is distributed as Student's T, plugging into the definition
297        // of entropy shows scaling affects the entropy by an additive constant `- ln σ`
298        let shift = -self.scale.ln();
299        let result = (self.freedom + 1.0) / 2.0
300            * (gamma::digamma((self.freedom + 1.0) / 2.0) - gamma::digamma(self.freedom / 2.0))
301            + (self.freedom.sqrt() * beta::beta(self.freedom / 2.0, 0.5)).ln();
302        Some(result + shift)
303    }
304    /// Returns the skewness of the student's t-distribution
305    ///
306    /// # None
307    ///
308    /// If `x <= 3.0`
309    ///
310    /// # Formula
311    ///
312    /// ```ignore
313    /// 0
314    /// ```
315    fn skewness(&self) -> Option<f64> {
316        if self.freedom <= 3.0 {
317            None
318        } else {
319            Some(0.0)
320        }
321    }
322}
323
324impl Median<f64> for StudentsT {
325    /// Returns the median of the student's t-distribution
326    ///
327    /// # Formula
328    ///
329    /// ```ignore
330    /// μ
331    /// ```
332    ///
333    /// where `μ` is the location
334    fn median(&self) -> f64 {
335        self.location
336    }
337}
338
339impl Mode<Option<f64>> for StudentsT {
340    /// Returns the mode of the student's t-distribution
341    ///
342    /// # Formula
343    ///
344    /// ```ignore
345    /// μ
346    /// ```
347    ///
348    /// where `μ` is the location
349    fn mode(&self) -> Option<f64> {
350        Some(self.location)
351    }
352}
353
354impl Continuous<f64, f64> for StudentsT {
355    /// Calculates the probability density function for the student's
356    /// t-distribution
357    /// at `x`
358    ///
359    /// # Formula
360    ///
361    /// ```ignore
362    /// Γ((v + 1) / 2) / (sqrt(vπ) * Γ(v / 2) * σ) * (1 + k^2 / v)^(-1 / 2 * (v
363    /// + 1))
364    /// ```
365    ///
366    /// where `k = (x - μ) / σ`, `μ` is the location, `σ` is the scale, `v` is
367    /// the freedom,
368    /// and `Γ` is the gamma function
369    fn pdf(&self, x: f64) -> f64 {
370        if x.is_infinite() {
371            0.0
372        } else if self.freedom >= 1e8 {
373            super::normal::pdf_unchecked(x, self.location, self.scale)
374        } else {
375            let d = (x - self.location) / self.scale;
376            (gamma::ln_gamma((self.freedom + 1.0) / 2.0) - gamma::ln_gamma(self.freedom / 2.0))
377                .exp()
378                * (1.0 + d * d / self.freedom).powf(-0.5 * (self.freedom + 1.0))
379                / (self.freedom * f64::consts::PI).sqrt()
380                / self.scale
381        }
382    }
383
384    /// Calculates the log probability density function for the student's
385    /// t-distribution
386    /// at `x`
387    ///
388    /// # Formula
389    ///
390    /// ```ignore
391    /// ln(Γ((v + 1) / 2) / (sqrt(vπ) * Γ(v / 2) * σ) * (1 + k^2 / v)^(-1 / 2 *
392    /// (v + 1)))
393    /// ```
394    ///
395    /// where `k = (x - μ) / σ`, `μ` is the location, `σ` is the scale, `v` is
396    /// the freedom,
397    /// and `Γ` is the gamma function
398    fn ln_pdf(&self, x: f64) -> f64 {
399        if x.is_infinite() {
400            f64::NEG_INFINITY
401        } else if self.freedom >= 1e8 {
402            super::normal::ln_pdf_unchecked(x, self.location, self.scale)
403        } else {
404            let d = (x - self.location) / self.scale;
405            gamma::ln_gamma((self.freedom + 1.0) / 2.0)
406                - 0.5 * ((self.freedom + 1.0) * (1.0 + d * d / self.freedom).ln())
407                - gamma::ln_gamma(self.freedom / 2.0)
408                - 0.5 * (self.freedom * f64::consts::PI).ln()
409                - self.scale.ln()
410        }
411    }
412}
413
414#[cfg(all(test, feature = "nightly"))]
415mod tests {
416    use crate::consts::ACC;
417    use crate::distribution::internal::*;
418    use crate::distribution::{Continuous, ContinuousCDF, StudentsT};
419    use crate::statistics::*;
420    use crate::testing_boiler;
421    use std::panic;
422
423    testing_boiler!((f64, f64, f64), StudentsT);
424
425    #[test]
426    fn test_create() {
427        try_create((0.0, 0.1, 1.0));
428        try_create((0.0, 1.0, 1.0));
429        try_create((-5.0, 1.0, 3.0));
430        try_create((10.0, 10.0, f64::INFINITY));
431    }
432
433    // #[test]
434    // fn foo() {
435    //     let dist = StudentsT::new(0.0,1.0,1.0).unwrap();
436    //     dbg!(dist.mean());
437    // }
438
439    #[test]
440    fn test_bad_create() {
441        bad_create_case((f64::NAN, 1.0, 1.0));
442        bad_create_case((0.0, f64::NAN, 1.0));
443        bad_create_case((0.0, 1.0, f64::NAN));
444        bad_create_case((0.0, -10.0, 1.0));
445        bad_create_case((0.0, 10.0, -1.0));
446    }
447
448    #[test]
449    fn test_mean() {
450        let mean = |x: StudentsT| x.mean().unwrap();
451        test_case((0.0, 1.0, 3.0), 0.0, mean);
452        test_case((0.0, 10.0, 2.0), 0.0, mean);
453        test_case((0.0, 10.0, f64::INFINITY), 0.0, mean);
454        test_case((-5.0, 100.0, 1.5), -5.0, mean);
455        let mean = |x: StudentsT| x.mean();
456        test_none((0.0, 1.0, 1.0), mean);
457        test_none((0.0, 0.1, 1.0), mean);
458        test_none((0.0, 10.0, 1.0), mean);
459        test_none((10.0, 1.0, 1.0), mean);
460        test_none((0.0, f64::INFINITY, 1.0), mean);
461    }
462
463    #[test]
464    #[should_panic]
465    fn test_mean_freedom_lte_1() {
466        let mean = |x: StudentsT| x.mean().unwrap();
467        get_value((1.0, 1.0, 0.5), mean);
468    }
469
470    #[test]
471    fn test_variance() {
472        let variance = |x: StudentsT| x.variance().unwrap();
473        test_case((0.0, 1.0, 3.0), 3.0, variance);
474        test_case((0.0, 10.0, 2.5), 500.0, variance);
475        test_case((10.0, 1.0, 2.5), 5.0, variance);
476        let variance = |x: StudentsT| x.variance();
477        test_none((0.0, 10.0, 2.0), variance);
478        test_none((0.0, 1.0, 1.0), variance);
479        test_none((0.0, 0.1, 1.0), variance);
480        test_none((0.0, 10.0, 1.0), variance);
481        test_none((10.0, 1.0, 1.0), variance);
482        test_none((-5.0, 100.0, 1.5), variance);
483        test_none((0.0, f64::INFINITY, 1.0), variance);
484    }
485
486    #[test]
487    #[should_panic]
488    fn test_variance_freedom_lte1() {
489        let variance = |x: StudentsT| x.variance().unwrap();
490        get_value((1.0, 1.0, 0.5), variance);
491    }
492
493    // TODO: valid skewness tests
494    #[test]
495    #[should_panic]
496    fn test_skewness_freedom_lte_3() {
497        let skewness = |x: StudentsT| x.skewness().unwrap();
498        get_value((1.0, 1.0, 1.0), skewness);
499    }
500
501    #[test]
502    fn test_mode() {
503        let mode = |x: StudentsT| x.mode().unwrap();
504        test_case((0.0, 1.0, 1.0), 0.0, mode);
505        test_case((0.0, 0.1, 1.0), 0.0, mode);
506        test_case((0.0, 1.0, 3.0), 0.0, mode);
507        test_case((0.0, 10.0, 1.0), 0.0, mode);
508        test_case((0.0, 10.0, 2.0), 0.0, mode);
509        test_case((0.0, 10.0, 2.5), 0.0, mode);
510        test_case((0.0, 10.0, f64::INFINITY), 0.0, mode);
511        test_case((10.0, 1.0, 1.0), 10.0, mode);
512        test_case((10.0, 1.0, 2.5), 10.0, mode);
513        test_case((-5.0, 100.0, 1.5), -5.0, mode);
514        test_case((0.0, f64::INFINITY, 1.0), 0.0, mode);
515    }
516
517    #[test]
518    fn test_median() {
519        let median = |x: StudentsT| x.median();
520        test_case((0.0, 1.0, 1.0), 0.0, median);
521        test_case((0.0, 0.1, 1.0), 0.0, median);
522        test_case((0.0, 1.0, 3.0), 0.0, median);
523        test_case((0.0, 10.0, 1.0), 0.0, median);
524        test_case((0.0, 10.0, 2.0), 0.0, median);
525        test_case((0.0, 10.0, 2.5), 0.0, median);
526        test_case((0.0, 10.0, f64::INFINITY), 0.0, median);
527        test_case((10.0, 1.0, 1.0), 10.0, median);
528        test_case((10.0, 1.0, 2.5), 10.0, median);
529        test_case((-5.0, 100.0, 1.5), -5.0, median);
530        test_case((0.0, f64::INFINITY, 1.0), 0.0, median);
531    }
532
533    #[test]
534    fn test_min_max() {
535        let min = |x: StudentsT| x.min();
536        let max = |x: StudentsT| x.max();
537        test_case((0.0, 1.0, 1.0), f64::NEG_INFINITY, min);
538        test_case((2.5, 100.0, 1.5), f64::NEG_INFINITY, min);
539        test_case((10.0, f64::INFINITY, 3.5), f64::NEG_INFINITY, min);
540        test_case((0.0, 1.0, 1.0), f64::INFINITY, max);
541        test_case((2.5, 100.0, 1.5), f64::INFINITY, max);
542        test_case((10.0, f64::INFINITY, 5.5), f64::INFINITY, max);
543    }
544
545    #[test]
546    fn test_pdf() {
547        let pdf = |arg: f64| move |x: StudentsT| x.pdf(arg);
548        test_case((0.0, 1.0, 1.0), 0.318309886183791, pdf(0.0));
549        test_case((0.0, 1.0, 1.0), 0.159154943091895, pdf(1.0));
550        test_case((0.0, 1.0, 1.0), 0.159154943091895, pdf(-1.0));
551        test_case((0.0, 1.0, 1.0), 0.063661977236758, pdf(2.0));
552        test_case((0.0, 1.0, 1.0), 0.063661977236758, pdf(-2.0));
553        test_case((0.0, 1.0, 2.0), 0.353553390593274, pdf(0.0));
554        test_case((0.0, 1.0, 2.0), 0.192450089729875, pdf(1.0));
555        test_case((0.0, 1.0, 2.0), 0.192450089729875, pdf(-1.0));
556        test_case((0.0, 1.0, 2.0), 0.068041381743977, pdf(2.0));
557        test_case((0.0, 1.0, 2.0), 0.068041381743977, pdf(-2.0));
558        test_case((0.0, 1.0, f64::INFINITY), 0.398942280401433, pdf(0.0));
559        test_case((0.0, 1.0, f64::INFINITY), 0.241970724519143, pdf(1.0));
560        test_case((0.0, 1.0, f64::INFINITY), 0.053990966513188, pdf(2.0));
561    }
562
563    #[test]
564    fn test_ln_pdf() {
565        let ln_pdf = |arg: f64| move |x: StudentsT| x.ln_pdf(arg);
566        test_case((0.0, 1.0, 1.0), -1.144729885849399, ln_pdf(0.0));
567        test_case((0.0, 1.0, 1.0), -1.837877066409348, ln_pdf(1.0));
568        test_case((0.0, 1.0, 1.0), -1.837877066409348, ln_pdf(-1.0));
569        test_case((0.0, 1.0, 1.0), -2.754167798283503, ln_pdf(2.0));
570        test_case((0.0, 1.0, 1.0), -2.754167798283503, ln_pdf(-2.0));
571        test_case((0.0, 1.0, 2.0), -1.039720770839917, ln_pdf(0.0));
572        test_case((0.0, 1.0, 2.0), -1.647918433002166, ln_pdf(1.0));
573        test_case((0.0, 1.0, 2.0), -1.647918433002166, ln_pdf(-1.0));
574        test_case((0.0, 1.0, 2.0), -2.687639203842085, ln_pdf(2.0));
575        test_case((0.0, 1.0, 2.0), -2.687639203842085, ln_pdf(-2.0));
576        test_case((0.0, 1.0, f64::INFINITY), -0.918938533204672, ln_pdf(0.0));
577        test_case((0.0, 1.0, f64::INFINITY), -1.418938533204674, ln_pdf(1.0));
578        test_case((0.0, 1.0, f64::INFINITY), -2.918938533204674, ln_pdf(2.0));
579    }
580
581    #[test]
582    fn test_cdf() {
583        let cdf = |arg: f64| move |x: StudentsT| x.cdf(arg);
584        test_case((0.0, 1.0, 1.0), 0.5, cdf(0.0));
585        test_case((0.0, 1.0, 1.0), 0.75, cdf(1.0));
586        test_case((0.0, 1.0, 1.0), 0.25, cdf(-1.0));
587        test_case((0.0, 1.0, 1.0), 0.852416382349567, cdf(2.0));
588        test_case((0.0, 1.0, 1.0), 0.147583617650433, cdf(-2.0));
589        test_case((0.0, 1.0, 2.0), 0.5, cdf(0.0));
590        test_case((0.0, 1.0, 2.0), 0.788675134594813, cdf(1.0));
591        test_case((0.0, 1.0, 2.0), 0.211324865405187, cdf(-1.0));
592        test_case((0.0, 1.0, 2.0), 0.908248290463863, cdf(2.0));
593        test_case((0.0, 1.0, 2.0), 0.091751709536137, cdf(-2.0));
594        test_case((0.0, 1.0, f64::INFINITY), 0.5, cdf(0.0));
595
596        // TODO: these are curiously low accuracy and should be re-examined
597        test_case((0.0, 1.0, f64::INFINITY), 0.841344746068543, cdf(1.0));
598        test_case((0.0, 1.0, f64::INFINITY), 0.977249868051821, cdf(2.0));
599    }
600
601
602    #[test]
603    fn test_sf() {
604        let sf = |arg: f64| move |x: StudentsT| x.sf(arg);
605        test_case((0.0, 1.0, 1.0), 0.5, sf(0.0));
606        test_case((0.0, 1.0, 1.0), 0.25, sf(1.0));
607        test_case((0.0, 1.0, 1.0), 0.75, sf(-1.0));
608        test_case((0.0, 1.0, 1.0), 0.147583617650433, sf(2.0));
609        test_case((0.0, 1.0, 1.0), 0.852416382349566, sf(-2.0));
610        test_case((0.0, 1.0, 2.0), 0.5, sf(0.0));
611        test_case((0.0, 1.0, 2.0), 0.211324865405186, sf(1.0));
612        test_case((0.0, 1.0, 2.0), 0.788675134594813, sf(-1.0));
613        test_case((0.0, 1.0, 2.0), 0.091751709536137, sf(2.0));
614        test_case((0.0, 1.0, 2.0), 0.908248290463862, sf(-2.0));
615        test_case((0.0, 1.0, f64::INFINITY), 0.5, sf(0.0));
616
617        // TODO: these are curiously low accuracy and should be re-examined
618        test_case((0.0, 1.0, f64::INFINITY), 0.158655253945057, sf(1.0));
619        test_case((0.0, 1.0, f64::INFINITY), 0.022750131947162, sf(2.0));
620    }
621
622    #[test]
623    fn test_continuous() {
624        test::check_continuous_distribution(&try_create((0.0, 1.0, 3.0)), -30.0, 30.0);
625        test::check_continuous_distribution(&try_create((0.0, 1.0, 10.0)), -10.0, 10.0);
626        test::check_continuous_distribution(&try_create((20.0, 0.5, 10.0)), 10.0, 30.0);
627    }
628
629    #[test]
630    fn test_inv_cdf() {
631        let test = |x: f64, freedom: f64, expected: f64| {
632            use approx::*;
633            let d = StudentsT::new(0., 1., freedom).unwrap();
634            // Checks that left == right to 4 significant figures, unlike
635            // test_almost() which uses decimal places
636            assert_relative_eq!(d.inverse_cdf(x), expected, max_relative = 0.001);
637        };
638
639        // This test checks our implementation against the whole t-table
640        // copied from https://en.wikipedia.org/wiki/Student's_t-distribution
641
642        test(0.75, 1.0, 1.000);
643        test(0.8, 1.0, 1.376);
644        test(0.85, 1.0, 1.963);
645        test(0.9, 1.0, 3.078);
646        test(0.95, 1.0, 6.314);
647        test(0.975, 1.0, 12.71);
648        test(0.99, 1.0, 31.82);
649        test(0.995, 1.0, 63.66);
650        test(0.9975, 1.0, 127.3);
651        test(0.999, 1.0, 318.3);
652        test(0.9995, 1.0, 636.6);
653
654        test(0.75, 002.0, 0.816);
655        // TODO: investigate
656        // test(0.8, 002.0, 1.080);  // We get 1.061 for some reason...
657        test(0.85, 002.0, 1.386);
658        test(0.9, 002.0, 1.886);
659        test(0.95, 002.0, 2.920);
660        test(0.975, 002.0, 4.303);
661        test(0.99, 002.0, 6.965);
662        test(0.995, 002.0, 9.925);
663        test(0.9975, 002.0, 14.09);
664        test(0.999, 002.0, 22.33);
665        test(0.9995, 002.0, 31.60);
666
667        test(0.75, 003.0, 0.765);
668        test(0.8, 003.0, 0.978);
669        test(0.85, 003.0, 1.250);
670        test(0.9, 003.0, 1.638);
671        test(0.95, 003.0, 2.353);
672        test(0.975, 003.0, 3.182);
673        test(0.99, 003.0, 4.541);
674        test(0.995, 003.0, 5.841);
675        test(0.9975, 003.0, 7.453);
676        test(0.999, 003.0, 10.21);
677        test(0.9995, 003.0, 12.92);
678
679        test(0.75, 004.0, 0.741);
680        test(0.8, 004.0, 0.941);
681        test(0.85, 004.0, 1.190);
682        test(0.9, 004.0, 1.533);
683        test(0.95, 004.0, 2.132);
684        test(0.975, 004.0, 2.776);
685        test(0.99, 004.0, 3.747);
686        test(0.995, 004.0, 4.604);
687        test(0.9975, 004.0, 5.598);
688        test(0.999, 004.0, 7.173);
689        test(0.9995, 004.0, 8.610);
690
691        test(0.75, 005.0, 0.727);
692        test(0.8, 005.0, 0.920);
693        test(0.85, 005.0, 1.156);
694        test(0.9, 005.0, 1.476);
695        test(0.95, 005.0, 2.015);
696        test(0.975, 005.0, 2.571);
697        test(0.99, 005.0, 3.365);
698        test(0.995, 005.0, 4.032);
699        test(0.9975, 005.0, 4.773);
700        test(0.999, 005.0, 5.893);
701        test(0.9995, 005.0, 6.869);
702
703        test(0.75, 006.0, 0.718);
704        test(0.8, 006.0, 0.906);
705        test(0.85, 006.0, 1.134);
706        test(0.9, 006.0, 1.440);
707        test(0.95, 006.0, 1.943);
708        test(0.975, 006.0, 2.447);
709        test(0.99, 006.0, 3.143);
710        test(0.995, 006.0, 3.707);
711        test(0.9975, 006.0, 4.317);
712        test(0.999, 006.0, 5.208);
713        test(0.9995, 006.0, 5.959);
714
715        test(0.75, 007.0, 0.711);
716        test(0.8, 007.0, 0.896);
717        test(0.85, 007.0, 1.119);
718        test(0.9, 007.0, 1.415);
719        test(0.95, 007.0, 1.895);
720        test(0.975, 007.0, 2.365);
721        test(0.99, 007.0, 2.998);
722        test(0.995, 007.0, 3.499);
723        test(0.9975, 007.0, 4.029);
724        test(0.999, 007.0, 4.785);
725        test(0.9995, 007.0, 5.408);
726
727        test(0.75, 008.0, 0.706);
728        test(0.8, 008.0, 0.889);
729        test(0.85, 008.0, 1.108);
730        test(0.9, 008.0, 1.397);
731        test(0.95, 008.0, 1.860);
732        test(0.975, 008.0, 2.306);
733        test(0.99, 008.0, 2.896);
734        test(0.995, 008.0, 3.355);
735        test(0.9975, 008.0, 3.833);
736        test(0.999, 008.0, 4.501);
737        test(0.9995, 008.0, 5.041);
738
739        test(0.75, 009.0, 0.703);
740        test(0.8, 009.0, 0.883);
741        test(0.85, 009.0, 1.100);
742        test(0.9, 009.0, 1.383);
743        test(0.95, 009.0, 1.833);
744        test(0.975, 009.0, 2.262);
745        test(0.99, 009.0, 2.821);
746        test(0.995, 009.0, 3.250);
747        test(0.9975, 009.0, 3.690);
748        test(0.999, 009.0, 4.297);
749        test(0.9995, 009.0, 4.781);
750
751        test(0.75, 010.0, 0.700);
752        test(0.8, 010.0, 0.879);
753        test(0.85, 010.0, 1.093);
754        test(0.9, 010.0, 1.372);
755        test(0.95, 010.0, 1.812);
756        test(0.975, 010.0, 2.228);
757        test(0.99, 010.0, 2.764);
758        test(0.995, 010.0, 3.169);
759        test(0.9975, 010.0, 3.581);
760        test(0.999, 010.0, 4.144);
761        test(0.9995, 010.0, 4.587);
762
763        test(0.75, 011.0, 0.697);
764        test(0.8, 011.0, 0.876);
765        test(0.85, 011.0, 1.088);
766        test(0.9, 011.0, 1.363);
767        test(0.95, 011.0, 1.796);
768        test(0.975, 011.0, 2.201);
769        test(0.99, 011.0, 2.718);
770        test(0.995, 011.0, 3.106);
771        test(0.9975, 011.0, 3.497);
772        test(0.999, 011.0, 4.025);
773        test(0.9995, 011.0, 4.437);
774
775        test(0.75, 012.0, 0.695);
776        test(0.8, 012.0, 0.873);
777        test(0.85, 012.0, 1.083);
778        test(0.9, 012.0, 1.356);
779        test(0.95, 012.0, 1.782);
780        test(0.975, 012.0, 2.179);
781        test(0.99, 012.0, 2.681);
782        test(0.995, 012.0, 3.055);
783        test(0.9975, 012.0, 3.428);
784        test(0.999, 012.0, 3.930);
785        test(0.9995, 012.0, 4.318);
786
787        test(0.75, 013.0, 0.694);
788        test(0.8, 013.0, 0.870);
789        test(0.85, 013.0, 1.079);
790        test(0.9, 013.0, 1.350);
791        test(0.95, 013.0, 1.771);
792        test(0.975, 013.0, 2.160);
793        test(0.99, 013.0, 2.650);
794        test(0.995, 013.0, 3.012);
795        test(0.9975, 013.0, 3.372);
796        test(0.999, 013.0, 3.852);
797        test(0.9995, 013.0, 4.221);
798
799        test(0.75, 014.0, 0.692);
800        test(0.8, 014.0, 0.868);
801        test(0.85, 014.0, 1.076);
802        test(0.9, 014.0, 1.345);
803        test(0.95, 014.0, 1.761);
804        test(0.975, 014.0, 2.145);
805        test(0.99, 014.0, 2.624);
806        test(0.995, 014.0, 2.977);
807        test(0.9975, 014.0, 3.326);
808        test(0.999, 014.0, 3.787);
809        test(0.9995, 014.0, 4.140);
810
811        test(0.75, 015.0, 0.691);
812        test(0.8, 015.0, 0.866);
813        test(0.85, 015.0, 1.074);
814        test(0.9, 015.0, 1.341);
815        test(0.95, 015.0, 1.753);
816        test(0.975, 015.0, 2.131);
817        test(0.99, 015.0, 2.602);
818        test(0.995, 015.0, 2.947);
819        test(0.9975, 015.0, 3.286);
820        test(0.999, 015.0, 3.733);
821        test(0.9995, 015.0, 4.073);
822
823        test(0.75, 016.0, 0.690);
824        test(0.8, 016.0, 0.865);
825        test(0.85, 016.0, 1.071);
826        test(0.9, 016.0, 1.337);
827        test(0.95, 016.0, 1.746);
828        test(0.975, 016.0, 2.120);
829        test(0.99, 016.0, 2.583);
830        test(0.995, 016.0, 2.921);
831        test(0.9975, 016.0, 3.252);
832        test(0.999, 016.0, 3.686);
833        test(0.9995, 016.0, 4.015);
834
835        test(0.75, 017.0, 0.689);
836        test(0.8, 017.0, 0.863);
837        test(0.85, 017.0, 1.069);
838        test(0.9, 017.0, 1.333);
839        test(0.95, 017.0, 1.740);
840        test(0.975, 017.0, 2.110);
841        test(0.99, 017.0, 2.567);
842        test(0.995, 017.0, 2.898);
843        test(0.9975, 017.0, 3.222);
844        test(0.999, 017.0, 3.646);
845        test(0.9995, 017.0, 3.965);
846
847        test(0.75, 018.0, 0.688);
848        test(0.8, 018.0, 0.862);
849        test(0.85, 018.0, 1.067);
850        test(0.9, 018.0, 1.330);
851        test(0.95, 018.0, 1.734);
852        test(0.975, 018.0, 2.101);
853        test(0.99, 018.0, 2.552);
854        test(0.995, 018.0, 2.878);
855        test(0.9975, 018.0, 3.197);
856        test(0.999, 018.0, 3.610);
857        test(0.9995, 018.0, 3.922);
858
859        test(0.75, 019.0, 0.688);
860        test(0.8, 019.0, 0.861);
861        test(0.85, 019.0, 1.066);
862        test(0.9, 019.0, 1.328);
863        test(0.95, 019.0, 1.729);
864        test(0.975, 019.0, 2.093);
865        test(0.99, 019.0, 2.539);
866        test(0.995, 019.0, 2.861);
867        test(0.9975, 019.0, 3.174);
868        test(0.999, 019.0, 3.579);
869        test(0.9995, 019.0, 3.883);
870
871        test(0.75, 020.0, 0.687);
872        test(0.8, 020.0, 0.860);
873        test(0.85, 020.0, 1.064);
874        test(0.9, 020.0, 1.325);
875        test(0.95, 020.0, 1.725);
876        test(0.975, 020.0, 2.086);
877        test(0.99, 020.0, 2.528);
878        test(0.995, 020.0, 2.845);
879        test(0.9975, 020.0, 3.153);
880        test(0.999, 020.0, 3.552);
881        test(0.9995, 020.0, 3.850);
882
883        test(0.75, 021.0, 0.686);
884        test(0.8, 021.0, 0.859);
885        test(0.85, 021.0, 1.063);
886        test(0.9, 021.0, 1.323);
887        test(0.95, 021.0, 1.721);
888        test(0.975, 021.0, 2.080);
889        test(0.99, 021.0, 2.518);
890        test(0.995, 021.0, 2.831);
891        test(0.9975, 021.0, 3.135);
892        test(0.999, 021.0, 3.527);
893        test(0.9995, 021.0, 3.819);
894
895        test(0.75, 022.0, 0.686);
896        test(0.8, 022.0, 0.858);
897        test(0.85, 022.0, 1.061);
898        test(0.9, 022.0, 1.321);
899        test(0.95, 022.0, 1.717);
900        test(0.975, 022.0, 2.074);
901        test(0.99, 022.0, 2.508);
902        test(0.995, 022.0, 2.819);
903        test(0.9975, 022.0, 3.119);
904        test(0.999, 022.0, 3.505);
905        test(0.9995, 022.0, 3.792);
906
907        test(0.75, 023.0, 0.685);
908        test(0.8, 023.0, 0.858);
909        test(0.85, 023.0, 1.060);
910        test(0.9, 023.0, 1.319);
911        test(0.95, 023.0, 1.714);
912        test(0.975, 023.0, 2.069);
913        test(0.99, 023.0, 2.500);
914        test(0.995, 023.0, 2.807);
915        test(0.9975, 023.0, 3.104);
916        test(0.999, 023.0, 3.485);
917        test(0.9995, 023.0, 3.767);
918
919        test(0.75, 024.0, 0.685);
920        test(0.8, 024.0, 0.857);
921        test(0.85, 024.0, 1.059);
922        test(0.9, 024.0, 1.318);
923        test(0.95, 024.0, 1.711);
924        test(0.975, 024.0, 2.064);
925        test(0.99, 024.0, 2.492);
926        test(0.995, 024.0, 2.797);
927        test(0.9975, 024.0, 3.091);
928        test(0.999, 024.0, 3.467);
929        test(0.9995, 024.0, 3.745);
930
931        test(0.75, 025.0, 0.684);
932        test(0.8, 025.0, 0.856);
933        test(0.85, 025.0, 1.058);
934        test(0.9, 025.0, 1.316);
935        test(0.95, 025.0, 1.708);
936        test(0.975, 025.0, 2.060);
937        test(0.99, 025.0, 2.485);
938        test(0.995, 025.0, 2.787);
939        test(0.9975, 025.0, 3.078);
940        test(0.999, 025.0, 3.450);
941        test(0.9995, 025.0, 3.725);
942
943        test(0.75, 026.0, 0.684);
944        test(0.8, 026.0, 0.856);
945        test(0.85, 026.0, 1.058);
946        test(0.9, 026.0, 1.315);
947        test(0.95, 026.0, 1.706);
948        test(0.975, 026.0, 2.056);
949        test(0.99, 026.0, 2.479);
950        test(0.995, 026.0, 2.779);
951        test(0.9975, 026.0, 3.067);
952        test(0.999, 026.0, 3.435);
953        test(0.9995, 026.0, 3.707);
954
955        test(0.75, 027.0, 0.684);
956        test(0.8, 027.0, 0.855);
957        test(0.85, 027.0, 1.057);
958        test(0.9, 027.0, 1.314);
959        test(0.95, 027.0, 1.703);
960        test(0.975, 027.0, 2.052);
961        test(0.99, 027.0, 2.473);
962        test(0.995, 027.0, 2.771);
963        test(0.9975, 027.0, 3.057);
964        test(0.999, 027.0, 3.421);
965        test(0.9995, 027.0, 3.690);
966
967        test(0.75, 028.0, 0.683);
968        test(0.8, 028.0, 0.855);
969        test(0.85, 028.0, 1.056);
970        test(0.9, 028.0, 1.313);
971        test(0.95, 028.0, 1.701);
972        test(0.975, 028.0, 2.048);
973        test(0.99, 028.0, 2.467);
974        test(0.995, 028.0, 2.763);
975        test(0.9975, 028.0, 3.047);
976        test(0.999, 028.0, 3.408);
977        test(0.9995, 028.0, 3.674);
978
979        test(0.75, 029.0, 0.683);
980        test(0.8, 029.0, 0.854);
981        test(0.85, 029.0, 1.055);
982        test(0.9, 029.0, 1.311);
983        test(0.95, 029.0, 1.699);
984        test(0.975, 029.0, 2.045);
985        test(0.99, 029.0, 2.462);
986        test(0.995, 029.0, 2.756);
987        test(0.9975, 029.0, 3.038);
988        test(0.999, 029.0, 3.396);
989        test(0.9995, 029.0, 3.659);
990
991        test(0.75, 030.0, 0.683);
992        test(0.8, 030.0, 0.854);
993        test(0.85, 030.0, 1.055);
994        test(0.9, 030.0, 1.310);
995        test(0.95, 030.0, 1.697);
996        test(0.975, 030.0, 2.042);
997        test(0.99, 030.0, 2.457);
998        test(0.995, 030.0, 2.750);
999        test(0.9975, 030.0, 3.030);
1000        test(0.999, 030.0, 3.385);
1001        test(0.9995, 030.0, 3.646);
1002
1003        test(0.75, 040.0, 0.681);
1004        test(0.8, 040.0, 0.851);
1005        test(0.85, 040.0, 1.050);
1006        test(0.9, 040.0, 1.303);
1007        test(0.95, 040.0, 1.684);
1008        test(0.975, 040.0, 2.021);
1009        test(0.99, 040.0, 2.423);
1010        test(0.995, 040.0, 2.704);
1011        test(0.9975, 040.0, 2.971);
1012        test(0.999, 040.0, 3.307);
1013        test(0.9995, 040.0, 3.551);
1014
1015        test(0.75, 050.0, 0.679);
1016        test(0.8, 050.0, 0.849);
1017        test(0.85, 050.0, 1.047);
1018        test(0.9, 050.0, 1.299);
1019        test(0.95, 050.0, 1.676);
1020        test(0.975, 050.0, 2.009);
1021        test(0.99, 050.0, 2.403);
1022        test(0.995, 050.0, 2.678);
1023        test(0.9975, 050.0, 2.937);
1024        test(0.999, 050.0, 3.261);
1025        test(0.9995, 050.0, 3.496);
1026
1027        test(0.75, 060.0, 0.679);
1028        test(0.8, 060.0, 0.848);
1029        test(0.85, 060.0, 1.045);
1030        test(0.9, 060.0, 1.296);
1031        test(0.95, 060.0, 1.671);
1032        test(0.975, 060.0, 2.000);
1033        test(0.99, 060.0, 2.390);
1034        test(0.995, 060.0, 2.660);
1035        test(0.9975, 060.0, 2.915);
1036        test(0.999, 060.0, 3.232);
1037        test(0.9995, 060.0, 3.460);
1038
1039        test(0.75, 080.0, 0.678);
1040        test(0.8, 080.0, 0.846);
1041        test(0.85, 080.0, 1.043);
1042        test(0.9, 080.0, 1.292);
1043        test(0.95, 080.0, 1.664);
1044        test(0.975, 080.0, 1.990);
1045        test(0.99, 080.0, 2.374);
1046        test(0.995, 080.0, 2.639);
1047        test(0.9975, 080.0, 2.887);
1048        test(0.999, 080.0, 3.195);
1049        test(0.9995, 080.0, 3.416);
1050
1051        test(0.75, 100.0, 0.677);
1052        test(0.8, 100.0, 0.845);
1053        test(0.85, 100.0, 1.042);
1054        test(0.9, 100.0, 1.290);
1055        test(0.95, 100.0, 1.660);
1056        test(0.975, 100.0, 1.984);
1057        test(0.99, 100.0, 2.364);
1058        test(0.995, 100.0, 2.626);
1059        test(0.9975, 100.0, 2.871);
1060        test(0.999, 100.0, 3.174);
1061        test(0.9995, 100.0, 3.390);
1062
1063        test(0.75, 120.0, 0.677);
1064        test(0.8, 120.0, 0.845);
1065        test(0.85, 120.0, 1.041);
1066        test(0.9, 120.0, 1.289);
1067        test(0.95, 120.0, 1.658);
1068        test(0.975, 120.0, 1.980);
1069        test(0.99, 120.0, 2.358);
1070        test(0.995, 120.0, 2.617);
1071        test(0.9975, 120.0, 2.860);
1072        test(0.999, 120.0, 3.160);
1073        test(0.9995, 120.0, 3.373);
1074    }
1075
1076    #[test]
1077    fn test_inv_cdf_high_precision() {
1078        let test = |x: f64, freedom: f64, expected: f64| {
1079            use approx::assert_relative_eq;
1080            let d = StudentsT::new(0., 1., freedom).unwrap();
1081            assert_relative_eq!(d.inverse_cdf(x), expected, max_relative = ACC);
1082        };
1083        // The data in this table of expected values was generated in
1084        // Python, using the mpsci package (based on mpmath):
1085        //
1086        //   import mpmath
1087        //   from mpsci.distributions import t
1088        //
1089        //   # Set the number of digits of precision
1090        //   mpmath.mp.dps = 200
1091        //
1092        //   ps = [0.001, 0.01, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45]
1093        //   dfs = [1.0, 10.0, 100.0]
1094        //
1095        //   for df in dfs:
1096        //       for p in ps:
1097        //           q = t.invcdf(p, df)
1098        //           print(f"({p:5.3f}, {df:5.1f}, {float(q)}),")
1099        //
1100        let invcdf_data = [
1101            // p       df    inverse_cdf(p, df)
1102            (0.001,   1.0, -318.30883898555044),
1103            (0.010,   1.0, -31.820515953773956),
1104            (0.100,   1.0, -3.077683537175253),
1105            (0.150,   1.0, -1.9626105055051506),
1106            (0.200,   1.0, -1.3763819204711734),
1107            (0.250,   1.0, -1.0),
1108            (0.300,   1.0, -0.7265425280053609),
1109            (0.350,   1.0, -0.5095254494944289),
1110            (0.400,   1.0, -0.32491969623290623),
1111            (0.450,   1.0, -0.15838444032453625),
1112            (0.001,  10.0, -4.143700494046589),
1113            (0.010,  10.0, -2.763769458112696),
1114            (0.100,  10.0, -1.3721836411103356),
1115            (0.150,  10.0, -1.093058073590526),
1116            (0.200,  10.0, -0.8790578285505887),
1117            (0.250,  10.0, -0.6998120613124317),
1118            (0.300,  10.0, -0.5415280387550157),
1119            (0.350,  10.0, -0.3965914937556218),
1120            (0.400,  10.0, -0.26018482949208016),
1121            (0.450,  10.0, -0.12889018929327375),
1122            (0.001, 100.0, -3.173739493738783),
1123            (0.010, 100.0, -2.364217366238482),
1124            (0.100, 100.0, -1.290074761346516),
1125            (0.150, 100.0, -1.041835900908347),
1126            (0.200, 100.0, -0.845230424491016),
1127            (0.250, 100.0, -0.6769510430114715),
1128            (0.300, 100.0, -0.5260762706003463),
1129            (0.350, 100.0, -0.3864289804076715),
1130            (0.400, 100.0, -0.2540221824582278),
1131            (0.450, 100.0, -0.12598088204153965),
1132        ];
1133        for (p, df, expected) in invcdf_data.iter() {
1134            test(*p, *df, *expected);
1135            test(1.0 - *p, *df, -*expected);
1136        }
1137    }
1138
1139    #[test]
1140    fn test_inv_cdf_midpoint() {
1141        for loc in [0.0, 1.0, -3.5] {
1142            let d = StudentsT::new(loc, 1.0, 12.0).unwrap();
1143            // inverse_cdf(p) is a floating point calculation, so using
1144            // assert_eq here is optimistic.  For the given location values,
1145            // the check passes, so let's use the optimistic check for now.
1146            assert_eq!(d.inverse_cdf(0.5), loc);
1147        }
1148    }
1149
1150    #[test]
1151    fn test_inv_cdf_p0() {
1152        let d = StudentsT::new(0.0, 1.0, 12.0).unwrap();
1153        assert_eq!(d.inverse_cdf(0.0), std::f64::NEG_INFINITY);
1154    }
1155
1156    #[test]
1157    fn test_inv_cdf_p1() {
1158        let d = StudentsT::new(0.0, 1.0, 12.0).unwrap();
1159        assert_eq!(d.inverse_cdf(1.0), std::f64::INFINITY);
1160    }
1161}
statrs/distribution/students_t.rs

statrs/distribution/
students_t.rs
