Cpp, CDF, Phi, Cook, erf(x)

Stackoverflow

-----------------------------------------------------------------------------------------------------------------------

https://stackoverflow.com/questions/2328258/cumulative-normal-distribution-function-in-c-c

https://en.wikipedia.org/wiki/Error_function

https://en.m.wikipedia.org/wiki/Normal_distribution

Theres is no straight function. But since the gaussian error function and its complementary function is related to the normal cumulative distribution function (see here, or here) we can use the implemented c-function erfc (complementary error function):

double normalCDF(double value)
{
   return 0.5 * erfc(-value * M_SQRT1_2);
}

Which considers the relation of erfc(x) = 1-erf(x) with M_SQRT1_2 = √0,5.

I use it for statistical calculations and it works great. No need for using coefficients.

------------------------------------------------------------------------------------------------

Jonh Cook

https://www.johndcook.com/blog/cpp_phi/

https://www.johndcook.com/erf_and_normal_cdf.pdf

http://poivs.tsput.ru/ru/Math/Functions/SpecialFunctions/ErrorFunction

The function Φ(x) is the cumulative density function (CDF) of a standard normal (Gaussian) random variable. It is closely related to the error function erf(x).

Stand-alone C++ implementation of Φ(x)

CDF

#include <cmath>

double phi(double x)
{
    // constants
    double a1 =  0.254829592;
    double a2 = -0.284496736;
    double a3 =  1.421413741;
    double a4 = -1.453152027;
    double a5 =  1.061405429;
    double p  =  0.3275911;

    // Save the sign of x
    int sign = 1;
    if (x < 0)
        sign = -1;
    x = fabs(x)/sqrt(2.0);

    // A&S formula 7.1.26
    double t = 1.0/(1.0 + p*x);
    double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);

    return 0.5*(1.0 + sign*y);
}

void testPhi()
{
    // Select a few input values
    double x[] = 
    {
        -3, 
        -1, 
        0.0, 
        0.5, 
        2.1 
    };

    // Output computed by Mathematica
    // y = Phi[x]
    double y[] = 
    { 
        0.00134989803163, 
        0.158655253931, 
        0.5, 
        0.691462461274, 
        0.982135579437 
    };

        int numTests = sizeof(x)/sizeof(double);

    double maxError = 0.0;
    for (int i = 0; i < numTests; ++i)
    {
        double error = fabs(y[i] - phi(x[i]));
        if (error > maxError)
            maxError = error;
    }

        std::cout << "Maximum error: " << maxError << "\n";
}