22 |
*/ |
*/ |
23 |
|
|
24 |
#include <config.h> |
#include <config.h> |
25 |
|
#include <math.h> |
26 |
#include <gsl/gsl_cdf.h> |
#include <gsl/gsl_cdf.h> |
27 |
#include <gsl/gsl_sf.h> |
#include <gsl/gsl_math.h> |
28 |
|
#include <gsl/gsl_sf_gamma.h> |
29 |
|
|
30 |
#define BIG_SHAPE 85 |
#define BIG_SHAPE 85 |
31 |
|
|
40 |
* This initial coding of the approximation is giving obviously |
* This initial coding of the approximation is giving obviously |
41 |
* incorrect errors. |
* incorrect errors. |
42 |
*/ |
*/ |
43 |
static double gsl_cdf_g(double x) |
|
44 |
|
static double |
45 |
|
norm_arg (double x, double shape) |
46 |
{ |
{ |
47 |
double val; |
double t; |
48 |
double tmp; |
double arg = x + (1.0 / 3.0) - shape - (0.02 / shape); |
49 |
|
double u = (shape - 0.5) / x; |
50 |
|
|
51 |
if (fabs(x-1.0) < GSL_DBL_EPSILON ) |
if (fabs (u - 1.0) < GSL_DBL_EPSILON) |
52 |
{ |
{ |
53 |
return 0.0; |
t = 0.0; |
54 |
} |
} |
55 |
else if ( fabs(x) < GSL_DBL_EPSILON) |
else if (fabs (u) < GSL_DBL_EPSILON) |
56 |
{ |
{ |
57 |
return 1.0; |
t = 1.0; |
58 |
} |
} |
59 |
else if ( x > 0.0 ) |
else if (u > 0.0) |
60 |
{ |
{ |
61 |
tmp = 1.0-x; |
double v = 1.0 - u; |
62 |
tmp *= tmp; |
t = (1.0 - u * u + 2 * u * log (u)) / (v * v); |
|
val = (1.0-x*x+2*x*log(x))/tmp; |
|
|
return val; |
|
63 |
} |
} |
64 |
else |
else |
65 |
{ |
{ |
66 |
return GSL_EDOM; |
t = GSL_NAN; |
|
} |
|
|
} |
|
|
static double gsl_cdf_norm_arg ( double x, double shape ) |
|
|
{ |
|
|
double val; |
|
|
double tmp; |
|
|
|
|
|
val = x + 1/3 - shape - 0.02/shape; |
|
|
tmp = gsl_cdf_g ( (shape-.5)/x); |
|
|
if (tmp != GSL_EDOM ) |
|
|
{ |
|
|
val *= sqrt ( (1+tmp)/x); |
|
|
return val; |
|
|
} |
|
|
else |
|
|
{ |
|
|
return GSL_EDOM; |
|
67 |
} |
} |
68 |
|
|
69 |
|
arg *= sqrt ((1 + t) / x); |
70 |
|
|
71 |
|
return arg; |
72 |
} |
} |
73 |
|
|
74 |
/* |
/* |
75 |
* Wrapper for the functions that do the work. |
* Wrapper for the functions that do the work. |
76 |
*/ |
*/ |
77 |
double gsl_cdf_gamma_P ( double x, double scale, double shape ) |
double |
78 |
|
gsl_cdf_gamma_P (double x, double scale, double shape) |
79 |
{ |
{ |
80 |
double val; |
double P; |
81 |
double y; |
double y = x / scale; |
82 |
double z; |
|
83 |
int rc; |
if (x <= 0.0) |
|
gsl_sf_result result; |
|
|
|
|
|
if( shape <= 0.0 ) |
|
|
{ |
|
|
return GSL_EDOM; |
|
|
} |
|
|
if ( x <= 0.0 ) |
|
84 |
{ |
{ |
85 |
return 0.0; |
return 0.0; |
86 |
} |
} |
87 |
y = x / scale; |
|
88 |
if( shape < BIG_SHAPE ) |
if (shape < BIG_SHAPE) |
89 |
{ |
{ |
90 |
rc = gsl_sf_gamma_inc_P_e ( shape, y, &result); |
P = gsl_sf_gamma_inc_P (shape, y); |
|
|
|
|
if( rc == GSL_SUCCESS ) |
|
|
{ |
|
|
return result.val; |
|
|
} |
|
91 |
} |
} |
92 |
else |
else |
93 |
{ |
{ |
94 |
/* |
/* Use Peizer and Pratt's normal approximation. */ |
95 |
* Use Peizer and Pratt's normal approximation above. |
|
96 |
*/ |
double z = norm_arg (y, shape); |
97 |
z = gsl_cdf_norm_arg ( y, shape ); |
P = gsl_cdf_ugaussian_P (z); |
|
val = gsl_cdf_gauss_P ( z ); |
|
|
return val; |
|
98 |
} |
} |
99 |
return GSL_FAILURE; |
|
100 |
|
return P; |
101 |
} |
} |
102 |
|
|
103 |
double gsl_cdf_gamma_Q ( double x, double scale, double shape ) |
double |
104 |
|
gsl_cdf_gamma_Q (double x, double scale, double shape) |
105 |
{ |
{ |
106 |
double val; |
double P; |
107 |
double y; |
double y = x / scale; |
|
double z; |
|
|
int rc; |
|
|
gsl_sf_result result; |
|
108 |
|
|
109 |
if ( shape <= 0.0 ) |
if (x <= 0.0) |
|
{ |
|
|
return GSL_EDOM; |
|
|
} |
|
|
if( x <= 0.0 ) |
|
110 |
{ |
{ |
111 |
return 1.0; |
return 1.0; |
112 |
} |
} |
113 |
y = x / scale; |
|
114 |
if ( shape < BIG_SHAPE ) |
if (shape < BIG_SHAPE) |
115 |
{ |
{ |
116 |
rc = gsl_sf_gamma_inc_Q_e ( scale, y, &result); |
P = gsl_sf_gamma_inc_Q (scale, y); |
|
if( rc == GSL_SUCCESS ) |
|
|
{ |
|
|
return result.val; |
|
|
} |
|
117 |
} |
} |
118 |
else |
else |
119 |
{ |
{ |
120 |
/* |
/* |
121 |
* Peizer and Pratt's approximation mentioned above. |
* Peizer and Pratt's approximation mentioned above. |
122 |
*/ |
*/ |
123 |
z = gsl_cdf_norm_arg ( y, shape ); |
double z = norm_arg (y, shape); |
124 |
val = gsl_cdf_gauss_Q ( z ); |
P = gsl_cdf_ugaussian_Q (z); |
|
return val; |
|
125 |
} |
} |
126 |
return GSL_FAILURE; |
|
127 |
|
return P; |
128 |
} |
} |