19 |
|
|
20 |
/* |
/* |
21 |
* Computes the cumulative distribution function for the Gaussian distribution |
* Computes the cumulative distribution function for the Gaussian distribution |
22 |
* using a rational function approximation. This approximation is accurate |
* using a rational function approximation. The computation is for the |
23 |
* to at least double precision. I verified the accuracy witha pari-gp |
* standard Normal distribution, i.e., mean 0 and standard |
24 |
* script and found the largest error to be about 1.4E-20. The coefficients |
* deviation 1. If you want to compute Pr(X < t) for a Gaussian random |
25 |
|
* variable X with non-zero mean m and standard deviation sd not equal |
26 |
|
* to 1, find gsl_cdf_gauss ( (t-m)/sd ). This approximation is accurate |
27 |
|
* to at least double precision. The accuracy was verified with a pari-gp |
28 |
|
* script. The largest error found was about 1.4E-20. The coefficients |
29 |
* were derived by Cody. |
* were derived by Cody. |
30 |
* |
* |
31 |
* References: |
* References: |
93 |
/* |
/* |
94 |
* Normal cdf for fabs(x) < 0.66291 |
* Normal cdf for fabs(x) < 0.66291 |
95 |
*/ |
*/ |
96 |
static double gauss_smallx(double x) |
static double gauss_smallx(const double x) |
97 |
{ |
{ |
98 |
double result = 0.0; |
double result = 0.0; |
99 |
double temp = 0.0; |
double temp = 0.0; |
132 |
/* |
/* |
133 |
* Normal cdf for 0.66291 < fabs(x) < sqrt(32). |
* Normal cdf for 0.66291 < fabs(x) < sqrt(32). |
134 |
*/ |
*/ |
135 |
static double gauss_mediumx ( double x ) |
static double gauss_mediumx ( const double x ) |
136 |
{ |
{ |
137 |
int n; |
int n; |
138 |
double temp = 0.0; |
double temp = 0.0; |
187 |
* Normal cdf for |
* Normal cdf for |
188 |
* {sqrt(32) < x < GAUSS_XUPPER} union { GAUSS_XLOWER < x < -sqrt(32) }. |
* {sqrt(32) < x < GAUSS_XUPPER} union { GAUSS_XLOWER < x < -sqrt(32) }. |
189 |
*/ |
*/ |
190 |
static double gauss_grandex ( double x ) |
static double gauss_grandex ( const double x ) |
191 |
{ |
{ |
192 |
int n; |
int n; |
193 |
double result; |
double result; |
235 |
return result; |
return result; |
236 |
} |
} |
237 |
|
|
238 |
double gsl_cdf_gauss(double x) |
double gsl_cdf_gauss(const double x) |
239 |
{ |
{ |
240 |
double retval; |
double retval; |
241 |
double absx; |
double absx; |