Cumulative Normal Distribution Function in C/C++ with parameters

Question

I want to implement equivalent of matlab normcdf function in C++, I have already found this useful post: Cumulative Normal Distribution Function in C/C++ pointing to this implementation http://www.johndcook.com/cpp_phi.html. But I want it with optional mu and sigma parameters as in matlab.

Is it ok, when I change this:

x = fabs(x)/sqrt(2.0);

to:

x = fabs(x - mu)/sqrt(2.0 * sigma * sigma);

or I am supposed to do something else?

Answer 1

Watch out - you want to save the sign of x-mu, not just of x:

int sign = 1;
if (x < mu)
    sign = -1;
x = fabs(x-mu)/sqrt(2.0*sigma*sigma);

Otherwise your scaling is correct.

Answer 2

You could instead add

x = (x-mu)/fabs(sigma);

as the first line of the function for the correct result without changing the sign check. Ideally you should throw an exception for sigma<=0.0 rather than use fabs here.

On an unrelated note, the implementation that you've linked to is actually a single precision approximation just having had float replaced with double (hence the max error being a whopping 7e-8 ;-).
If you're comfortable translating from VB, West gives an implementation of Hart's much more accurate approximation in figure 2 of http://www.wilmott.com/pdfs/090721_west.pdf.
Unfortunately, as a PDF, you can't copy and paste all the magic numbers so you'll have to be very careful when checking that you've copied them correctly!