Cumulative Normal Distribution Function in C/C++

Question

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.

Answer 1

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.

Answer 2

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"; }