23 |
#include <gsl/gsl_math.h> |
#include <gsl/gsl_math.h> |
24 |
#include <gsl/gsl_sf_gamma.h> |
#include <gsl/gsl_sf_gamma.h> |
25 |
|
|
26 |
|
#include <stdio.h> |
27 |
|
|
28 |
double |
double |
29 |
gsl_cdf_gamma_Pinv (double p, double a, double b) |
gsl_cdf_gamma_Pinv (double p, double a, double b) |
30 |
{ |
{ |
31 |
double x; |
double x, xc; |
32 |
|
|
33 |
if (p == 1.0) |
if (p == 1.0) |
34 |
{ |
{ |
39 |
return 0.0; |
return 0.0; |
40 |
} |
} |
41 |
|
|
|
if (p < 0.05) |
|
|
{ |
|
|
double x0 = exp((lgamma(b) + log(p))/a); |
|
|
x = x0; |
|
|
} |
|
|
else if (p > 0.95) |
|
|
{ |
|
|
double x0 = exp(-log(1-p)-lgamma(b)); |
|
|
x = x0; |
|
|
} |
|
|
else |
|
|
{ |
|
|
double xg = gsl_cdf_ugaussian_Pinv (p); |
|
|
double x0 = (xg < -sqrt(b)) ? b : sqrt (b) * xg + b; |
|
|
x = x0; |
|
|
} |
|
|
|
|
|
{ |
|
|
double dp = p - gsl_cdf_gamma_P (x, 1.0, b); |
|
|
double u = log(x); |
|
|
|
|
|
while (fabs (dp) > 1e-10 * p) |
|
|
{ |
|
|
double dp1, step = dp / exp((b-1.0)*u); |
|
|
|
|
|
do |
|
|
{ |
|
|
u = u + step; |
|
|
x = exp(u); |
|
|
dp1 = p - gsl_cdf_gamma_P (x, 1.0, b); |
|
|
step *= 2.0; |
|
|
} |
|
|
while((dp > 0 && dp1 > 0) || (dp < 0 && dp1 < 0)); |
|
|
|
|
|
dp = dp1; |
|
|
} |
|
|
} |
|
42 |
|
|
43 |
return a * x; |
return a * x; |
44 |
} |
} |