|
Packit |
67cb25 |
/* cdf/gauss.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2002, 2004 Jason H. Stover.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software Foundation,
|
|
Packit |
67cb25 |
* Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* Computes the cumulative distribution function for the Gaussian
|
|
Packit |
67cb25 |
* distribution using a rational function approximation. The
|
|
Packit |
67cb25 |
* computation is for the standard Normal distribution, i.e., mean 0
|
|
Packit |
67cb25 |
* and standard deviation 1. If you want to compute Pr(X < t) for a
|
|
Packit |
67cb25 |
* Gaussian random variable X with non-zero mean m and standard
|
|
Packit |
67cb25 |
* deviation sd not equal to 1, find gsl_cdf_ugaussian ((t-m)/sd).
|
|
Packit |
67cb25 |
* This approximation is accurate to at least double precision. The
|
|
Packit |
67cb25 |
* accuracy was verified with a pari-gp script. The largest error
|
|
Packit |
67cb25 |
* found was about 1.4E-20. The coefficients were derived by Cody.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* References:
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* W.J. Cody. "Rational Chebyshev Approximations for the Error
|
|
Packit |
67cb25 |
* Function," Mathematics of Computation, v23 n107 1969, 631-637.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* W. Fraser, J.F Hart. "On the Computation of Rational Approximations
|
|
Packit |
67cb25 |
* to Continuous Functions," Communications of the ACM, v5 1962.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* W.J. Kennedy Jr., J.E. Gentle. "Statistical Computing." Marcel Dekker. 1980.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_cdf.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#ifndef M_1_SQRT2PI
|
|
Packit |
67cb25 |
#define M_1_SQRT2PI (M_2_SQRTPI * M_SQRT1_2 / 2.0)
|
|
Packit |
67cb25 |
#endif
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#define SQRT32 (4.0 * M_SQRT2)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* IEEE double precision dependent constants.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* GAUSS_EPSILON: Smallest positive value such that
|
|
Packit |
67cb25 |
* gsl_cdf_gaussian(x) > 0.5.
|
|
Packit |
67cb25 |
* GAUSS_XUPPER: Largest value x such that gsl_cdf_gaussian(x) < 1.0.
|
|
Packit |
67cb25 |
* GAUSS_XLOWER: Smallest value x such that gsl_cdf_gaussian(x) > 0.0.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#define GAUSS_EPSILON (GSL_DBL_EPSILON / 2)
|
|
Packit |
67cb25 |
#define GAUSS_XUPPER (8.572)
|
|
Packit |
67cb25 |
#define GAUSS_XLOWER (-37.519)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#define GAUSS_SCALE (16.0)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double
|
|
Packit |
67cb25 |
get_del (double x, double rational)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xsq = 0.0;
|
|
Packit |
67cb25 |
double del = 0.0;
|
|
Packit |
67cb25 |
double result = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
xsq = floor (x * GAUSS_SCALE) / GAUSS_SCALE;
|
|
Packit |
67cb25 |
del = (x - xsq) * (x + xsq);
|
|
Packit |
67cb25 |
del *= 0.5;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result = exp (-0.5 * xsq * xsq) * exp (-1.0 * del) * rational;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* Normal cdf for fabs(x) < 0.66291
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static double
|
|
Packit |
67cb25 |
gauss_small (const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
unsigned int i;
|
|
Packit |
67cb25 |
double result = 0.0;
|
|
Packit |
67cb25 |
double xsq;
|
|
Packit |
67cb25 |
double xnum;
|
|
Packit |
67cb25 |
double xden;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const double a[5] = {
|
|
Packit |
67cb25 |
2.2352520354606839287,
|
|
Packit |
67cb25 |
161.02823106855587881,
|
|
Packit |
67cb25 |
1067.6894854603709582,
|
|
Packit |
67cb25 |
18154.981253343561249,
|
|
Packit |
67cb25 |
0.065682337918207449113
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
const double b[4] = {
|
|
Packit |
67cb25 |
47.20258190468824187,
|
|
Packit |
67cb25 |
976.09855173777669322,
|
|
Packit |
67cb25 |
10260.932208618978205,
|
|
Packit |
67cb25 |
45507.789335026729956
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
xsq = x * x;
|
|
Packit |
67cb25 |
xnum = a[4] * xsq;
|
|
Packit |
67cb25 |
xden = xsq;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < 3; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
xnum = (xnum + a[i]) * xsq;
|
|
Packit |
67cb25 |
xden = (xden + b[i]) * xsq;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result = x * (xnum + a[3]) / (xden + b[3]);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* Normal cdf for 0.66291 < fabs(x) < sqrt(32).
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static double
|
|
Packit |
67cb25 |
gauss_medium (const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
unsigned int i;
|
|
Packit |
67cb25 |
double temp = 0.0;
|
|
Packit |
67cb25 |
double result = 0.0;
|
|
Packit |
67cb25 |
double xnum;
|
|
Packit |
67cb25 |
double xden;
|
|
Packit |
67cb25 |
double absx;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const double c[9] = {
|
|
Packit |
67cb25 |
0.39894151208813466764,
|
|
Packit |
67cb25 |
8.8831497943883759412,
|
|
Packit |
67cb25 |
93.506656132177855979,
|
|
Packit |
67cb25 |
597.27027639480026226,
|
|
Packit |
67cb25 |
2494.5375852903726711,
|
|
Packit |
67cb25 |
6848.1904505362823326,
|
|
Packit |
67cb25 |
11602.651437647350124,
|
|
Packit |
67cb25 |
9842.7148383839780218,
|
|
Packit |
67cb25 |
1.0765576773720192317e-8
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
const double d[8] = {
|
|
Packit |
67cb25 |
22.266688044328115691,
|
|
Packit |
67cb25 |
235.38790178262499861,
|
|
Packit |
67cb25 |
1519.377599407554805,
|
|
Packit |
67cb25 |
6485.558298266760755,
|
|
Packit |
67cb25 |
18615.571640885098091,
|
|
Packit |
67cb25 |
34900.952721145977266,
|
|
Packit |
67cb25 |
38912.003286093271411,
|
|
Packit |
67cb25 |
19685.429676859990727
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
absx = fabs (x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
xnum = c[8] * absx;
|
|
Packit |
67cb25 |
xden = absx;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < 7; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
xnum = (xnum + c[i]) * absx;
|
|
Packit |
67cb25 |
xden = (xden + d[i]) * absx;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
temp = (xnum + c[7]) / (xden + d[7]);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result = get_del (x, temp);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* Normal cdf for
|
|
Packit |
67cb25 |
* {sqrt(32) < x < GAUSS_XUPPER} union { GAUSS_XLOWER < x < -sqrt(32) }.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static double
|
|
Packit |
67cb25 |
gauss_large (const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int i;
|
|
Packit |
67cb25 |
double result;
|
|
Packit |
67cb25 |
double xsq;
|
|
Packit |
67cb25 |
double temp;
|
|
Packit |
67cb25 |
double xnum;
|
|
Packit |
67cb25 |
double xden;
|
|
Packit |
67cb25 |
double absx;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const double p[6] = {
|
|
Packit |
67cb25 |
0.21589853405795699,
|
|
Packit |
67cb25 |
0.1274011611602473639,
|
|
Packit |
67cb25 |
0.022235277870649807,
|
|
Packit |
67cb25 |
0.001421619193227893466,
|
|
Packit |
67cb25 |
2.9112874951168792e-5,
|
|
Packit |
67cb25 |
0.02307344176494017303
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
const double q[5] = {
|
|
Packit |
67cb25 |
1.28426009614491121,
|
|
Packit |
67cb25 |
0.468238212480865118,
|
|
Packit |
67cb25 |
0.0659881378689285515,
|
|
Packit |
67cb25 |
0.00378239633202758244,
|
|
Packit |
67cb25 |
7.29751555083966205e-5
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
absx = fabs (x);
|
|
Packit |
67cb25 |
xsq = 1.0 / (x * x);
|
|
Packit |
67cb25 |
xnum = p[5] * xsq;
|
|
Packit |
67cb25 |
xden = xsq;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < 4; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
xnum = (xnum + p[i]) * xsq;
|
|
Packit |
67cb25 |
xden = (xden + q[i]) * xsq;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
temp = xsq * (xnum + p[4]) / (xden + q[4]);
|
|
Packit |
67cb25 |
temp = (M_1_SQRT2PI - temp) / absx;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result = get_del (x, temp);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_cdf_ugaussian_P (const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double result;
|
|
Packit |
67cb25 |
double absx = fabs (x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (absx < GAUSS_EPSILON)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result = 0.5;
|
|
Packit |
67cb25 |
return result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (absx < 0.66291)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result = 0.5 + gauss_small (x);
|
|
Packit |
67cb25 |
return result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (absx < SQRT32)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result = gauss_medium (x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (x > 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result = 1.0 - result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (x > GAUSS_XUPPER)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result = 1.0;
|
|
Packit |
67cb25 |
return result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (x < GAUSS_XLOWER)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result = 0.0;
|
|
Packit |
67cb25 |
return result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result = gauss_large (x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (x > 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result = 1.0 - result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_cdf_ugaussian_Q (const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double result;
|
|
Packit |
67cb25 |
double absx = fabs (x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (absx < GAUSS_EPSILON)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result = 0.5;
|
|
Packit |
67cb25 |
return result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (absx < 0.66291)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result = gauss_small (x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (x < 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result = fabs (result) + 0.5;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result = 0.5 - result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (absx < SQRT32)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result = gauss_medium (x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (x < 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result = 1.0 - result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (x > -(GAUSS_XLOWER))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result = 0.0;
|
|
Packit |
67cb25 |
return result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (x < -(GAUSS_XUPPER))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result = 1.0;
|
|
Packit |
67cb25 |
return result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result = gauss_large (x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (x < 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result = 1.0 - result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_cdf_gaussian_P (const double x, const double sigma)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return gsl_cdf_ugaussian_P (x / sigma);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_cdf_gaussian_Q (const double x, const double sigma)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return gsl_cdf_ugaussian_Q (x / sigma);
|
|
Packit |
67cb25 |
}
|