44 |
|
|
45 |
#include "beta_inc.c" |
#include "beta_inc.c" |
46 |
|
|
47 |
|
#if 0 |
48 |
/* |
/* |
49 |
* First approximation is also found in Abramowitz |
* First approximation is also found in Abramowitz |
50 |
* and Stegun. This is used only for small values of t |
* and Stegun. This is used only for small values of t |
53 |
static double |
static double |
54 |
t_small (const double t, const double nu) |
t_small (const double t, const double nu) |
55 |
{ |
{ |
56 |
/* Final result is always a small amount to be added to or |
/* Compute beta_inc when t^2 < n. Final result is always a small |
57 |
subtracted from 0.5 so we only need absolute accuracy of |
amount to be added to or subtracted from 0.5 so we only need |
58 |
O(DBL_EPS) -- this information could be used to optimise inside |
absolute accuracy of O(DBL_EPS) -- this information could be used |
59 |
beta_inc (currently it isn't) */ |
to optimise inside beta_inc (currently it isn't) */ |
60 |
|
|
61 |
double u = (t * t) / nu; |
double u = (t * t) / nu; |
62 |
double eps = u / (1 + u); |
double eps = u / (1 + u); |
68 |
static double |
static double |
69 |
t_large (const double t, const double nu) |
t_large (const double t, const double nu) |
70 |
{ |
{ |
71 |
/* Final result may be needed to full relative precision, depending |
/* Compute beta_inc when t^2 > n. Final result may be needed to full |
72 |
on whether we are computing the upper or lower tail. Again this |
relative precision, depending on whether we are computing the |
73 |
information could be used to optimise inside beta_inc. */ |
upper or lower tail. Again this information could be used to |
74 |
|
optimise inside beta_inc. */ |
75 |
|
|
76 |
double v = nu / (t * t); |
double v = nu / (t * t); |
77 |
double eps = v / (1 + v); |
double eps = v / (1 + v); |
80 |
return q; |
return q; |
81 |
} |
} |
82 |
|
|
83 |
|
#endif |
84 |
|
|
85 |
static double |
static double |
86 |
poly_eval (const double c[], unsigned int n, double x) |
poly_eval (const double c[], unsigned int n, double x) |
87 |
{ |
{ |
227 |
|
|
228 |
double x2 = x * x; |
double x2 = x * x; |
229 |
|
|
230 |
if (nu <= 20.0) |
if (1 || nu <= 20.0) |
231 |
{ |
{ |
232 |
if (x2 < nu) |
if (x2 < nu) |
233 |
{ |
{ |
234 |
double q = t_small (fabs (x), nu); |
double u = x2 / nu; |
235 |
|
double eps = u / (1 + u); |
236 |
|
|
237 |
if (x >= 0.0) |
if (x >= 0) |
238 |
{ |
{ |
239 |
P = 0.5 * (1.0 + q); |
P = beta_inc_AXPY (0.5, 0.5, 0.5, nu / 2.0, eps); |
240 |
} |
} |
241 |
else |
else |
242 |
{ |
{ |
243 |
P = 0.5 * (1.0 - q); |
P = beta_inc_AXPY (-0.5, 0.5, 0.5, nu / 2.0, eps); |
244 |
} |
} |
245 |
} |
} |
246 |
else |
else |
247 |
{ |
{ |
248 |
double q = t_large (fabs (x), nu); |
double v = nu / (x * x); |
249 |
|
double eps = v / (1 + v); |
250 |
|
|
251 |
|
if (x >= 0) |
252 |
|
{ |
253 |
|
P = beta_inc_AXPY (-0.5, 1.0, nu / 2.0, 0.5, eps); |
254 |
|
} |
255 |
|
else |
256 |
|
{ |
257 |
|
P = beta_inc_AXPY (0.5, 0.0, nu / 2.0, 0.5, eps); |
258 |
|
} |
259 |
|
} |
260 |
|
} |
261 |
|
#if 0 |
262 |
|
else |
263 |
|
{ |
264 |
|
if (0 && x2 < nu) |
265 |
|
{ |
266 |
|
double q = t_cornish_fisher (x, nu); |
267 |
|
|
268 |
|
if (x <= 0.0) |
269 |
|
{ |
270 |
|
P = q; |
271 |
|
} |
272 |
|
else |
273 |
|
{ |
274 |
|
P = 1.0 - q; |
275 |
|
} |
276 |
|
} |
277 |
|
else if (x2 > nu) |
278 |
|
{ |
279 |
|
double q = t_large (fabs (x), nu); |
280 |
|
|
281 |
if (x >= 0.0) |
if (x >= 0.0) |
282 |
{ |
{ |
283 |
P = 1 - 0.5 * q; |
P = 1 - 0.5 * q; |
288 |
} |
} |
289 |
} |
} |
290 |
} |
} |
291 |
else |
#endif |
|
{ |
|
|
double q = t_cornish_fisher (x, nu); |
|
|
|
|
|
if (x <= 0.0) |
|
|
{ |
|
|
P = q; |
|
|
} |
|
|
else |
|
|
{ |
|
|
P = 1.0 - q; |
|
|
} |
|
|
} |
|
292 |
#if 0 |
#if 0 |
293 |
{ |
{ |
294 |
P = normal_approx (x, nu); |
P = normal_approx (x, nu); |
303 |
return P; |
return P; |
304 |
} |
} |
305 |
|
|
306 |
|
#if 0 |
307 |
double |
double |
308 |
gsl_cdf_t_Q (const double x, const double nu) |
gsl_cdf_t_Q (const double x, const double nu) |
309 |
{ |
{ |
520 |
printf ("result = %f\n", result); |
printf ("result = %f\n", result); |
521 |
return result; |
return result; |
522 |
} |
} |
523 |
|
#endif |