Cumulative Normal Distribution Function in C/C++

user_216314 0
Viewed 76343

I was wondering if there were statistics functions built into math libraries that are part of the standard C++ libraries like cmath. If not, can you guys recommend a good stats library that would have a cumulative normal distribution function? Thanks in advance.

More specifically, I am looking to use/create a cumulative distribution function.

distribution math statistics c c++
7 Answers

Here's a stand-alone C++ implementation of the cumulative normal distribution in 14 lines of code.

http://www.johndcook.com/cpp_phi.html

#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";
}
John_D_Cook_251880
user_1538650

From https://en.cppreference.com/w/cpp/numeric/math/erfc

Normal CDF can be calculated as below:

#include <iostream>
#include <cmath>
#include <iomanip>
using namespace std;

double normalCDF(double x) // Phi(-∞, x) aka N(x)
{
    return erfc(-x / sqrt(2))/2;
}

Using 2.0 instead of 2 in the denominator helps in getting decimals instead of integers.

Hope that helps.

user_9844710
Related