17 |
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA. |
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA. |
18 |
*/ |
*/ |
19 |
|
|
|
#include<gsl/gsl_sf_gamma.h> |
|
|
#include<gsl/gsl_sf_result.h> |
|
20 |
#include <gsl/gsl_cdf.h> |
#include <gsl/gsl_cdf.h> |
21 |
#include <gsl/gsl_cdf_result.h> |
#include <gsl/gsl_cdf_result.h> |
22 |
#include <gsl/gsl_errno.h> |
#include <gsl/gsl_errno.h> |
23 |
|
|
24 |
|
#include "gamma.h" |
25 |
|
|
26 |
/* |
/* |
27 |
* Lower tail. |
* Lower tail. |
28 |
*/ |
*/ |
29 |
static int gsl_F_ratio_cdf_P(double x, double nu1, double nu2, |
static int gsl_F_ratio_cdf_P(double x, double nu1, double nu2, |
30 |
gsl_sf_result * result) |
gsl_cdf_result * result) |
31 |
{ |
{ |
32 |
|
|
33 |
int rc; |
int rc; |
37 |
nu1_2 = nu1 / 2.0; |
nu1_2 = nu1 / 2.0; |
38 |
nu2_2 = nu2 / 2.0; |
nu2_2 = nu2 / 2.0; |
39 |
|
|
40 |
rc = gsl_sf_beta_inc_e(nu2_2, nu1_2, F, result); |
rc = gsl_cdf_beta_inc_e(nu2_2, nu1_2, F, result); |
41 |
|
|
42 |
|
/* |
43 |
|
* I don't like subtracting from 1.0 here. It could cause |
44 |
|
* too much cancellation error. A direct approximation for |
45 |
|
* the lower tail would be better. |
46 |
|
*/ |
47 |
if(rc == GSL_SUCCESS) |
if(rc == GSL_SUCCESS) |
48 |
{ |
{ |
49 |
result->val = 1.0 - result->val; |
result->val = 1.0 - result->val; |
56 |
* Upper tail. |
* Upper tail. |
57 |
*/ |
*/ |
58 |
static int gsl_F_ratio_cdf_Q(double x, double nu1, double nu2, |
static int gsl_F_ratio_cdf_Q(double x, double nu1, double nu2, |
59 |
gsl_sf_result * result) |
gsl_cdf_result * result) |
60 |
{ |
{ |
61 |
|
|
62 |
int rc; |
int rc; |
66 |
nu1_2 = nu1 / 2.0; |
nu1_2 = nu1 / 2.0; |
67 |
nu2_2 = nu2 / 2.0; |
nu2_2 = nu2 / 2.0; |
68 |
|
|
69 |
rc = gsl_sf_beta_inc_e(nu2_2, nu1_2, F, &result); |
rc = gsl_cdf_beta_inc_e(nu2_2, nu1_2, F, result); |
70 |
|
|
71 |
return rc; |
return rc; |
72 |
|
|
80 |
gsl_cdf_result * result, gsl_cdf_tail_t tail) |
gsl_cdf_result * result, gsl_cdf_tail_t tail) |
81 |
{ |
{ |
82 |
int rc; |
int rc; |
83 |
gsl_sf_result * sf_result; |
gsl_cdf_result * result; |
84 |
|
|
85 |
if( tail == GSL_CDF_UPPER ) |
if( tail == GSL_CDF_UPPER ) |
86 |
{ |
{ |
90 |
result->err = 0.0; |
result->err = 0.0; |
91 |
return GSL_SUCCESS; |
return GSL_SUCCESS; |
92 |
} |
} |
93 |
rc = gsl_F_ratio_cdf_P( x, nu1, nu2, sf_result); |
else |
94 |
|
{ |
95 |
|
rc = gsl_F_ratio_cdf_P( x, nu1, nu2, result); |
96 |
|
} |
97 |
} |
} |
98 |
else if (tail == GSL_CDF_LOWER ) |
else if (tail == GSL_CDF_LOWER ) |
99 |
{ |
{ |
103 |
result->err = 0.0; |
result->err = 0.0; |
104 |
return rc; |
return rc; |
105 |
} |
} |
106 |
rc = gsl_F_ratio_cdf_Q (x, nu1, nu2, sf_result); |
else |
107 |
|
{ |
108 |
|
rc = gsl_F_ratio_cdf_Q (x, nu1, nu2, result); |
109 |
|
} |
110 |
} |
} |
111 |
else |
else |
112 |
{ |
{ |
113 |
return GSL_EINVAL; |
return GSL_EINVAL; |
114 |
} |
} |
115 |
if ( rc == GSL_SUCCESS ) |
|
|
{ |
|
|
result->val = sf_result->val; |
|
|
result->err = sf_result->err; |
|
|
} |
|
116 |
return rc; |
return rc; |
117 |
} |
} |