447 |
} |
} |
448 |
return val; |
return val; |
449 |
} |
} |
450 |
|
/* |
451 |
|
* Invert the T distribution. Uses a method as shown in |
452 |
|
* Statistical Computing, 5.4.2. This method uses an initial |
453 |
|
* approximation based on Hill's method above, then improves |
454 |
|
* the initial method with a Taylor series. |
455 |
|
*/ |
456 |
|
double gsl_cdf_ut_P_inv (double prob, double nu) |
457 |
|
{ |
458 |
|
double initial_result; |
459 |
|
double result; |
460 |
|
double d; |
461 |
|
double tmp; |
462 |
|
double tmp2; |
463 |
|
double z; |
464 |
|
double zz; |
465 |
|
double zzz; |
466 |
|
double y; |
467 |
|
double tsqn; |
468 |
|
double a; |
469 |
|
double method_test; |
470 |
|
double psi; |
471 |
|
double psi_prime; |
472 |
|
double d_psi_dt; |
473 |
|
double d_psiprime_dt; |
474 |
|
double tcdf; |
475 |
|
double x; |
476 |
|
double w; |
477 |
|
double c2; |
478 |
|
double c3; |
479 |
|
double c4; |
480 |
|
|
481 |
|
if ( prob < 0.0 ) |
482 |
|
{ |
483 |
|
return GSL_EDOM; |
484 |
|
} |
485 |
|
if ( prob > 1.0 ) |
486 |
|
{ |
487 |
|
return GSL_EDOM; |
488 |
|
} |
489 |
|
if ( nu < 0.0 ) |
490 |
|
{ |
491 |
|
return GSL_EDOM; |
492 |
|
} |
493 |
|
if ( fabs(prob) < GSL_DBL_EPSILON ) |
494 |
|
{ |
495 |
|
return GSL_POSINF; |
496 |
|
} |
497 |
|
if ( fabs ( 1.0 - prob ) < GSL_DBL_EPSILON ) |
498 |
|
{ |
499 |
|
return GSL_NEGINF; |
500 |
|
} |
501 |
|
printf("prob is %f\t",prob); |
502 |
|
tmp = nu / 2.0; |
503 |
|
d = tmp * M_SQRTPI * gsl_sf_gamma ( tmp ) |
504 |
|
/ gsl_sf_gamma ( tmp + 0.5 ); |
505 |
|
method_test = gsl_max ( d*prob, 2.0 / nu); |
506 |
|
/* |
507 |
|
* There are two possible initial approximations. |
508 |
|
* Which is used depends on prob and nu. |
509 |
|
*/ |
510 |
|
if ( method_test > 0.05 ) |
511 |
|
{ |
512 |
|
tmp = prob / 2.0; |
513 |
|
a = nu - 0.5; |
514 |
|
/* gsl_cdf_ugaussian_P_inv(tmp) ? */ |
515 |
|
x = gsl_cdf_ugaussian_Q_inv (tmp ); |
516 |
|
y = cornish_fisher ( x, nu ); |
517 |
|
tsqn = -1.0 + exp ( a*y*y); |
518 |
|
} |
519 |
|
else |
520 |
|
{ |
521 |
|
z = pow ( prob * d, 1 / tmp ); |
522 |
|
zz = z * z; |
523 |
|
zzz = z*z*z; |
524 |
|
tsqn = 1/z + (nu + 1.0)*(-1.0 + z/(2.0 * (nu + 4.0)) + |
525 |
|
nu * zz / (3.0 * (nu+2.0)*(nu+6.0)) |
526 |
|
+ nu*(nu+3.0)*(2.0*nu*nu+9.0*nu-2.0)*zzz |
527 |
|
/(8.0 * (nu + 2.0) * (nu + 2.0) |
528 |
|
* (nu + 4.0) * (nu + 4.0) * |
529 |
|
(nu + 8.0))) |
530 |
|
/(nu + 2); |
531 |
|
|
532 |
|
} |
533 |
|
initial_result = sqrt(tsqn * nu); |
534 |
|
printf("initial_result = %f \t",initial_result); |
535 |
|
tcdf = gsl_cdf_t_Q( initial_result, nu ); |
536 |
|
tmp = 0.5 * ( tcdf - prob ); |
537 |
|
tmp2 = tdist_pdf ( initial_result, nu ); |
538 |
|
w = tmp / tmp2; |
539 |
|
psi = initial_result * (nu+1.0) / |
540 |
|
(nu + initial_result * initial_result); |
541 |
|
psi_prime = (nu+1.0) * (nu - initial_result * initial_result) / |
542 |
|
((nu + initial_result*initial_result) * |
543 |
|
(nu + initial_result * initial_result)); |
544 |
|
d_psi_dt = 2 * initial_result * nu * (nu + 1.0) / |
545 |
|
((nu + initial_result * initial_result) * |
546 |
|
(nu + initial_result * initial_result)); |
547 |
|
d_psiprime_dt = -2.0 * initial_result * (nu + 1.0) * |
548 |
|
(3.0 * nu - initial_result * initial_result) / |
549 |
|
((nu + initial_result * initial_result) * |
550 |
|
(nu + initial_result * initial_result) * |
551 |
|
(nu + initial_result * initial_result)); |
552 |
|
c2 = psi/2.0; |
553 |
|
c3 = (2.0 * psi * psi + psi_prime) / 6.0; |
554 |
|
c4 = 3.0 * psi * (2.0 * psi * psi + psi_prime) + |
555 |
|
4.0 * psi * d_psi_dt + d_psiprime_dt; |
556 |
|
c4 /= 24.0; |
557 |
|
result = initial_result + w + c2 * w * w + |
558 |
|
c3 * w * w * w + c4 * w * w * w * w; |
559 |
|
printf("result = %f\n",result); |
560 |
|
return result; |
561 |
|
} |
562 |
|
|