statrs/distribution/
students_t.rs

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

statrs/distribution/
students_t.rs
