1 /* stdtr.c
2 *
3 * Student's t distribution
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double t, stdtr();
10 * short k;
11 *
12 * y = stdtr( k, t );
13 *
14 *
15 * DESCRIPTION:
16 *
17 * Computes the integral from minus infinity to t of the Student
18 * t distribution with integer k > 0 degrees of freedom:
19 *
20 * t
21 * -
22 * | |
23 * - | 2 -(k+1)/2
24 * | ( (k+1)/2 ) | ( x )
25 * ---------------------- | ( 1 + --- ) dx
26 * - | ( k )
27 * sqrt( k pi ) | ( k/2 ) |
28 * | |
29 * -
30 * -inf.
31 *
32 * Relation to incomplete beta integral:
33 *
34 * 1 - stdtr(k,t) = 0.5 * incbet( k/2, 1/2, z )
35 * where
36 * z = k/(k + t**2).
37 *
38 * For t < -2, this is the method of computation. For higher t,
39 * a direct method is derived from integration by parts.
40 * Since the function is symmetric about t=0, the area under the
41 * right tail of the density is found by calling the function
42 * with -t instead of t.
43 *
44 * ACCURACY:
45 *
46 * Tested at random 1 <= k <= 25. The "domain" refers to t.
47 * Relative error:
48 * arithmetic domain # trials peak rms
49 * IEEE -100,-2 50000 5.9e-15 1.4e-15
50 * IEEE -2,100 500000 2.7e-15 4.9e-17
51 */
52
53 /* stdtri.c
54 *
55 * Functional inverse of Student's t distribution
56 *
57 *
58 *
59 * SYNOPSIS:
60 *
61 * double p, t, stdtri();
62 * int k;
63 *
64 * t = stdtri( k, p );
65 *
66 *
67 * DESCRIPTION:
68 *
69 * Given probability p, finds the argument t such that stdtr(k,t)
70 * is equal to p.
71 *
72 * ACCURACY:
73 *
74 * Tested at random 1 <= k <= 100. The "domain" refers to p:
75 * Relative error:
76 * arithmetic domain # trials peak rms
77 * IEEE .001,.999 25000 5.7e-15 8.0e-16
78 * IEEE 10^-6,.001 25000 2.0e-12 2.9e-14
79 */
80
81
82 /*
83 * Cephes Math Library Release 2.3: March, 1995
84 * Copyright 1984, 1987, 1995 by Stephen L. Moshier
85 */
86
87 #include "mconf.h"
88 #include <float.h>
89
90 extern double MACHEP;
91
stdtr(k,t)92 double stdtr(k, t)
93 int k;
94 double t;
95 {
96 double x, rk, z, f, tz, p, xsqk;
97 int j;
98
99 if (k <= 0) {
100 sf_error("stdtr", SF_ERROR_DOMAIN, NULL);
101 return (NPY_NAN);
102 }
103
104 if (t == 0)
105 return (0.5);
106
107 if (t < -2.0) {
108 rk = k;
109 z = rk / (rk + t * t);
110 p = 0.5 * incbet(0.5 * rk, 0.5, z);
111 return (p);
112 }
113
114 /* compute integral from -t to + t */
115
116 if (t < 0)
117 x = -t;
118 else
119 x = t;
120
121 rk = k; /* degrees of freedom */
122 z = 1.0 + (x * x) / rk;
123
124 /* test if k is odd or even */
125 if ((k & 1) != 0) {
126
127 /* computation for odd k */
128
129 xsqk = x / sqrt(rk);
130 p = atan(xsqk);
131 if (k > 1) {
132 f = 1.0;
133 tz = 1.0;
134 j = 3;
135 while ((j <= (k - 2)) && ((tz / f) > MACHEP)) {
136 tz *= (j - 1) / (z * j);
137 f += tz;
138 j += 2;
139 }
140 p += f * xsqk / z;
141 }
142 p *= 2.0 / NPY_PI;
143 }
144
145
146 else {
147
148 /* computation for even k */
149
150 f = 1.0;
151 tz = 1.0;
152 j = 2;
153
154 while ((j <= (k - 2)) && ((tz / f) > MACHEP)) {
155 tz *= (j - 1) / (z * j);
156 f += tz;
157 j += 2;
158 }
159 p = f * x / sqrt(z * rk);
160 }
161
162 /* common exit */
163
164
165 if (t < 0)
166 p = -p; /* note destruction of relative accuracy */
167
168 p = 0.5 + 0.5 * p;
169 return (p);
170 }
171
stdtri(k,p)172 double stdtri(k, p)
173 int k;
174 double p;
175 {
176 double t, rk, z;
177 int rflg;
178
179 if (k <= 0 || p <= 0.0 || p >= 1.0) {
180 sf_error("stdtri", SF_ERROR_DOMAIN, NULL);
181 return (NPY_NAN);
182 }
183
184 rk = k;
185
186 if (p > 0.25 && p < 0.75) {
187 if (p == 0.5)
188 return (0.0);
189 z = 1.0 - 2.0 * p;
190 z = incbi(0.5, 0.5 * rk, fabs(z));
191 t = sqrt(rk * z / (1.0 - z));
192 if (p < 0.5)
193 t = -t;
194 return (t);
195 }
196 rflg = -1;
197 if (p >= 0.5) {
198 p = 1.0 - p;
199 rflg = 1;
200 }
201 z = incbi(0.5 * rk, 0.5, 2.0 * p);
202
203 if (DBL_MAX * z < rk)
204 return (rflg * NPY_INFINITY);
205 t = sqrt(rk / z - rk);
206 return (rflg * t);
207 }
208