faster normcdf function

Based on Marsaglia "Evaluating the Normal Distribution", J. of Stat Soft 11, 4 (July 2004).

Description

/*
 * Based on:
 * Marsaglia, George "Evaluating the Normal Distribution",
 * Journal of Statistical Software 11, 4 (July 2004).
 * http://www.jstatsoft.org/
 */

http://www.mathworks.com/matlabcentral/newsreader/view_thread/252933#655297

• MATLAB release: MATLAB 7.6 (R2008a)

[joe]

Fantastic! Thank you very much!

[Takehisa]

Thank you very much!
I think I found a small bug.
The closing brace in line 40 should be before line 39, shouldn't it ?

[Misha Koshelev]

Thank you for your comment. You are right! I uploaded new version. Thx

[Markus Hahn]

 

[Markus Hahn]

It's a nice mex file.
On my machine it's about 4 times faster than normcdf.
However ...
... there is (for the sake of speed) no error checkin: phi(nan) crashes Matlab (at least on my machine).
So I would prefer to use simply 0.5 * erfc(-z ./ sqrt(2)), which is as fast as phi, and much more robust.

[Misha Koshelev]

Thank you very much Dr Hahn for your wonderful comment.

You are correct no error checking is currently done for the sake of speed. In our application of Bayesian clustering if we are computing phi(nan) somewhere then there is indeed something seriously wrong with our code. I agree, a crash may not be the best way to alert a user to such a problem, but at the very least on our application we have not encountered such a case.

I am also quite curious about your speed comment with regards to erfc being just as fast. I cannot vouch for this one hundred percent as it was some time ago but I am fairly sure that, at the time of writing phi, I did find a speed difference between phi and use of erfc as you suggested on my Ubuntu 8.04 amd64 machine. However, as we are computing marginal likelihoods over many cycles (~100,000) of a Gibbs sampler, we are calling this function quite a bit and so even a very small speed increase becomes significant.

Maybe MATLAB has improved their erfc implementation or perhaps you are on a different platform? In any case thank you for your valuable feedback!

Misha