/* cdf/cdf_gamma.c * * Copyright (C) 2003 Jason Stover. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or (at * your option) any later version. * * This program is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA. */ /* * Author: J. Stover */ #include #include #include #define BIG_SHAPE 85 /* * Use the normal approximation as defined in * * D.B. Peizer and J.W. Pratt. "A Normal Approximation * for Binomial, F, Beta, and Other Common, Related Tail * Probabilities, I." Journal of the American Statistical * Association, volume 63, issue 324, Dec. 1968. pp 1416-1456. * * This initial coding of the approximation is giving obviously * incorrect errors. */ static double gsl_cdf_g(double x) { double val; double tmp; if (fabs(x-1.0) < GSL_DBL_EPSILON ) { return 0.0; } else if ( fabs(x) < GSL_DBL_EPSILON) { return 1.0; } else if ( x > 0.0 ) { tmp = 1.0-x; tmp *= tmp; val = (1.0-x*x+2*x*log(x))/tmp; return val; } else { return GSL_EDOM; } } 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; } } /* * Wrapper for the functions that do the work. */ double gsl_cdf_gamma_P ( double x, double scale, double shape ) { double val; double y; double z; int rc; gsl_sf_result result; if( shape <= 0.0 ) { return GSL_EDOM; } if ( x <= 0.0 ) { return 0.0; } y = x / scale; if( shape < BIG_SHAPE ) { rc = gsl_sf_gamma_inc_P_e ( shape, y, &result); if( rc == GSL_SUCCESS ) { return result.val; } } else { /* * Use Peizer and Pratt's normal approximation above. */ z = gsl_cdf_norm_arg ( y, shape ); val = gsl_cdf_gauss_P ( z ); return val; } return GSL_FAILURE; } double gsl_cdf_gamma_Q ( double x, double scale, double shape ) { double val; double y; double z; int rc; gsl_sf_result result; if ( shape <= 0.0 ) { return GSL_EDOM; } if( x <= 0.0 ) { return 1.0; } y = x / scale; if ( shape < BIG_SHAPE ) { rc = gsl_sf_gamma_inc_Q_e ( scale, y, &result); if( rc == GSL_SUCCESS ) { return result.val; } } else { /* * Peizer and Pratt's approximation mentioned above. */ z = gsl_cdf_norm_arg ( y, shape ); val = gsl_cdf_gauss_Q ( z ); return val; } return GSL_FAILURE; }