19 |
|
|
20 |
#include <config.h> |
#include <config.h> |
21 |
#include <gsl/gsl_cdf.h> |
#include <gsl/gsl_cdf.h> |
22 |
#include <gsl/gsl_sf.h> |
#include <gsl/gsl_sf_gamma.h> |
|
#include <gsl/gsl_sf_result.h> |
|
23 |
#include <gsl/gsl_errno.h> |
#include <gsl/gsl_errno.h> |
24 |
|
|
25 |
/* |
/* |
26 |
* Lower tail. |
* Lower tail. |
27 |
*/ |
*/ |
28 |
double |
double |
29 |
gsl_cdf_F_ratio_P(double x, double nu1, double nu2) |
gsl_cdf_F_ratio_P (double x, double nu1, double nu2) |
30 |
{ |
{ |
31 |
|
double P; |
32 |
int rc; |
double r = nu2 / nu1; |
33 |
double F, nu1_2, nu2_2; |
|
34 |
gsl_sf_result result; |
if (x < r) |
35 |
|
{ |
36 |
F = nu2 / (nu2 + nu1 * x); |
double u = x / (r + x); |
37 |
nu1_2 = nu1 / 2.0; |
|
38 |
nu2_2 = nu2 / 2.0; |
P = gsl_sf_beta_inc (nu1 / 2.0, nu2 / 2.0, u); |
39 |
|
} |
40 |
rc = gsl_sf_beta_inc_e(nu2_2, nu1_2, F, &result); |
else |
41 |
|
{ |
42 |
/* |
double u = r / (r + x); |
43 |
* I don't like subtracting from 1.0 here. It could cause |
|
44 |
* too much cancellation error. A direct approximation for |
P = 1 - gsl_sf_beta_inc (nu2 / 2.0, nu1 / 2.0, u); |
|
* the lower tail would be better. |
|
|
*/ |
|
|
if(rc == GSL_SUCCESS) |
|
|
{ |
|
|
result.val = 1.0 - result.val; |
|
45 |
} |
} |
46 |
|
|
47 |
return result.val; |
return P; |
48 |
} |
} |
49 |
|
|
50 |
/* |
/* |
51 |
* Upper tail. |
* Upper tail. |
52 |
*/ |
*/ |
53 |
double |
double |
54 |
gsl_cdf_F_ratio_Q(double x, double nu1, double nu2) |
gsl_cdf_F_ratio_Q (double x, double nu1, double nu2) |
55 |
{ |
{ |
56 |
int rc; |
double P; |
57 |
double F, nu1_2, nu2_2; |
double r = nu2 / nu1; |
|
gsl_sf_result result; |
|
|
|
|
|
F = nu2 / (nu2 + nu1 * x); |
|
|
nu1_2 = nu1 / 2.0; |
|
|
nu2_2 = nu2 / 2.0; |
|
|
|
|
|
rc = gsl_sf_beta_inc_e(nu2_2, nu1_2, F, &result); |
|
58 |
|
|
59 |
return result.val; |
if (x < r) |
60 |
|
{ |
61 |
|
double u = x / (r + x); |
62 |
|
|
63 |
} |
P = 1 - gsl_sf_beta_inc (nu1 / 2.0, nu2 / 2.0, u); |
64 |
|
} |
65 |
|
else |
66 |
|
{ |
67 |
|
double u = r / (r + x); |
68 |
|
|
69 |
/* |
P = gsl_sf_beta_inc (nu2 / 2.0, nu1 / 2.0, u); |
70 |
* Compute the F distribution function by transforming |
} |
|
* to the incomplete beta function. |
|
|
*/ |
|
|
/* int gsl_cdf_F_e ( double x, double nu1, double nu2, */ |
|
|
/* gsl_cdf_result * result, gsl_cdf_tail_t tail) */ |
|
|
/* { */ |
|
|
/* int rc; */ |
|
|
|
|
|
/* if( tail == GSL_CDF_LOWER ) */ |
|
|
/* { */ |
|
|
/* if( x <= 0.0 ) */ |
|
|
/* { */ |
|
|
/* result->val = 0.0; */ |
|
|
/* result->err = 0.0; */ |
|
|
/* return GSL_SUCCESS; */ |
|
|
/* } */ |
|
|
/* else */ |
|
|
/* { */ |
|
|
/* rc = gsl_F_ratio_cdf_P( x, nu1, nu2, result); */ |
|
|
/* } */ |
|
|
/* } */ |
|
|
/* else if (tail == GSL_CDF_UPPER ) */ |
|
|
/* { */ |
|
|
/* if( x <= 0.0 ) */ |
|
|
/* { */ |
|
|
/* result->val = 1.0; */ |
|
|
/* result->err = 0.0; */ |
|
|
/* return GSL_SUCCESS; */ |
|
|
/* } */ |
|
|
/* else */ |
|
|
/* { */ |
|
|
/* rc = gsl_F_ratio_cdf_Q (x, nu1, nu2, result); */ |
|
|
/* } */ |
|
|
/* } */ |
|
|
/* else */ |
|
|
/* { */ |
|
|
/* return GSL_EINVAL; */ |
|
|
/* } */ |
|
71 |
|
|
72 |
/* return rc; */ |
return P; |
73 |
/* } */ |
} |