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.8: June, 2000
84 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
85 */
86
87 #include "mconf.h"
88
stdtr(double rk,double t)89 double stdtr (double rk, double t)
90 {
91 double x, z, f, tz, p, xsqk;
92 int k, j, is_int;
93
94 if (rk <= 0) {
95 mtherr("stdtr", CEPHES_DOMAIN);
96 return 0.0;
97 } else {
98 k = (int) rk;
99 is_int = (rk - k) == 0.0;
100 }
101
102 if (t == 0) {
103 return 0.5;
104 }
105
106 if (t < -2.0) {
107 z = rk / (rk + t * t);
108 p = 0.5 * incbet(0.5*rk, 0.5, z);
109 return p;
110 }
111
112 if (is_int) {
113 /* integer df */
114
115 /* compute integral from -t to +t */
116
117 if (t < 0) {
118 x = -t;
119 } else {
120 x = t;
121 }
122
123 z = 1.0 + (x * x)/rk;
124
125 /* test if k is odd or even */
126 if ((k & 1) != 0) {
127
128 /* computation for odd k */
129
130 xsqk = x/sqrt(rk);
131 p = atan(xsqk);
132 if (k > 1) {
133 f = 1.0;
134 tz = 1.0;
135 j = 3;
136 while ((j <= (k-2)) && ((tz/f) > MACHEP)) {
137 tz *= (j - 1)/(z * j);
138 f += tz;
139 j += 2;
140 }
141 p += f * xsqk/z;
142 }
143 p *= 2.0/PI;
144 } else {
145
146 /* computation for even k */
147
148 f = 1.0;
149 tz = 1.0;
150 j = 2;
151
152 while ((j <= (k-2)) && ((tz/f) > MACHEP)) {
153 tz *= (j - 1)/(z * j);
154 f += tz;
155 j += 2;
156 }
157 p = f * x/sqrt(z*rk);
158 }
159
160 /* common exit */
161
162 if (t < 0) {
163 p = -p; /* note destruction of relative accuracy */
164 }
165
166 p = 0.5 + 0.5 * p;
167
168 } else {
169 /* non-integer df */
170 z = rk / (rk + t * t);
171 p = 0.5 * incbet(0.5*rk, 0.5, z);
172 if (t > 0) {
173 p = 1 - p;
174 }
175 }
176
177
178 return p;
179 }
180
stdtri(double rk,double p)181 double stdtri (double rk, double p)
182 {
183 double t, z;
184 int rflg;
185
186 if (rk <= 0 || p <= 0.0 || p >= 1.0) {
187 mtherr("stdtri", CEPHES_DOMAIN);
188 return 0.0;
189 }
190
191 if (p > 0.25 && p < 0.75) {
192 if (p == 0.5) {
193 return 0.0;
194 }
195 z = 1.0 - 2.0 * p;
196 z = incbi(0.5, 0.5*rk, fabs(z));
197 t = sqrt(rk * z/(1.0-z));
198 if (p < 0.5) {
199 t = -t;
200 }
201 return t;
202 }
203
204 rflg = -1;
205
206 if (p >= 0.5) {
207 p = 1.0 - p;
208 rflg = 1;
209 }
210
211 z = incbi(0.5*rk, 0.5, 2.0*p);
212
213 if (MAXNUM * z < rk) {
214 return rflg * MAXNUM;
215 }
216
217 t = sqrt(rk/z - rk);
218
219 return rflg * t;
220 }
221