spec_math/cephes64/
stdtr.rs

1/*
2* Cephes Math Library Release 2.3:  March, 1995
3* Copyright 1984, 1987, 1995 by Stephen L. Moshier
4*/
5
6use crate::cephes64::consts::{M_PI, MACHEP};
7use crate::cephes64::incbet::incbet;
8use crate::cephes64::incbi::incbi;
9
10pub fn stdtr(k: isize, t: f64) -> f64 {
11    //! Student's t distribution
12    //!
13    //! ## DESCRIPTION:
14    //!
15    //! Computes the integral from minus infinity to t of the Student
16    //! t distribution with integer k > 0 degrees of freedom:
17    //!
18    #![doc=include_str!("stdtr.svg")]
19    //!
20    //! Relation to incomplete beta integral:
21    //!
22    //! `1 - stdtr(k,t) = 0.5 * incbet( k/2, 1/2, z )`
23    //!
24    //! where
25    //!
26    //! `z = k/(k + t**2)`
27    //!
28    //! For t < -2, this is the method of computation.  For higher t,
29    //! a direct method is derived from integration by parts.
30    //! Since the function is symmetric about t=0, the area under the
31    //! right tail of the density is found by calling the function
32    //! with -t instead of t.
33    //!
34    //! ## ACCURACY:
35    //!
36    //! Tested at random 1 <= k <= 25.  The "domain" refers to t.
37    //!
38    //! Relative error:
39    //!
40    //!<table>
41    //! <tr>
42    //!     <th>Arithmetic</th>
43    //!     <th>Domain</th>
44    //!     <th># Trials</th>
45    //!     <th>Peak</th>
46    //!     <th>RMS</th>
47    //! </tr>
48    //! <tr>
49    //!     <td>IEEE</td>
50    //!     <td>-100, -2</td>
51    //!     <td>50000</td>
52    //!     <td>5.9e-15</td>
53    //!     <td>1.4e-15</td>
54    //! </tr>
55    //! <tr>
56    //!     <td>IEEE</td>
57    //!     <td>-2, 100</td>
58    //!     <td>500000</td>
59    //!     <td>2.7e-15</td>
60    //!     <td>4.9e-17</td>
61    //! </tr>
62    //!</table>
63
64    if k <= 0 {
65        //sf_error("stdtr", SF_ERROR_DOMAIN, NULL);
66        f64::NAN
67    } else if t == 0.0 {
68        0.5
69    } else if t < -2.0 {
70        let rk = k as f64;
71        let z = rk / (rk + t * t);
72        0.5 * incbet(0.5 * rk, 0.5, z)
73    } else {
74        // compute integral from -t to + t
75
76        let x = t.abs();
77
78        let rk = k as f64; // degrees of freedom
79        let z = 1.0 + (x * x) / rk;
80
81        let mut p: f64;
82
83        // test if k is odd or even
84        if k & 1 != 0 {
85            // computation for odd k
86            let xsqk = x / rk.sqrt();
87            p = xsqk.atan();
88            if k > 1 {
89                let mut f = 1.0;
90                let mut tz = 1.0;
91                let mut j = 3;
92                while j <= k - 2 && tz / f > MACHEP {
93                    tz *= (j - 1) as f64 / (z * j as f64);
94                    f += tz;
95                    j += 2;
96                }
97                p += f * xsqk / z;
98            }
99            p *= 2.0 / M_PI;
100
101        } else {
102            // computation for even k
103
104            let mut f = 1.0;
105            let mut tz = 1.0;
106            let mut j = 2;
107
108            while j <= k - 2 && tz / f > MACHEP {
109                tz *= (j - 1) as f64 / (z * j as f64);
110                f += tz;
111                j += 2;
112            }
113            p = f * x / (z * rk).sqrt();
114        }
115
116        // common exit
117
118
119        if t < 0.0 {
120            p = -p; // note destruction of relative accuracy
121        }
122
123        0.5 + 0.5 * p
124    }
125}
126
127pub fn stdtri(k: isize, p: f64) -> f64 {
128    //! Functional inverse of Student's t distribution
129    //!
130    //! ## DESCRIPTION:
131    //!
132    //! Given probability `p`, finds the argument `t` such that `stdtr(k,t)`
133    //! is equal to `p`.
134    //!
135    //! ## ACCURACY:
136    //!
137    //! Tested at random `1 <= k <= 100`.  The "domain" refers to `p`:
138    //!
139    //! Relative error:
140    //!
141    //!<table>
142    //! <tr>
143    //!     <th>Arithmetic</th>
144    //!     <th>Domain</th>
145    //!     <th># Trials</th>
146    //!     <th>Peak</th>
147    //!     <th>RMS</th>
148    //! </tr>
149    //! <tr>
150    //!     <td>IEEE</td>
151    //!     <td>0.001, 0.999</td>
152    //!     <td>25000</td>
153    //!     <td>5.7e-15</td>
154    //!     <td>8.0e-16</td>
155    //! </tr>
156    //! <tr>
157    //!     <td>IEEE</td>
158    //!     <td>1e-6, 0.001</td>
159    //!     <td>25000</td>
160    //!     <td>2.0e-12</td>
161    //!     <td>2.9e-146</td>
162    //! </tr>
163    //!</table>
164
165    if k <= 0 || !(0.0..=1.0).contains(&p) {
166        //sf_error("stdtri", SF_ERROR_DOMAIN, NULL);
167        f64::NAN
168    } else if p == 0.0 { //M. Romanowicz: p = 0 and 1 should return inf, not nan
169        -f64::INFINITY
170    } else if p == 1.0 {
171        f64::INFINITY
172    } else {
173
174        let rk = k as f64;
175
176        if p > 0.25 && p < 0.75 {
177            if p == 0.5 {
178                return 0.0;
179            }
180            let z = 1.0 - 2.0 * p;
181            let z = incbi(0.5, 0.5 * rk, z.abs());
182            let t = (rk * z / (1.0 - z)).sqrt();
183            if p < 0.5 {
184                -t
185            } else {
186                t
187            }
188        } else {
189            let mut rflg = -1;
190            let mut p = p;
191            if p >= 0.5 {
192                p = 1.0 - p;
193                rflg = 1;
194            }
195            let z = incbi(0.5 * rk, 0.5, 2.0 * p);
196
197            if f64::MAX * z < rk {
198                return rflg as f64 * f64::INFINITY;
199            } 
200            let t = (rk / z - rk).sqrt();
201            rflg as f64 * t
202        }
203    }
204}
205
206#[cfg(test)]
207mod stdtr_tests {
208    use super::*;
209
210    #[test]
211    fn stdtr_trivials() {
212        assert!(stdtr(0, 0.0).is_nan());
213        assert!(stdtr(-1, 0.5).is_nan());
214        assert_eq!(stdtr(1, 0.0), 0.5);
215        assert_eq!(stdtr(1, f64::INFINITY), 1.0);
216        assert_eq!(stdtr(1, -f64::INFINITY), 0.0);
217    }
218
219    #[test]
220    fn stdtr_odd_k() {
221        assert_eq!(stdtr(1, -2.0), 0.14758361765043326);
222        assert_eq!(stdtr(1, -1e-10), 0.499999999968169);
223        assert_eq!(stdtr(1, 1e-10), 0.500000000031831);
224        assert_eq!(stdtr(1, 1.5), 0.8128329581890013);
225        assert_eq!(stdtr(1, 10.0), 0.9682744825694465);
226        assert_eq!(stdtr(1, 1e5), 0.9999968169011383);
227        assert_eq!(stdtr(1, 1e10), 0.999999999968169);
228
229        assert_eq!(stdtr(5, -2.0), 0.05096973941492916);
230        assert_eq!(stdtr(5, -1e-10), 0.4999999999620393);
231        assert_eq!(stdtr(5, 1e-10), 0.5000000000379606);
232        assert_eq!(stdtr(5, 1.5), 0.9030481598787634);
233        assert_eq!(stdtr(5, 10.0), 0.9999145262121285);
234        assert_eq!(stdtr(5, 5e2), 0.9999999999996964);
235
236        assert_eq!(stdtr(101, -2.0), 0.02409260560369547);
237        assert_eq!(stdtr(101, -1e-10), 0.4999999999602044);
238        assert_eq!(stdtr(101, 1e-10), 0.5000000000397956);
239        assert_eq!(stdtr(101, 1.5), 0.9316330369201067);
240        assert_eq!(stdtr(101, 9.0), 0.9999999999999928);
241
242        // TODO: This takes forever and is inaccurate
243        //assert_eq!(stdtr(2147483647, -2.0), 0.022750132011032955);
244    }
245
246    #[test]
247    fn stdtr_even_k() {
248        assert_eq!(stdtr(2, -2.0), 0.09175170953613693);
249        assert_eq!(stdtr(2, -1e-10), 0.49999999996464467);
250        assert_eq!(stdtr(2, 1e-10), 0.5000000000353554);
251        assert_eq!(stdtr(2, 1.5), 0.8638034375544994);
252        assert_eq!(stdtr(2, 10.0), 0.9950737714883372);
253        assert_eq!(stdtr(2, 1e5), 0.99999999995);
254        assert_eq!(stdtr(2, 1e6), 0.9999999999995);
255
256        assert_eq!(stdtr(6, -2.0), 0.046213155765837566);
257        assert_eq!(stdtr(6, -1e-10), 0.4999999999617267);
258        assert_eq!(stdtr(6, 1e-10), 0.5000000000382733);
259        assert_eq!(stdtr(6, 1.5), 0.9078596319292591);
260        assert_eq!(stdtr(6, 10.0), 0.9999710400862252);
261        assert_eq!(stdtr(6, 5e2), 0.9999999999999979);
262
263        assert_eq!(stdtr(100, -2.0), 0.024106089365566685);
264        assert_eq!(stdtr(100, -1e-10), 0.4999999999602054);
265        assert_eq!(stdtr(100, 1e-10), 0.5000000000397946);
266        assert_eq!(stdtr(100, 1.5), 0.9316174709376557);
267        assert_eq!(stdtr(100, 9.0), 0.9999999999999927);
268
269    }
270
271    #[test]
272    fn stdtr_t_lt_neg2() {
273        assert_eq!(stdtr(1, -1e150), 3.1830988618379074e-151);
274        assert_eq!(stdtr(1, -10.0), 0.031725517430553574);
275        assert_eq!(stdtr(1, -2.0 - 1e-10), 0.1475836176440671);
276
277        assert_eq!(stdtr(10, -1e20), 1.2304687499999997e-196);
278        assert_eq!(stdtr(10, -10.0), 7.947765877982062e-7);
279        assert_eq!(stdtr(10, -2.0 - 1e-10), 0.03669401737925559);
280
281        assert_eq!(stdtr(100, -20.0), 4.997133930668492e-37);
282        assert_eq!(stdtr(100, -10.0), 4.95084449229707e-17);
283        assert_eq!(stdtr(100, -2.0 - 1e-10), 0.024106089360076077);
284
285        // TODO: This is inaccurate
286        // assert_eq!(stdtr(2147483647, -20.0), 2.753675665465448e-89);
287        // assert_eq!(stdtr(2147483647, -10.0), 7.61986207143429e-24);
288        // assert_eq!(stdtr(2147483647, -2.0 - 1e-10), 0.022750132005633864);
289        // assert_eq!(stdtr(2147483647, -2.0), 0.022750132011032955);
290    }
291}
292
293#[cfg(test)]
294mod stdtri_tests {
295    use super::*;
296
297    #[test]
298    fn stdtri_trivials() {
299        assert!(stdtri(0, 0.5).is_nan());
300        assert!(stdtri(-1, 0.5).is_nan());
301        assert!(stdtri(1, -1e-10).is_nan());
302        assert!(stdtri(1, 1.0 + 1e10).is_nan());
303        assert_eq!(stdtri(1, 0.0), -f64::INFINITY);
304        assert_eq!(stdtri(1, 1.0), f64::INFINITY);
305        assert_eq!(stdtri(100, 0.5), 0.0);
306    }
307
308    #[test]
309    fn stdtri_p_medium() { // 0.25 < p < 0.75
310
311        // Note: These results do not match scipy's stdtrit, but they appear
312        // to be more accurate
313        assert_eq!(stdtri(1, 0.25 + 1e-10), -0.9999999993716808);
314        assert_eq!(stdtri(1, 0.35), -0.5095254494944288);
315        assert_eq!(stdtri(1, 0.5 - 1e-10), -3.141592913526334e-10);
316        assert_eq!(stdtri(1, 0.5 + 1e-10), 3.141592913526334e-10);
317        assert_eq!(stdtri(1, 0.65), 0.5095254494944288);
318        assert_eq!(stdtri(1, 0.75 - 1e-10), 0.9999999993716808);
319
320        assert_eq!(stdtri(10, 0.25 + 1e-10), -0.6998120609781328);
321        assert_eq!(stdtri(10, 0.35), -0.39659149375562175);
322        assert_eq!(stdtri(10, 0.5 - 1e-10), -2.5699782475714284e-10);
323        assert_eq!(stdtri(10, 0.5 + 1e-10), 2.5699782475714284e-10);
324        assert_eq!(stdtri(10, 0.65), 0.39659149375562175);
325        assert_eq!(stdtri(10, 0.75 - 1e-10), 0.6998120609781328);
326
327        assert_eq!(stdtri(100, 0.25 + 1e-10), -0.6769510426949146);
328        assert_eq!(stdtri(100, 0.35), -0.38642898040767143);
329        assert_eq!(stdtri(100, 0.5 - 1e-10), -2.5129027882894715e-10);
330        assert_eq!(stdtri(100, 0.5 + 1e-10), 2.5129027882894715e-10);
331        assert_eq!(stdtri(100, 0.65), 0.38642898040767143);
332        assert_eq!(stdtri(100, 0.75 - 1e-10), 0.6769510426949146);
333    }
334
335    #[test]
336    fn stdtri_p_large() { // p <= 0.25 or 0.75 <= p
337
338        // Note: These results do not match scipy's stdtrit, but they appear
339        // to be more accurate
340        assert_eq!(stdtri(1, 1e-10), -3183098861.8379073);
341        assert_eq!(stdtri(1, 0.15), -1.9626105055051508);
342        assert_eq!(stdtri(1, 0.25), -1.0000000000000007);
343        assert_eq!(stdtri(1, 0.75), 1.0000000000000007);
344        assert_eq!(stdtri(1, 0.85), 1.9626105055051506);
345        assert_eq!(stdtri(1, 1.0 - 1e-10), 3183098598.467149);
346
347        assert_eq!(stdtri(10, 1e-10), -25.46600802169773);
348        assert_eq!(stdtri(10, 0.15), -1.0930580735905255);
349        assert_eq!(stdtri(10, 0.25), -0.6998120613124323);
350        assert_eq!(stdtri(10, 0.75), 0.6998120613124323);
351        assert_eq!(stdtri(10, 0.85), 1.0930580735905255);
352        assert_eq!(stdtri(10, 1.0 - 1e-10), 25.466007808016013);
353
354        assert_eq!(stdtri(100, 1e-10), -7.083375481400722);
355        assert_eq!(stdtri(100, 0.15), -1.0418359009083464);
356        assert_eq!(stdtri(100, 0.25), -0.6769510430114689);
357        assert_eq!(stdtri(100, 0.75), 0.6769510430114689);
358        assert_eq!(stdtri(100, 0.85), 1.0418359009083464);
359        assert_eq!(stdtri(100, 1.0 - 1e-10), 7.083375464183688);
360    }
361}
spec_math/cephes64/stdtr.rs

spec_math/cephes64/
stdtr.rs
