/[gsl]/gsl/cdf/cdf_gamma.c
ViewVC logotype

Diff of /gsl/cdf/cdf_gamma.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1.2.6 by jstover, Wed Jan 15 22:29:51 2003 UTC revision 1.1.2.7 by jstover, Mon Feb 17 23:00:08 2003 UTC
# Line 25  Line 25 
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  }  }

Legend:
Removed from v.1.1.2.6  
changed lines
  Added in v.1.1.2.7

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26