1 |
/* cdf/t.c |
/* cdf/tdist.c |
2 |
* |
* |
3 |
* Copyright (C) 2002 Jason H. Stover. |
* Copyright (C) 2002 Jason H. Stover. |
4 |
* |
* |
263 |
|
|
264 |
|
|
265 |
|
|
|
#if 0 |
|
|
/* |
|
|
* Invert the T distribution. Uses a method as shown in |
|
|
* Statistical Computing, 5.4.2. This method uses an initial |
|
|
* approximation based on Hill's method above, then improves |
|
|
* the initial method with a Taylor series or a Cornish-Fisher |
|
|
* expansion. |
|
|
*/ |
|
|
static double |
|
|
inv_cornish_fisher (double z, double n) |
|
|
{ |
|
|
double ret_val = 0.0; |
|
|
double b; |
|
|
double zsq; |
|
|
double tmp; |
|
|
double c; |
|
|
double u; |
|
|
double d; |
|
|
double zcube; |
|
|
double zfour; |
|
|
double zfive; |
|
|
double zseven; |
|
|
int i; |
|
|
|
|
|
tmp = n - 0.5; |
|
|
b = 48.0 * tmp * tmp; |
|
|
if (n > 5.0) |
|
|
{ |
|
|
c = |
|
|
96.36 - 16.0 / tmp - 98.0 / (tmp * tmp) + |
|
|
20700.0 / (tmp * tmp * tmp * b); |
|
|
} |
|
|
else |
|
|
{ |
|
|
c = 0.3 * (n - 4.5) * (z + 0.6); |
|
|
} |
|
|
zsq = z * z; |
|
|
zcube = zsq * z; |
|
|
zfour = zcube * z; |
|
|
zfive = zfour * z; |
|
|
zseven = zfour * zcube; |
|
|
d = n * M_SQRTPI * gsl_sf_gamma (tmp) / (2.0 * gsl_sf_gamma (tmp + 0.5)); |
|
|
u = |
|
|
10.0 * b * (b + c - 2.0 * z - 7.0 * zsq - 5.0 * zcube + 0.05 * d * zfour); |
|
|
tmp = 4.0 * zseven + 63.0 * zfive + 360.0 * zcube + 945.0 * z; |
|
|
ret_val = z - (zcube + 3.0 * z) / b + tmp / u; |
|
|
return ret_val; |
|
|
} |
|
|
|
|
|
static double |
|
|
tdist_pdf (const double x, const double nu) |
|
|
{ |
|
|
double p; |
|
|
|
|
|
double lg1 = gsl_sf_lngamma (nu / 2); |
|
|
double lg2 = gsl_sf_lngamma ((nu + 1) / 2); |
|
|
|
|
|
p = ((exp (lg2 - lg1) / sqrt (M_PI * nu)) |
|
|
* pow ((1 + x * x / nu), -(nu + 1) / 2)); |
|
|
return p; |
|
|
} |
|
|
|
|
|
double |
|
|
gsl_cdf_ut_Pinv (double prob, double nu) |
|
|
{ |
|
|
double initial_result; |
|
|
double result; |
|
|
double d; |
|
|
double tmp; |
|
|
double tmp2; |
|
|
double z; |
|
|
double zz; |
|
|
double zzz; |
|
|
double y; |
|
|
double tsqn; |
|
|
double a; |
|
|
double method_test; |
|
|
double psi; |
|
|
double psi_prime; |
|
|
double d_psi_dt; |
|
|
double d_psiprime_dt; |
|
|
double tcdf; |
|
|
double x; |
|
|
double w; |
|
|
double c2; |
|
|
double c3; |
|
|
double c4; |
|
|
|
|
|
if (prob < 0.0) |
|
|
{ |
|
|
return GSL_EDOM; |
|
|
} |
|
|
if (prob > 1.0) |
|
|
{ |
|
|
return GSL_EDOM; |
|
|
} |
|
|
if (nu < 0.0) |
|
|
{ |
|
|
return GSL_EDOM; |
|
|
} |
|
|
if (fabs (prob) < GSL_DBL_EPSILON) |
|
|
{ |
|
|
return GSL_POSINF; |
|
|
} |
|
|
if (fabs (1.0 - prob) < GSL_DBL_EPSILON) |
|
|
{ |
|
|
return GSL_NEGINF; |
|
|
} |
|
|
printf ("prob is %f\t", prob); |
|
|
tmp = nu / 2.0; |
|
|
d = tmp * M_SQRTPI * gsl_sf_gamma (tmp) / gsl_sf_gamma (tmp + 0.5); |
|
|
method_test = gsl_max (d * prob, 2.0 / nu); |
|
|
/* |
|
|
* There are two possible initial approximations. |
|
|
* Which is used depends on prob and nu. |
|
|
*/ |
|
|
if (method_test > 0.05) |
|
|
{ |
|
|
tmp = prob; |
|
|
a = nu - 0.5; |
|
|
/* gsl_cdf_ugaussian_Pinv(tmp) ? */ |
|
|
x = gsl_cdf_ugaussian_Pinv (tmp); |
|
|
y = inv_cornish_fisher (x, nu); |
|
|
tsqn = -1.0 + exp (a * y * y); |
|
|
} |
|
|
else |
|
|
{ |
|
|
z = pow (prob * d, 1 / tmp); |
|
|
zz = z * z; |
|
|
zzz = z * z * z; |
|
|
tsqn = 1 / z + (nu + 1.0) * (-1.0 + z / (2.0 * (nu + 4.0)) + |
|
|
nu * zz / (3.0 * (nu + 2.0) * (nu + 6.0)) |
|
|
+ nu * (nu + 3.0) * (2.0 * nu * nu + |
|
|
9.0 * nu - |
|
|
2.0) * zzz / (8.0 * |
|
|
(nu + |
|
|
2.0) * |
|
|
(nu + |
|
|
2.0) * |
|
|
(nu + |
|
|
4.0) * |
|
|
(nu + |
|
|
4.0) * |
|
|
(nu + |
|
|
8.0))) |
|
|
/ (nu + 2); |
|
|
|
|
|
} |
|
|
initial_result = sqrt (tsqn * nu); |
|
|
printf ("initial_result = %f \t", initial_result); |
|
|
tcdf = gsl_cdf_tdist_Q (initial_result, nu); |
|
|
tmp = 0.5 * (tcdf - prob); |
|
|
tmp2 = tdist_pdf (initial_result, nu); |
|
|
w = tmp / tmp2; |
|
|
psi = initial_result * (nu + 1.0) / (nu + initial_result * initial_result); |
|
|
psi_prime = (nu + 1.0) * (nu - initial_result * initial_result) / |
|
|
((nu + initial_result * initial_result) * |
|
|
(nu + initial_result * initial_result)); |
|
|
d_psi_dt = 2 * initial_result * nu * (nu + 1.0) / |
|
|
((nu + initial_result * initial_result) * |
|
|
(nu + initial_result * initial_result)); |
|
|
d_psiprime_dt = -2.0 * initial_result * (nu + 1.0) * |
|
|
(3.0 * nu - initial_result * initial_result) / |
|
|
((nu + initial_result * initial_result) * |
|
|
(nu + initial_result * initial_result) * |
|
|
(nu + initial_result * initial_result)); |
|
|
c2 = psi / 2.0; |
|
|
c3 = (2.0 * psi * psi + psi_prime) / 6.0; |
|
|
c4 = 3.0 * psi * (2.0 * psi * psi + psi_prime) + |
|
|
4.0 * psi * d_psi_dt + d_psiprime_dt; |
|
|
c4 /= 24.0; |
|
|
result = initial_result + w + c2 * w * w + |
|
|
c3 * w * w * w + c4 * w * w * w * w; |
|
|
printf ("result = %f\n", result); |
|
|
return result; |
|
|
} |
|
|
#endif |
|