25 |
|
|
26 |
#include <config.h> |
#include <config.h> |
27 |
#include <gsl/gsl_cdf.h> |
#include <gsl/gsl_cdf.h> |
28 |
#include "erf.h" |
#include <gsl/gsl_sf.h> |
|
#include "exp.h" |
|
|
#include "log.h" |
|
|
#include "gamma.h" |
|
|
#include "error.h" |
|
29 |
|
|
30 |
/* The dominant part, |
#define BIG_SHAPE 85 |
|
* D(a,x) := x^a e^(-x) / Gamma(a+1) |
|
|
*/ |
|
|
static |
|
|
int |
|
|
gamma_inc_D(const double a, const double x, gsl_cdf_result * result) |
|
|
{ |
|
|
if(a < 10.0) { |
|
|
double lnr; |
|
|
gsl_cdf_result lg; |
|
|
gsl_cdf_lngamma_e(a+1.0, &lg); |
|
|
lnr = a * log(x) - x - lg.val; |
|
|
result->val = exp(lnr); |
|
|
result->err = 2.0 * GSL_DBL_EPSILON * (fabs(lnr) + 1.0) * fabs(result->val); |
|
|
return GSL_SUCCESS; |
|
|
} |
|
|
else { |
|
|
double mu = (x-a)/a; |
|
|
double term1; |
|
|
gsl_cdf_result gstar; |
|
|
gsl_cdf_result ln_term; |
|
|
gsl_cdf_log_1plusx_mx_e(mu, &ln_term); /* log(1+mu) - mu */ |
|
|
gsl_cdf_gammastar_e(a, &gstar); |
|
|
term1 = exp(a*ln_term.val)/sqrt(2.0*M_PI*a); |
|
|
result->val = term1/gstar.val; |
|
|
result->err = 2.0 * GSL_DBL_EPSILON * (fabs(a*ln_term.val) + 1.0) * fabs(result->val); |
|
|
result->err += gstar.err/fabs(gstar.val) * fabs(result->val); |
|
|
return GSL_SUCCESS; |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
/* P series representation. |
|
|
*/ |
|
|
static |
|
|
int |
|
|
gamma_inc_P_series(const double a, const double x, gsl_cdf_result * result) |
|
|
{ |
|
|
const int nmax = 5000; |
|
|
|
|
|
gsl_cdf_result D; |
|
|
int stat_D = gamma_inc_D(a, x, &D); |
|
|
|
|
|
double sum = 1.0; |
|
|
double term = 1.0; |
|
|
int n; |
|
|
for(n=1; n<nmax; n++) { |
|
|
term *= x/(a+n); |
|
|
sum += term; |
|
|
if(fabs(term/sum) < GSL_DBL_EPSILON) break; |
|
|
} |
|
|
|
|
|
result->val = D.val * sum; |
|
|
result->err = D.err * fabs(sum); |
|
|
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
|
|
|
if(n == nmax) |
|
|
GSL_ERROR ("error", GSL_EMAXITER); |
|
|
else |
|
|
return stat_D; |
|
|
} |
|
|
|
|
|
|
|
|
/* Q large x asymptotic |
|
|
*/ |
|
|
static |
|
|
int |
|
|
gamma_inc_Q_large_x(const double a, const double x, gsl_cdf_result * result) |
|
|
{ |
|
|
const int nmax = 5000; |
|
|
|
|
|
gsl_cdf_result D; |
|
|
const int stat_D = gamma_inc_D(a, x, &D); |
|
|
|
|
|
double sum = 1.0; |
|
|
double term = 1.0; |
|
|
double last = 1.0; |
|
|
int n; |
|
|
for(n=1; n<nmax; n++) { |
|
|
term *= (a-n)/x; |
|
|
if(fabs(term/last) > 1.0) break; |
|
|
if(fabs(term/sum) < GSL_DBL_EPSILON) break; |
|
|
sum += term; |
|
|
last = term; |
|
|
} |
|
|
|
|
|
result->val = D.val * (a/x) * sum; |
|
|
result->err = D.err * fabs((a/x) * sum); |
|
|
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
|
|
|
if(n == nmax) |
|
|
GSL_ERROR ("error", GSL_EMAXITER); |
|
|
else |
|
|
return stat_D; |
|
|
} |
|
|
|
|
|
|
|
|
/* Uniform asymptotic for x near a, a and x large. |
|
|
* See [Temme, p. 285] |
|
|
* FIXME: need c1 coefficient |
|
|
*/ |
|
|
static |
|
|
int |
|
|
gamma_inc_Q_asymp_unif(const double a, const double x, gsl_cdf_result * result) |
|
|
{ |
|
|
const double rta = sqrt(a); |
|
|
const double eps = (x-a)/a; |
|
|
|
|
|
gsl_cdf_result ln_term; |
|
|
const int stat_ln = gsl_cdf_log_1plusx_mx_e(eps, &ln_term); /* log(1+eps) - eps */ |
|
|
const double eta = eps * sqrt(-2.0*ln_term.val/(eps*eps)); |
|
31 |
|
|
32 |
gsl_cdf_result erfc; |
/* |
33 |
|
* Use the normal approximation as defined in |
|
double R; |
|
|
double c0, c1; |
|
|
|
|
|
gsl_cdf_erfc_e(eta*M_SQRT2*rta, &erfc); |
|
|
|
|
|
if(fabs(eps) < GSL_ROOT5_DBL_EPSILON) { |
|
|
c0 = -1.0/3.0 + eps*(1.0/12.0 - eps*(23.0/540.0 - eps*(353.0/12960.0 - eps*589.0/30240.0))); |
|
|
c1 = 0.0; |
|
|
} |
|
|
else { |
|
|
double rt_term; |
|
|
rt_term = sqrt(-2.0 * ln_term.val/(eps*eps)); |
|
|
c0 = (1.0 - 1.0/rt_term)/eps; |
|
|
c1 = 0.0; |
|
|
} |
|
|
|
|
|
R = exp(-0.5*a*eta*eta)/(M_SQRT2*M_SQRTPI*rta) * (c0 + c1/a); |
|
|
|
|
|
result->val = 0.5 * erfc.val + R; |
|
|
result->err = GSL_DBL_EPSILON * fabs(R * 0.5 * a*eta*eta) + 0.5 * erfc.err; |
|
|
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
|
|
|
return stat_ln; |
|
|
} |
|
|
|
|
|
|
|
|
/* Continued fraction for Q. |
|
|
* |
|
|
* Q(a,x) = D(a,x) a/x F(a,x) |
|
|
* 1 (1-a)/x 1/x (2-a)/x 2/x (3-a)/x |
|
|
* F(a,x) = ---- ------- ----- -------- ----- -------- ... |
|
|
* 1 + 1 + 1 + 1 + 1 + 1 + |
|
|
* |
|
|
* Hans E. Plesser, 2002-01-22 (hans dot plesser at itf dot nlh dot no): |
|
|
* |
|
|
* Since the Gautschi equivalent series method for CF evaluation may lead |
|
|
* to singularities, I have replaced it with the modified Lentz algorithm |
|
|
* given in |
|
|
* |
|
|
* I J Thompson and A R Barnett |
|
|
* Coulomb and Bessel Functions of Complex Arguments and Order |
|
|
* J Computational Physics 64:490-509 (1986) |
|
|
* |
|
|
* In consequence, gamma_inc_Q_CF_protected() is now obsolete and has been |
|
|
* removed. |
|
|
* |
|
|
* Identification of terms between the above equation for F(a, x) and |
|
|
* the first equation in the appendix of Thompson&Barnett is as follows: |
|
|
* |
|
|
* b_0 = 0, b_n = 1 for all n > 0 |
|
|
* |
|
|
* a_1 = 1 |
|
|
* a_n = (n/2-a)/x for n even |
|
|
* a_n = (n-1)/(2x) for n odd |
|
34 |
* |
* |
35 |
|
* D.B. Peizer and J.W. Pratt. "A Normal Approximation |
36 |
|
* for Binomial, F, Beta, and Other Common, Related Tail |
37 |
|
* Probabilities, I." Journal of the American Statistical |
38 |
|
* Association, volume 63, issue 324, Dec. 1968. pp 1416-1456. |
39 |
*/ |
*/ |
40 |
|
static double gsl_cdf_g(double x) |
|
static |
|
|
int |
|
|
gamma_inc_Q_CF(const double a, const double x, gsl_cdf_result * result) |
|
41 |
{ |
{ |
42 |
const int nmax = 5000; |
double val; |
43 |
const double small = gsl_pow_3 (GSL_DBL_EPSILON); |
double tmp; |
|
|
|
|
gsl_cdf_result D; |
|
|
const int stat_D = gamma_inc_D(a, x, &D); |
|
|
|
|
|
double hn = 1.0; /* convergent */ |
|
|
double Cn = 1.0 / small; |
|
|
double Dn = 1.0; |
|
|
int n; |
|
|
|
|
|
/* n == 1 has a_1, b_1, b_0 independent of a,x, |
|
|
so that has been done by hand */ |
|
|
for ( n = 2 ; n < nmax ; n++ ) { |
|
|
double an; |
|
|
double delta; |
|
|
|
|
|
if(GSL_IS_ODD(n)) |
|
|
an = 0.5*(n-1)/x; |
|
|
else |
|
|
an = (0.5*n-a)/x; |
|
|
|
|
|
Dn = 1.0 + an * Dn; |
|
|
if ( fabs(Dn) < small ) |
|
|
Dn = small; |
|
|
Cn = 1.0 + an/Cn; |
|
|
if ( fabs(Cn) < small ) |
|
|
Cn = small; |
|
|
Dn = 1.0 / Dn; |
|
|
delta = Cn * Dn; |
|
|
hn *= delta; |
|
|
if(fabs(delta-1) < GSL_DBL_EPSILON) break; |
|
|
} |
|
|
|
|
|
result->val = D.val * (a/x) * hn; |
|
|
result->err = D.err * fabs((a/x) * hn); |
|
|
result->err += GSL_DBL_EPSILON * (2.0 + 0.5*n) * fabs(result->val); |
|
|
|
|
|
if(n == nmax) |
|
|
GSL_ERROR ("error", GSL_EMAXITER); |
|
|
else |
|
|
return stat_D; |
|
|
} |
|
44 |
|
|
45 |
/* Useful for small a and x. Handles the subtraction analytically. |
if (fabs(x-1.0) < GSL_DBL_EPSILON ) |
46 |
*/ |
{ |
47 |
static |
return 0.0; |
|
int |
|
|
gamma_inc_Q_series(const double a, const double x, gsl_cdf_result * result) |
|
|
{ |
|
|
double term1; /* 1 - x^a/Gamma(a+1) */ |
|
|
double sum; /* 1 + (a+1)/(a+2)(-x)/2! + (a+1)/(a+3)(-x)^2/3! + ... */ |
|
|
int stat_sum; |
|
|
double term2; /* a temporary variable used at the end */ |
|
|
|
|
|
{ |
|
|
/* Evaluate series for 1 - x^a/Gamma(a+1), small a |
|
|
*/ |
|
|
const double pg21 = -2.404113806319188570799476; /* PolyGamma[2,1] */ |
|
|
const double lnx = log(x); |
|
|
const double el = M_EULER+lnx; |
|
|
const double c1 = -el; |
|
|
const double c2 = M_PI*M_PI/12.0 - 0.5*el*el; |
|
|
const double c3 = el*(M_PI*M_PI/12.0 - el*el/6.0) + pg21/6.0; |
|
|
const double c4 = -0.04166666666666666667 |
|
|
* (-1.758243446661483480 + lnx) |
|
|
* (-0.764428657272716373 + lnx) |
|
|
* ( 0.723980571623507657 + lnx) |
|
|
* ( 4.107554191916823640 + lnx); |
|
|
const double c5 = -0.0083333333333333333 |
|
|
* (-2.06563396085715900 + lnx) |
|
|
* (-1.28459889470864700 + lnx) |
|
|
* (-0.27583535756454143 + lnx) |
|
|
* ( 1.33677371336239618 + lnx) |
|
|
* ( 5.17537282427561550 + lnx); |
|
|
const double c6 = -0.0013888888888888889 |
|
|
* (-2.30814336454783200 + lnx) |
|
|
* (-1.65846557706987300 + lnx) |
|
|
* (-0.88768082560020400 + lnx) |
|
|
* ( 0.17043847751371778 + lnx) |
|
|
* ( 1.92135970115863890 + lnx) |
|
|
* ( 6.22578557795474900 + lnx); |
|
|
const double c7 = -0.00019841269841269841 |
|
|
* (-2.5078657901291800 + lnx) |
|
|
* (-1.9478900888958200 + lnx) |
|
|
* (-1.3194837322612730 + lnx) |
|
|
* (-0.5281322700249279 + lnx) |
|
|
* ( 0.5913834939078759 + lnx) |
|
|
* ( 2.4876819633378140 + lnx) |
|
|
* ( 7.2648160783762400 + lnx); |
|
|
const double c8 = -0.00002480158730158730 |
|
|
* (-2.677341544966400 + lnx) |
|
|
* (-2.182810448271700 + lnx) |
|
|
* (-1.649350342277400 + lnx) |
|
|
* (-1.014099048290790 + lnx) |
|
|
* (-0.191366955370652 + lnx) |
|
|
* ( 0.995403817918724 + lnx) |
|
|
* ( 3.041323283529310 + lnx) |
|
|
* ( 8.295966556941250 + lnx); |
|
|
const double c9 = -2.75573192239859e-6 |
|
|
* (-2.8243487670469080 + lnx) |
|
|
* (-2.3798494322701120 + lnx) |
|
|
* (-1.9143674728689960 + lnx) |
|
|
* (-1.3814529102920370 + lnx) |
|
|
* (-0.7294312810261694 + lnx) |
|
|
* ( 0.1299079285269565 + lnx) |
|
|
* ( 1.3873333251885240 + lnx) |
|
|
* ( 3.5857258865210760 + lnx) |
|
|
* ( 9.3214237073814600 + lnx); |
|
|
const double c10 = -2.75573192239859e-7 |
|
|
* (-2.9540329644556910 + lnx) |
|
|
* (-2.5491366926991850 + lnx) |
|
|
* (-2.1348279229279880 + lnx) |
|
|
* (-1.6741881076349450 + lnx) |
|
|
* (-1.1325949616098420 + lnx) |
|
|
* (-0.4590034650618494 + lnx) |
|
|
* ( 0.4399352987435699 + lnx) |
|
|
* ( 1.7702236517651670 + lnx) |
|
|
* ( 4.1231539047474080 + lnx) |
|
|
* ( 10.342627908148680 + lnx); |
|
|
|
|
|
term1 = a*(c1+a*(c2+a*(c3+a*(c4+a*(c5+a*(c6+a*(c7+a*(c8+a*(c9+a*c10))))))))); |
|
|
} |
|
|
|
|
|
{ |
|
|
/* Evaluate the sum. |
|
|
*/ |
|
|
const int nmax = 5000; |
|
|
double t = 1.0; |
|
|
int n; |
|
|
sum = 1.0; |
|
|
|
|
|
for(n=1; n<nmax; n++) { |
|
|
t *= -x/(n+1.0); |
|
|
sum += (a+1.0)/(a+n+1.0)*t; |
|
|
if(fabs(t/sum) < GSL_DBL_EPSILON) break; |
|
|
} |
|
|
|
|
|
if(n == nmax) |
|
|
stat_sum = GSL_EMAXITER; |
|
|
else |
|
|
stat_sum = GSL_SUCCESS; |
|
|
} |
|
|
|
|
|
term2 = (1.0 - term1) * a/(a+1.0) * x * sum; |
|
|
result->val = term1 + term2; |
|
|
result->err = GSL_DBL_EPSILON * (fabs(term1) + 2.0*fabs(term2)); |
|
|
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
return stat_sum; |
|
|
} |
|
|
|
|
|
|
|
|
/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ |
|
|
|
|
|
/* |
|
|
* Normalized Incomplete Gamma Function |
|
|
* |
|
|
* Q(a,x) = 1/Gamma(a) Integral[ t^(a-1) e^(-t), {t,x,Infinity} ] |
|
|
* |
|
|
* a > 0, x >= 0 |
|
|
* |
|
|
* exceptions: GSL_EDOM |
|
|
*/ |
|
|
static int gamma_inc_Q_e(const double a, const double x, gsl_cdf_result * result) |
|
|
{ |
|
|
if(a <= 0.0) { |
|
|
DOMAIN_ERROR(result); |
|
|
} |
|
|
else if(x <= 0.0) { |
|
|
result->val = 1.0; |
|
|
result->err = 0.0; |
|
|
return GSL_SUCCESS; |
|
|
} |
|
|
else if(x <= 0.5*a) { |
|
|
/* If the series is quick, do that. It is |
|
|
* robust and simple. |
|
|
*/ |
|
|
gsl_cdf_result P; |
|
|
int stat_P = gamma_inc_P_series(a, x, &P); |
|
|
result->val = 1.0 - P.val; |
|
|
result->err = P.err; |
|
|
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
return stat_P; |
|
|
} |
|
|
else if(a >= 1.0e+06 && (x-a)*(x-a) < a) { |
|
|
/* Then try the difficult asymptotic regime. |
|
|
* This is the only way to do this region. |
|
|
*/ |
|
|
return gamma_inc_Q_asymp_unif(a, x, result); |
|
|
} |
|
|
else if(a < 0.2 && x < 5.0) { |
|
|
/* Cancellations at small a must be handled |
|
|
* analytically; x should not be too big |
|
|
* either since the series terms grow |
|
|
* with x and log(x). |
|
|
*/ |
|
|
return gamma_inc_Q_series(a, x, result); |
|
|
} |
|
|
else if(a <= x) { |
|
|
if(x <= 1.0e+06) { |
|
|
/* Continued fraction is excellent for x >~ a. |
|
|
* We do not let x be too large when x > a since |
|
|
* it is somewhat pointless to try this there; |
|
|
* the function is rapidly decreasing for |
|
|
* x large and x > a, and it will just |
|
|
* underflow in that region anyway. We |
|
|
* catch that case in the standard |
|
|
* large-x method. |
|
|
*/ |
|
|
return gamma_inc_Q_CF(a, x, result); |
|
48 |
} |
} |
49 |
else { |
else if ( fabs(x) < GSL_DBL_EPSILON) |
50 |
return gamma_inc_Q_large_x(a, x, result); |
{ |
51 |
|
return 1.0; |
52 |
} |
} |
53 |
} |
else if ( x > 0.0 ) |
54 |
else { |
{ |
55 |
if(a < 0.8*x) { |
tmp = 1.0-x*x; |
56 |
/* Continued fraction again. The convergence |
val = (tmp+2*x*log(x))/(tmp*tmp); |
57 |
* is a little slower here, but that is fine. |
return val; |
|
* We have to trade that off against the slow |
|
|
* convergence of the series, which is the |
|
|
* only other option. |
|
|
*/ |
|
|
return gamma_inc_Q_CF(a, x, result); |
|
58 |
} |
} |
59 |
else { |
else |
60 |
gsl_cdf_result P; |
{ |
61 |
int stat_P = gamma_inc_P_series(a, x, &P); |
return GSL_EDOM; |
|
result->val = 1.0 - P.val; |
|
|
result->err = P.err; |
|
|
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
return stat_P; |
|
62 |
} |
} |
|
} |
|
63 |
} |
} |
64 |
|
static double gsl_cdf_norm_arg ( double x, double y ) |
|
|
|
|
/* |
|
|
* Complementary Normalized Incomplete Gamma Function |
|
|
* |
|
|
* P(a,x) = 1/Gamma(a) Integral[ t^(a-1) e^(-t), {t,0,x} ] |
|
|
* |
|
|
* a > 0, x >= 0 |
|
|
* |
|
|
* exceptions: GSL_EDOM |
|
|
*/ |
|
|
static int gamma_inc_P_e(const double a, const double x, gsl_cdf_result * result) |
|
65 |
{ |
{ |
66 |
if(a <= 0.0 ) { |
double val; |
67 |
DOMAIN_ERROR(result); |
double tmp; |
68 |
} |
|
69 |
else if(x <= 0.0) { |
val = x + 1/3 - y - 0.2/y; |
70 |
result->val = 0.0; |
tmp = gsl_cdf_g ( (y-.5)/x); |
71 |
result->err = 0.0; |
if (tmp != GSL_EDOM ) |
72 |
return GSL_SUCCESS; |
{ |
73 |
} |
val *= sqrt ( (1+tmp)/x); |
74 |
else if(x < 20.0 || x < 0.5*a) { |
return val; |
|
/* Do the easy series cases. Robust and quick. |
|
|
*/ |
|
|
return gamma_inc_P_series(a, x, result); |
|
|
} |
|
|
else if(a > 1.0e+06 && (x-a)*(x-a) < a) { |
|
|
/* Crossover region. Note that Q and P are |
|
|
* roughly the same order of magnitude here, |
|
|
* so the subtraction is stable. |
|
|
*/ |
|
|
gsl_cdf_result Q; |
|
|
int stat_Q = gamma_inc_Q_asymp_unif(a, x, &Q); |
|
|
result->val = 1.0 - Q.val; |
|
|
result->err = Q.err; |
|
|
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
return stat_Q; |
|
|
} |
|
|
else if(a <= x) { |
|
|
/* Q <~ P in this area, so the |
|
|
* subtractions are stable. |
|
|
*/ |
|
|
gsl_cdf_result Q; |
|
|
int stat_Q; |
|
|
if(a > 0.2*x) { |
|
|
stat_Q = gamma_inc_Q_CF(a, x, &Q); |
|
|
} |
|
|
else { |
|
|
stat_Q = gamma_inc_Q_large_x(a, x, &Q); |
|
|
} |
|
|
result->val = 1.0 - Q.val; |
|
|
result->err = Q.err; |
|
|
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
return stat_Q; |
|
|
} |
|
|
else { |
|
|
if((x-a)*(x-a) < a) { |
|
|
/* This condition is meant to insure |
|
|
* that Q is not very close to 1, |
|
|
* so the subtraction is stable. |
|
|
*/ |
|
|
gsl_cdf_result Q; |
|
|
int stat_Q = gamma_inc_Q_CF(a, x, &Q); |
|
|
result->val = 1.0 - Q.val; |
|
|
result->err = Q.err; |
|
|
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
return stat_Q; |
|
75 |
} |
} |
76 |
else { |
else |
77 |
return gamma_inc_P_series(a, x, result); |
{ |
78 |
|
return GSL_EDOM; |
79 |
} |
} |
|
} |
|
|
} |
|
|
|
|
|
/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ |
|
|
|
|
|
#include "eval.h" |
|
|
|
|
|
double gsl_cdf_gamma_inc_P(const double a, const double x) |
|
|
{ |
|
|
EVAL_RESULT(gamma_inc_P_e(a, x, &result)); |
|
80 |
} |
} |
|
|
|
|
double gsl_cdf_gamma_inc_Q(const double a, const double x) |
|
|
{ |
|
|
EVAL_RESULT(gamma_inc_Q_e(a, x, &result)); |
|
|
} |
|
|
|
|
81 |
/* |
/* |
82 |
* Wrapper for the functions that do the work. |
* Wrapper for the functions that do the work. |
83 |
*/ |
*/ |
84 |
int gsl_cdf_gamma_e ( double x, double scale, double shape, |
double gsl_cdf_gamma_P ( double x, double scale, double shape ) |
|
gsl_cdf_tail_t tail, gsl_cdf_result *result) |
|
85 |
{ |
{ |
86 |
int rc = 0; |
double val; |
87 |
|
double y; |
88 |
|
double z; |
89 |
|
int rc; |
90 |
|
gsl_sf_result result; |
91 |
|
|
92 |
if( shape <= 0.0 ) |
if( shape <= 0.0 ) |
93 |
{ |
{ |
94 |
DOMAIN_ERROR(result); |
return GSL_EDOM; |
95 |
} |
} |
96 |
else if ( tail == GSL_CDF_LOWER ) |
y = x / scale; |
97 |
|
if( shape < BIG_SHAPE ) |
98 |
{ |
{ |
99 |
rc = gamma_inc_P_e ( scale, (x/shape), result); |
rc = gsl_sf_gamma_inc_P_e ( shape, y, &result); |
100 |
|
|
101 |
|
if( rc == GSL_SUCCESS ) |
102 |
|
{ |
103 |
|
return result.val; |
104 |
|
} |
105 |
} |
} |
106 |
else if (tail == GSL_CDF_UPPER ) |
else |
107 |
{ |
{ |
108 |
rc = gamma_inc_Q_e ( scale, (x/shape), result); |
/* |
109 |
|
* Use Peizer and Pratt's normal approximation above. |
110 |
|
*/ |
111 |
|
z = norm_arg ( y, shape ); |
112 |
|
val = gsl_cdf_gauss_P ( z ); |
113 |
|
return val; |
114 |
} |
} |
115 |
return rc; |
return GSL_FAILURE; |
116 |
} |
} |
117 |
|
|
118 |
/* |
double gsl_cdf_gamma_Q ( double x, double scale, double shape ) |
|
* No error code. |
|
|
*/ |
|
|
double gsl_cdf_gamma ( double x, double scale, |
|
|
double shape, gsl_cdf_tail_t tail ) |
|
119 |
{ |
{ |
120 |
EVAL_RESULT ( gsl_cdf_gamma_e ( x, scale, shape, tail, &result)); |
double val; |
121 |
|
double y; |
122 |
|
double z; |
123 |
|
int rc; |
124 |
|
gsl_sf_result result; |
125 |
|
|
126 |
|
if ( shape <= 0.0 ) |
127 |
|
{ |
128 |
|
return GSL_EDOM; |
129 |
|
} |
130 |
|
y = x / scale; |
131 |
|
if ( shape < BIG_SHAPE ) |
132 |
|
{ |
133 |
|
rc = gsl_sf_gamma_inc_Q_e ( scale, y, &result); |
134 |
|
if( rc == GSL_SUCCESS ) |
135 |
|
{ |
136 |
|
return result.val; |
137 |
|
} |
138 |
|
} |
139 |
|
else |
140 |
|
{ |
141 |
|
/* |
142 |
|
* Peizer and Pratt's approximation mentioned above. |
143 |
|
*/ |
144 |
|
z = norm_arg ( y, shape ); |
145 |
|
val = gsl_cdf_gauss_Q ( z ); |
146 |
|
return val; |
147 |
|
} |
148 |
|
return GSL_FAILURE; |
149 |
} |
} |