27 |
#include <gsl/gsl_math.h> |
#include <gsl/gsl_math.h> |
28 |
#include <gsl/gsl_sf_gamma.h> |
#include <gsl/gsl_sf_gamma.h> |
29 |
|
|
30 |
#define LARGE_B 85 |
#define LARGE_A 85 |
31 |
|
|
32 |
/* |
/* |
33 |
* Use the normal approximation as defined in |
* Use the normal approximation as defined in |
42 |
*/ |
*/ |
43 |
|
|
44 |
static double |
static double |
45 |
norm_arg (double x, double b) |
norm_arg (double x, double a) |
46 |
{ |
{ |
47 |
double t; |
double t; |
48 |
double arg = x + (1.0 / 3.0) - b - (0.02 / b); |
double arg = x + (1.0 / 3.0) - a - (0.02 / a); |
49 |
double u = (b - 0.5) / x; |
double u = (a - 0.5) / x; |
50 |
|
|
51 |
if (fabs (u - 1.0) < GSL_DBL_EPSILON) |
if (fabs (u - 1.0) < GSL_DBL_EPSILON) |
52 |
{ |
{ |
78 |
gsl_cdf_gamma_P (double x, double a, double b) |
gsl_cdf_gamma_P (double x, double a, double b) |
79 |
{ |
{ |
80 |
double P; |
double P; |
81 |
double y = x / a; |
double y = x / b; |
82 |
|
|
83 |
if (x <= 0.0) |
if (x <= 0.0) |
84 |
{ |
{ |
85 |
return 0.0; |
return 0.0; |
86 |
} |
} |
87 |
|
|
88 |
if (b < LARGE_B) |
if (b < LARGE_A) |
89 |
{ |
{ |
90 |
P = gsl_sf_gamma_inc_P (b, y); |
P = gsl_sf_gamma_inc_P (a, y); |
91 |
} |
} |
92 |
else |
else |
93 |
{ |
{ |
94 |
/* Use Peizer and Pratt's normal approximation. */ |
/* Use Peizer and Pratt's normal approximation. */ |
95 |
|
|
96 |
double z = norm_arg (y, b); |
double z = norm_arg (y, a); |
97 |
P = gsl_cdf_ugaussian_P (z); |
P = gsl_cdf_ugaussian_P (z); |
98 |
} |
} |
99 |
|
|
104 |
gsl_cdf_gamma_Q (double x, double a, double b) |
gsl_cdf_gamma_Q (double x, double a, double b) |
105 |
{ |
{ |
106 |
double P; |
double P; |
107 |
double y = x / a; |
double y = x / b; |
108 |
|
|
109 |
if (x <= 0.0) |
if (x <= 0.0) |
110 |
{ |
{ |
111 |
return 1.0; |
return 1.0; |
112 |
} |
} |
113 |
|
|
114 |
if (b < LARGE_B) |
if (a < LARGE_A) |
115 |
{ |
{ |
116 |
P = gsl_sf_gamma_inc_Q (b, y); |
P = gsl_sf_gamma_inc_Q (a, y); |
117 |
} |
} |
118 |
else |
else |
119 |
{ |
{ |
120 |
/* |
/* |
121 |
* Peizer and Pratt's approximation mentioned above. |
* Peizer and Pratt's approximation mentioned above. |
122 |
*/ |
*/ |
123 |
double z = norm_arg (y, b); |
double z = norm_arg (y, a); |
124 |
P = gsl_cdf_ugaussian_Q (z); |
P = gsl_cdf_ugaussian_Q (z); |
125 |
} |
} |
126 |
|
|