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

Diff of /gsl/cdf/gamma.c

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

revision 1.1.2.1 by bjg, Mon Apr 7 09:52:13 2003 UTC revision 1.1.2.2 by bjg, Mon Jun 23 18:19:31 2003 UTC
# Line 22  Line 22 
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    
# Line 38  Line 40 
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  }  }

Legend:
Removed from v.1.1.2.1  
changed lines
  Added in v.1.1.2.2

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