/[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.2 by bjg, Mon Jun 23 18:19:31 2003 UTC revision 1.1.2.3 by bjg, Tue Jun 24 14:13:15 2003 UTC
# Line 27  Line 27 
27  #include <gsl/gsl_math.h>  #include <gsl/gsl_math.h>
28  #include <gsl/gsl_sf_gamma.h>  #include <gsl/gsl_sf_gamma.h>
29    
30  #define BIG_SHAPE 85  #define LARGE_B 85
31    
32  /*  /*
33   * Use the normal approximation as defined in   * Use the normal approximation as defined in
# Line 42  Line 42 
42   */   */
43    
44  static double  static double
45  norm_arg (double x, double shape)  norm_arg (double x, double b)
46  {  {
47    double t;    double t;
48    double arg = x + (1.0 / 3.0) - shape - (0.02 / shape);    double arg = x + (1.0 / 3.0) - b - (0.02 / b);
49    double u = (shape - 0.5) / x;    double u = (b - 0.5) / x;
50    
51    if (fabs (u - 1.0) < GSL_DBL_EPSILON)    if (fabs (u - 1.0) < GSL_DBL_EPSILON)
52      {      {
# Line 75  norm_arg (double x, double shape) Line 75  norm_arg (double x, double shape)
75   * Wrapper for the functions that do the work.   * Wrapper for the functions that do the work.
76   */   */
77  double  double
78  gsl_cdf_gamma_P (double x, double scale, double shape)  gsl_cdf_gamma_P (double x, double a, double b)
79  {  {
80    double P;    double P;
81    double y = x / scale;    double y = x / a;
82    
83    if (x <= 0.0)    if (x <= 0.0)
84      {      {
85        return 0.0;        return 0.0;
86      }      }
87    
88    if (shape < BIG_SHAPE)    if (b < LARGE_B)
89      {      {
90        P = gsl_sf_gamma_inc_P (shape, y);        P = gsl_sf_gamma_inc_P (b, y);
91      }      }
92    else    else
93      {      {
94        /* Use Peizer and Pratt's normal approximation. */        /* Use Peizer and Pratt's normal approximation. */
95    
96        double z = norm_arg (y, shape);        double z = norm_arg (y, b);
97        P = gsl_cdf_ugaussian_P (z);        P = gsl_cdf_ugaussian_P (z);
98      }      }
99    
# Line 101  gsl_cdf_gamma_P (double x, double scale, Line 101  gsl_cdf_gamma_P (double x, double scale,
101  }  }
102    
103  double  double
104  gsl_cdf_gamma_Q (double x, double scale, double shape)  gsl_cdf_gamma_Q (double x, double a, double b)
105  {  {
106    double P;    double P;
107    double y = x / scale;    double y = x / a;
108    
109    if (x <= 0.0)    if (x <= 0.0)
110      {      {
111        return 1.0;        return 1.0;
112      }      }
113    
114    if (shape < BIG_SHAPE)    if (b < LARGE_B)
115      {      {
116        P = gsl_sf_gamma_inc_Q (scale, y);        P = gsl_sf_gamma_inc_Q (a, y);
117      }      }
118    else    else
119      {      {
120        /*        /*
121         * Peizer and Pratt's approximation mentioned above.         * Peizer and Pratt's approximation mentioned above.
122         */         */
123        double z = norm_arg (y, shape);        double z = norm_arg (y, b);
124        P = gsl_cdf_ugaussian_Q (z);        P = gsl_cdf_ugaussian_Q (z);
125      }      }
126    

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

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