60 |
return y; |
return y; |
61 |
} |
} |
62 |
|
|
63 |
|
/* |
64 |
|
* Use the Cornish-Fisher asymptotic expansion to find a point u such |
65 |
|
* that gsl_cdf_gauss(y) = tcdf(t). |
66 |
|
* |
67 |
|
*/ |
68 |
|
|
69 |
static double |
static double |
70 |
cornish_fisher (double t, double n) |
cornish_fisher (double t, double n) |
71 |
{ |
{ |
100 |
8518.5 |
8518.5 |
101 |
}; |
}; |
102 |
const double coeffs3[4] = { |
const double coeffs3[4] = { |
103 |
.4, |
0.4, |
104 |
3.3, |
3.3, |
105 |
24.0, |
24.0, |
106 |
85.5 |
85.5 |
107 |
}; |
}; |
108 |
|
|
109 |
double m = n - 0.5; |
double a = n - 0.5; |
110 |
double b = 48.0 * m * m; |
double b = 48.0 * a * a; |
111 |
|
|
112 |
double z2 = m * log (1 + (t * t / n)); |
double z2 = a * log (1 + (t * t / n)); |
113 |
double z = sqrt (z2); |
double z = sqrt (z2); |
114 |
|
|
115 |
double p5 = z * poly_eval (coeffs6, 9, z2); |
double p5 = z * poly_eval (coeffs6, 9, z2); |
116 |
double p4 = z * poly_eval (coeffs5, 7, z2); |
double p4 = -z * poly_eval (coeffs5, 7, z2); |
117 |
double p3 = z * poly_eval (coeffs4, 5, z2); |
double p3 = z * poly_eval (coeffs4, 5, z2); |
118 |
double p2 = z * poly_eval (coeffs3, 3, z2); |
double p2 = -z * poly_eval (coeffs3, 3, z2); |
119 |
double p1 = z * (z2 + 3.0); |
double p1 = z * (z2 + 3.0); |
120 |
double p0 = z; |
double p0 = z; |
121 |
|
|
126 |
y = (y / b) + p1; |
y = (y / b) + p1; |
127 |
y = (y / b) + p0; |
y = (y / b) + p0; |
128 |
|
|
129 |
return y; |
if (t < 0) y *= -1; |
|
} |
|
|
|
|
130 |
|
|
131 |
/* |
return y; |
|
* Use the Cornish-Fisher expansion to find |
|
|
* a point u such that gsl_cdf_gauss(u) = tcdf(t). |
|
|
* Approximate the t cdf with gsl_cdf_gauss. |
|
|
* |
|
|
* Valid for: ???? |
|
|
* |
|
|
*/ |
|
|
static double |
|
|
t_cornish_fisher (const double x, const double nu) |
|
|
{ |
|
|
double u = cornish_fisher (x, nu); |
|
|
double q = gsl_cdf_ugaussian_Q (u); |
|
|
|
|
|
return q; |
|
132 |
} |
} |
133 |
|
|
134 |
|
#if 0 |
135 |
/* |
/* |
136 |
* Series approximation for t > 4.0. This needs to be fixed; |
* Series approximation for t > 4.0. This needs to be fixed; |
137 |
* it shouldn't subtract the result from 1.0. A better way is |
* it shouldn't subtract the result from 1.0. A better way is |
172 |
|
|
173 |
return q; |
return q; |
174 |
} |
} |
175 |
|
#endif |
176 |
|
|
177 |
double |
double |
178 |
gsl_cdf_tdist_P (const double x, const double nu) |
gsl_cdf_tdist_P (const double x, const double nu) |
181 |
|
|
182 |
double x2 = x * x; |
double x2 = x * x; |
183 |
|
|
184 |
if (x2 < nu) |
if (nu > 30 && x2 < 10 * nu) |
185 |
|
{ |
186 |
|
double u = cornish_fisher (x, nu); |
187 |
|
P = gsl_cdf_ugaussian_P (u); |
188 |
|
|
189 |
|
return P; |
190 |
|
} |
191 |
|
|
192 |
|
if (x2 < nu) |
193 |
{ |
{ |
194 |
double u = x2 / nu; |
double u = x2 / nu; |
195 |
double eps = u / (1 + u); |
double eps = u / (1 + u); |
264 |
|
|
265 |
|
|
266 |
#if 0 |
#if 0 |
|
else |
|
|
{ |
|
|
if (0 && x2 < nu) |
|
|
{ |
|
|
double q = t_cornish_fisher (x, nu); |
|
|
|
|
|
if (x <= 0.0) |
|
|
{ |
|
|
P = q; |
|
|
} |
|
|
else |
|
|
{ |
|
|
P = 1.0 - q; |
|
|
} |
|
|
} |
|
|
else if (x2 > nu) |
|
|
{ |
|
|
double q = t_large (fabs (x), nu); |
|
|
|
|
|
if (x >= 0.0) |
|
|
{ |
|
|
P = 1 - 0.5 * q; |
|
|
} |
|
|
else |
|
|
{ |
|
|
P = 0.5 * q; |
|
|
} |
|
|
} |
|
|
} |
|
|
#endif |
|
|
|
|
|
#if 0 |
|
|
{ |
|
|
P = normal_approx (x, nu); |
|
|
|
|
|
if (x > 0.0) |
|
|
{ |
|
|
P = 1.0 - P; |
|
|
} |
|
|
} |
|
|
|
|
|
return P; |
|
|
} |
|
|
#endif |
|
|
|
|
|
|
|
|
|
|
|
#if 0 |
|
267 |
/* |
/* |
268 |
* Invert the T distribution. Uses a method as shown in |
* Invert the T distribution. Uses a method as shown in |
269 |
* Statistical Computing, 5.4.2. This method uses an initial |
* Statistical Computing, 5.4.2. This method uses an initial |
326 |
} |
} |
327 |
|
|
328 |
double |
double |
329 |
gsl_cdf_ut_P_inv (double prob, double nu) |
gsl_cdf_ut_Pinv (double prob, double nu) |
330 |
{ |
{ |
331 |
double initial_result; |
double initial_result; |
332 |
double result; |
double result; |
383 |
{ |
{ |
384 |
tmp = prob; |
tmp = prob; |
385 |
a = nu - 0.5; |
a = nu - 0.5; |
386 |
/* gsl_cdf_ugaussian_P_inv(tmp) ? */ |
/* gsl_cdf_ugaussian_Pinv(tmp) ? */ |
387 |
x = gsl_cdf_ugaussian_P_inv (tmp); |
x = gsl_cdf_ugaussian_Pinv (tmp); |
388 |
y = inv_cornish_fisher (x, nu); |
y = inv_cornish_fisher (x, nu); |
389 |
tsqn = -1.0 + exp (a * y * y); |
tsqn = -1.0 + exp (a * y * y); |
390 |
} |
} |