160 |
* this to be true, the method below can be replaced with that |
* this to be true, the method below can be replaced with that |
161 |
* of the incomplete beta function. |
* of the incomplete beta function. |
162 |
*/ |
*/ |
163 |
int gsl_cdf_t_e ( double t, double nu, gsl_cdf_result * result, |
/* int gsl_cdf_t_e ( double t, double nu) */ |
164 |
gsl_cdf_tail_t tail ) |
/* { */ |
165 |
|
/* int rc; */ |
166 |
|
/* int stopval; */ |
167 |
|
/* double b; */ |
168 |
|
/* double ck; */ |
169 |
|
/* double q; */ |
170 |
|
/* double y; */ |
171 |
|
/* double diff; */ |
172 |
|
/* double ckp2; */ |
173 |
|
/* gsl_sf_result lg1; */ |
174 |
|
/* gsl_sf_result lg2; */ |
175 |
|
/* double num; */ |
176 |
|
/* int k; */ |
177 |
|
/* double u; */ |
178 |
|
/* double p; */ |
179 |
|
/* int i; */ |
180 |
|
/* int inu = (int) nu; */ |
181 |
|
|
182 |
|
/* /\* */ |
183 |
|
/* * First approximation is also found in Abramowitz */ |
184 |
|
/* * and Stegun. This is used only for small values of t */ |
185 |
|
/* * and small degrees of freedom. */ |
186 |
|
/* *\/ */ |
187 |
|
/* y = t / sqrt(nu); */ |
188 |
|
/* if ( fabs(t) < 4.0 ) */ |
189 |
|
/* { */ |
190 |
|
/* diff = fabs(nu-1.0); */ |
191 |
|
/* if ( diff < GSL_DBL_EPSILON ) */ |
192 |
|
/* { */ |
193 |
|
/* q = M_PI_2 * atan(y); */ |
194 |
|
/* } */ |
195 |
|
/* else if ( (nu < 21.0) && (nu > 1.0)) */ |
196 |
|
/* { */ |
197 |
|
/* ckp2 = 1.0; */ |
198 |
|
/* b = 1.0 + t*t/nu; */ |
199 |
|
|
200 |
|
/* for ( k = inu-2; k > 1; k-=2) */ |
201 |
|
/* { */ |
202 |
|
/* ck = 1 + (ckp2 * (k-1)) / ((double)k*b); */ |
203 |
|
/* ckp2 = ck; */ |
204 |
|
/* } */ |
205 |
|
/* if ( (inu%2) == 0 ) */ |
206 |
|
/* { */ |
207 |
|
/* q = ck * y / sqrt(b); */ |
208 |
|
/* } */ |
209 |
|
/* else */ |
210 |
|
/* { */ |
211 |
|
/* q = M_PI_2 * (atan(y) + ck * y / b); */ |
212 |
|
/* } */ |
213 |
|
/* if ( ((tail == GSL_CDF_UPPER ) && (t > 0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */ |
214 |
|
/* { */ |
215 |
|
/* result->val = (1.0 - q)/2; */ |
216 |
|
/* } */ |
217 |
|
/* else */ |
218 |
|
/* { */ |
219 |
|
/* result->val = (1.0+q)/2; */ |
220 |
|
/* } */ |
221 |
|
/* } */ |
222 |
|
/* else if ( nu >= 21.0 ) */ |
223 |
|
/* { */ |
224 |
|
/* |
225 |
|
* Use the Cornish-Fisher expansion to find |
226 |
|
* a point u such that gsl_cdf_gauss(u) = tcdf(t). |
227 |
|
* Approximate the t cdf with gsl_cdf_gauss. |
228 |
|
*/ |
229 |
|
/* u = cornish_fisher ( t, nu ); */ |
230 |
|
/* q = gsl_cdf_gauss_Q ( u ); */ |
231 |
|
/* if ( ((tail == GSL_CDF_UPPER )&&(t>0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */ |
232 |
|
/* { */ |
233 |
|
/* result->val = q; */ |
234 |
|
/* } */ |
235 |
|
/* else */ |
236 |
|
/* { */ |
237 |
|
/* result->val = 1.0 - q; */ |
238 |
|
/* } */ |
239 |
|
/* } */ |
240 |
|
/* else */ |
241 |
|
/* { */ |
242 |
|
/* GSL_ERROR_VAL ("degrees of freedom are less than 1", */ |
243 |
|
/* GSL_EINVAL, nu ); */ |
244 |
|
/* } */ |
245 |
|
/* } */ |
246 |
|
/* return rc; */ |
247 |
|
/* } */ |
248 |
|
|
249 |
|
/* |
250 |
|
* First approximation is also found in Abramowitz |
251 |
|
* and Stegun. This is used only for small values of t |
252 |
|
* and small degrees of freedom. |
253 |
|
*/ |
254 |
|
static double gsl_cdf_t_smallx ( const double t, const double nu ) |
255 |
{ |
{ |
256 |
int rc; |
double y; |
257 |
int stopval; |
double q; |
258 |
|
int i; |
259 |
double b; |
double b; |
260 |
|
double diff; |
261 |
|
double ckp2; |
262 |
double ck; |
double ck; |
263 |
|
int k; |
264 |
|
int inu = (int) nu; |
265 |
|
|
266 |
|
y = t / sqrt(nu); |
267 |
|
diff = fabs(nu-1.0); |
268 |
|
if ( diff < GSL_DBL_EPSILON ) |
269 |
|
{ |
270 |
|
q = M_PI_2 * atan(y); |
271 |
|
} |
272 |
|
else if ( (nu < 21.0) && (nu > 1.0)) |
273 |
|
{ |
274 |
|
ckp2 = 1.0; |
275 |
|
b = 1.0 + t*t/nu; |
276 |
|
|
277 |
|
for ( k = inu-2; k > 1; k-=2) |
278 |
|
{ |
279 |
|
ck = 1.0 + (ckp2 * (k-1)) / ((double)k*b); |
280 |
|
ckp2 = ck; |
281 |
|
} |
282 |
|
if ( (inu%2) == 0 ) |
283 |
|
{ |
284 |
|
q = ck * y / sqrt(b); |
285 |
|
} |
286 |
|
else |
287 |
|
{ |
288 |
|
q = M_PI_2 * (atan(y) + ck * y / b); |
289 |
|
} |
290 |
|
/* if ( ((tail == GSL_CDF_UPPER ) && (t > 0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */ |
291 |
|
/* { */ |
292 |
|
/* result->val = (1.0 - q)/2; */ |
293 |
|
/* } */ |
294 |
|
/* else */ |
295 |
|
/* { */ |
296 |
|
/* result->val = (1.0+q)/2; */ |
297 |
|
/* } */ |
298 |
|
} |
299 |
|
return q; |
300 |
|
} |
301 |
|
/* |
302 |
|
* Use the Cornish-Fisher expansion to find |
303 |
|
* a point u such that gsl_cdf_gauss(u) = tcdf(t). |
304 |
|
* Approximate the t cdf with gsl_cdf_gauss. |
305 |
|
*/ |
306 |
|
static double t_cornish_fisher ( const double x, const double nu ) |
307 |
|
{ |
308 |
|
double y; |
309 |
double q; |
double q; |
310 |
|
double u; |
311 |
|
|
312 |
|
y = x / sqrt(nu); |
313 |
|
u = cornish_fisher ( x, nu ); |
314 |
|
q = gsl_cdf_gauss_Q ( u); |
315 |
|
/* if ( ((tail == GSL_CDF_UPPER )&&(t>0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */ |
316 |
|
/* { */ |
317 |
|
/* result->val = q; */ |
318 |
|
/* } */ |
319 |
|
/* else */ |
320 |
|
/* { */ |
321 |
|
/* result->val = 1.0 - q; */ |
322 |
|
/* } */ |
323 |
|
return q; |
324 |
|
} |
325 |
|
/* |
326 |
|
* Series approximation for t > 4.0. This needs to be fixed; |
327 |
|
* it shouldn't subtract the result from 1.0. A better way is |
328 |
|
* to use two different series expansions. Figuring this out |
329 |
|
* means rummaging through Fisher's paper in Metron, v5, 1926, |
330 |
|
* "Expansion of Student's integral in powers of n^{-1}." |
331 |
|
*/ |
332 |
|
static double normal_approx ( const double x, const double nu ) |
333 |
|
{ |
334 |
double y; |
double y; |
335 |
|
double num; |
336 |
double diff; |
double diff; |
337 |
double ckp2; |
double q; |
338 |
|
int i; |
339 |
|
int rc; |
340 |
|
|
341 |
gsl_sf_result lg1; |
gsl_sf_result lg1; |
342 |
gsl_sf_result lg2; |
gsl_sf_result lg2; |
|
double num; |
|
|
int k; |
|
|
double u; |
|
|
double p; |
|
|
int i; |
|
|
int inu = (int) nu; |
|
343 |
|
|
344 |
/* |
y = 1/sqrt(1+x*x/nu); |
345 |
* First approximation is also found in Abramowitz |
num = 1.0; |
346 |
* and Stegun. This is used only for small values of t |
q = 0.0; |
347 |
* and small degrees of freedom. |
diff = 2*GSL_DBL_EPSILON; |
348 |
*/ |
for ( i = 2; (i < MAXI)&&(diff>GSL_DBL_EPSILON); i+=2) |
349 |
y = t / sqrt(nu); |
{ |
350 |
if ( fabs(t) < 4.0 ) |
diff = q; |
351 |
|
num *= y*y*(i-1)/i; |
352 |
|
q += num / (nu+i); |
353 |
|
diff = q - diff; |
354 |
|
} |
355 |
|
q += 1/nu; |
356 |
|
rc = gsl_sf_lngamma_e ((nu / 2.0), &lg1); |
357 |
|
if ( rc == GSL_SUCCESS ) |
358 |
{ |
{ |
359 |
diff = fabs(nu-1.0); |
rc = gsl_sf_lngamma_e (((nu + 1) / 2.0), &lg2); |
360 |
if ( diff < GSL_DBL_EPSILON ) |
if ( rc != GSL_SUCCESS ) |
361 |
{ |
{ |
362 |
q = M_PI_2 * atan(y); |
return rc; |
363 |
} |
} |
364 |
else if ( (nu < 21.0) && (nu > 1.0)) |
} |
365 |
|
diff = (lg2.val) - (lg1.val); |
366 |
|
q *= pow(y,nu) * exp (diff) / sqrt (M_PI); |
367 |
|
|
368 |
|
/* if ( ((tail == GSL_CDF_UPPER )&&(t>0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) */ |
369 |
|
/* { */ |
370 |
|
/* result->val = q; */ |
371 |
|
/* } */ |
372 |
|
/* else if ( tail == GSL_CDF_LOWER ) */ |
373 |
|
/* { */ |
374 |
|
/* result->val = 1.0 - q; */ |
375 |
|
/* } */ |
376 |
|
|
377 |
|
return q; |
378 |
|
} |
379 |
|
|
380 |
|
double gsl_cdf_t_P ( const double x, const double nu ) |
381 |
|
{ |
382 |
|
double val = 0.0; |
383 |
|
double q; |
384 |
|
double absx = fabs(x); |
385 |
|
|
386 |
|
if ( fabs(x) < 4.0 ) |
387 |
|
{ |
388 |
|
q = gsl_cdf_t_smallx (x, nu ); |
389 |
|
if ( x <= 0.0 ) |
390 |
{ |
{ |
391 |
ckp2 = 1.0; |
val = (1.0 - q) / 2.0 ; |
|
b = 1.0 + t*t/nu; |
|
|
|
|
|
for ( k = inu-2; k > 1; k-=2) |
|
|
{ |
|
|
ck = 1 + (ckp2 * (k-1)) / ((double)k*b); |
|
|
ckp2 = ck; |
|
|
} |
|
|
if ( (inu%2) == 0 ) |
|
|
{ |
|
|
q = ck * y / sqrt(b); |
|
|
} |
|
|
else |
|
|
{ |
|
|
q = M_PI_2 * (atan(y) + ck * y / b); |
|
|
} |
|
|
if ( ((tail == GSL_CDF_UPPER ) && (t > 0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) |
|
|
{ |
|
|
result->val = (1.0 - q)/2; |
|
|
} |
|
|
else |
|
|
{ |
|
|
result->val = (1.0+q)/2; |
|
|
} |
|
392 |
} |
} |
393 |
else if ( nu >= 21.0 ) |
else |
394 |
{ |
{ |
395 |
/* |
val = (1.0 + q ) / 2.0; |
396 |
* Use the Cornish-Fisher expansion to find |
} |
397 |
* a point u such that gsl_cdf_gauss(u) = tcdf(t). |
} |
398 |
* Approximate the t cdf with gsl_cdf_gauss. |
else if ( nu >= 21.0 ) |
399 |
*/ |
{ |
400 |
u = cornish_fisher ( t, nu ); |
q = t_cornish_fisher ( x, nu ); |
401 |
q = gsl_cdf_gauss ( u, GSL_CDF_UPPER ); |
if ( x <= 0.0) |
402 |
if ( ((tail == GSL_CDF_UPPER )&&(t>0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) |
{ |
403 |
{ |
val = q; |
|
result->val = q; |
|
|
} |
|
|
else |
|
|
{ |
|
|
result->val = 1.0 - q; |
|
|
} |
|
404 |
} |
} |
405 |
else |
else |
406 |
{ |
{ |
407 |
GSL_ERROR_VAL ("degrees of freedom are less than 1", |
val = 1.0 - q; |
|
GSL_EINVAL, nu ); |
|
408 |
} |
} |
409 |
} |
} |
410 |
else |
else |
411 |
{ |
{ |
412 |
/* |
val = normal_approx ( x, nu ); |
413 |
* Series approximation for t > 4.0. This needs to be fixed; |
if ( x > 0.0 ) |
414 |
* it shouldn't subtract the result from 1.0. A better way is |
{ |
415 |
* to use two different series expansions. Figuring this out |
val = 1.0 - val; |
416 |
* means rummaging through Fisher's paper in Metron, v5, 1926, |
} |
417 |
* "Expansion of Student's integral in powers of n^{-1}." |
} |
418 |
*/ |
return val; |
419 |
|
} |
420 |
|
double gsl_cdf_t_Q ( const double x, const double nu ) |
421 |
y = 1/sqrt(1+t*t/nu); |
{ |
422 |
num = 1.0; |
double val = 0.0; |
423 |
p = 0.0; |
double q; |
424 |
diff = 2*GSL_DBL_EPSILON; |
double absx = fabs(x); |
425 |
for ( i = 2; (i < MAXI)&&(diff>GSL_DBL_EPSILON); i+=2) |
|
426 |
{ |
if ( fabs(x) < 4.0 ) |
427 |
diff = p; |
{ |
428 |
num *= y*y*(i-1)/i; |
q = gsl_cdf_t_smallx ( x, nu ); |
429 |
p += num / (nu+i); |
if ( x > 0.0 ) |
430 |
diff = p - diff; |
{ |
431 |
} |
val = (1.0 - q) / 2.0 ; |
|
p += 1/nu; |
|
|
rc = gsl_sf_lngamma_e ((nu / 2.0), &lg1); |
|
|
if ( rc == GSL_SUCCESS ) |
|
|
{ |
|
|
rc = gsl_sf_lngamma_e (((nu + 1) / 2.0), &lg2); |
|
|
if ( rc != GSL_SUCCESS ) |
|
|
{ |
|
|
return rc; |
|
|
} |
|
432 |
} |
} |
433 |
else |
else |
434 |
{ |
{ |
435 |
return rc; |
val = (1.0 + q ) / 2.0; |
436 |
} |
} |
437 |
diff = (lg2.val) - (lg1.val); |
} |
438 |
p *= pow(y,nu) * exp (diff) / sqrt (M_PI); |
else if ( nu >= 21.0 ) |
439 |
|
{ |
440 |
if ( ((tail == GSL_CDF_UPPER )&&(t>0.0))||((tail==GSL_CDF_LOWER)&&(t<=0.0))) |
val = t_cornish_fisher ( x, nu ); |
441 |
|
if ( x <= 0.0) |
442 |
{ |
{ |
443 |
result->val = p; |
val = 1.0 - val; |
444 |
} |
} |
445 |
else if ( tail == GSL_CDF_LOWER ) |
} |
446 |
|
else |
447 |
|
{ |
448 |
|
val = normal_approx ( x, nu ); |
449 |
|
if ( x <= 0.0 ) |
450 |
{ |
{ |
451 |
result->val = 1.0 - p; |
val = 1.0 - val; |
452 |
} |
} |
453 |
} |
} |
454 |
return rc; |
return val; |
455 |
} |
} |