451 |
* Invert the T distribution. Uses a method as shown in |
* Invert the T distribution. Uses a method as shown in |
452 |
* Statistical Computing, 5.4.2. This method uses an initial |
* Statistical Computing, 5.4.2. This method uses an initial |
453 |
* approximation based on Hill's method above, then improves |
* approximation based on Hill's method above, then improves |
454 |
* the initial method with a Taylor series. |
* the initial method with a Taylor series or a Cornish-Fisher |
455 |
|
* expansion. |
456 |
*/ |
*/ |
457 |
|
static double inv_cornish_fisher ( double z, double n ) |
458 |
|
{ |
459 |
|
double ret_val = 0.0; |
460 |
|
double b; |
461 |
|
double zsq; |
462 |
|
double tmp; |
463 |
|
double c; |
464 |
|
double u; |
465 |
|
double d; |
466 |
|
double zcube; |
467 |
|
double zfour; |
468 |
|
double zfive; |
469 |
|
double zseven; |
470 |
|
int i; |
471 |
|
|
472 |
|
tmp = n-0.5; |
473 |
|
b = 48.0 * tmp * tmp; |
474 |
|
if ( n > 5.0 ) |
475 |
|
{ |
476 |
|
c = 96.36 - 16.0 / tmp - 98.0 / (tmp*tmp) + 20700.0 / (tmp*tmp*tmp*b); |
477 |
|
} |
478 |
|
else |
479 |
|
{ |
480 |
|
c = 0.3 * (n - 4.5) * (z + 0.6); |
481 |
|
} |
482 |
|
zsq = z*z; |
483 |
|
zcube = zsq*z; |
484 |
|
zfour = zcube*z; |
485 |
|
zfive = zfour*z; |
486 |
|
zseven = zfour*zcube; |
487 |
|
d = n * M_SQRTPI * gsl_sf_gamma ( tmp ) |
488 |
|
/ (2.0*gsl_sf_gamma ( tmp + 0.5 )); |
489 |
|
u = 10.0 * b * ( b + c - 2.0*z - 7.0 * zsq - 5.0*zcube + 0.05*d*zfour); |
490 |
|
tmp = 4.0*zseven + 63.0*zfive + 360.0*zcube + 945.0*z; |
491 |
|
ret_val = z - (zcube+3.0*z)/b + tmp/u; |
492 |
|
return ret_val; |
493 |
|
} |
494 |
|
|
495 |
double gsl_cdf_ut_P_inv (double prob, double nu) |
double gsl_cdf_ut_P_inv (double prob, double nu) |
496 |
{ |
{ |
497 |
double initial_result; |
double initial_result; |
548 |
*/ |
*/ |
549 |
if ( method_test > 0.05 ) |
if ( method_test > 0.05 ) |
550 |
{ |
{ |
551 |
tmp = prob / 2.0; |
tmp = prob; |
552 |
a = nu - 0.5; |
a = nu - 0.5; |
553 |
/* gsl_cdf_ugaussian_P_inv(tmp) ? */ |
/* gsl_cdf_ugaussian_P_inv(tmp) ? */ |
554 |
x = gsl_cdf_ugaussian_Q_inv (tmp ); |
x = gsl_cdf_ugaussian_P_inv (tmp ); |
555 |
y = cornish_fisher ( x, nu ); |
y = inv_cornish_fisher ( x, nu ); |
556 |
tsqn = -1.0 + exp ( a*y*y); |
tsqn = -1.0 + exp ( a*y*y); |
557 |
} |
} |
558 |
else |
else |
598 |
printf("result = %f\n",result); |
printf("result = %f\n",result); |
599 |
return result; |
return result; |
600 |
} |
} |
|
|
|