The algorithm is based on formula 7.1.26 from Handbook of Mathematical Functions by Abramowitz and Stegun, affectionately known as "A&S." The algorithm given in A&S is for computing the error function erf(x). The error function is related to the the function Phi by a simple transformation. (The corresponding function in SciPy is scipy.special.erf.)The code he uses in the article is available for download.
Because this code has no dependencies other than the standard math module, it should run anywhere Python runs. It has been tested on Python 2.5, Python 3.0, and IronPython. The code should be accurate to seven decimal places.
Now that we have an implementation of phi(x), we can write a simple wrapper around it to compute normal probabilities.
The main reason John needed to write this stand alone code is because the normal libraries for this kind of numerical computing in Python are SciPy and Numpy. These have various parts written in speed, and therefore don't work with IronPython which is written in C# rather than C.
However, there is one way that these extensions can be made to work with IronPython. This is through Ironclad, a compatibility layer that mimics the Python 2.5 C API and allows Python extensions written in C to be used from IronPython 2. Ironclad is an open source project created by Resolver Systems, to allow the use of numpy with Resolver One - a spreadsheet system written in IronPython and programmable through Python.
As of recent releases of Ironclad it successfully runs thousands of the SciPy and Numpy tests, and many other C extensions can also be used. John invited William Reade, the lead developer of Ironclad, to write a guest post explaining how Ironclad works and what it does:
The upshot of my recent work is that you can now download Ironclad, type ‘import ironclad; import scipy‘ in an IronPython console, and it will Just Work. I am programmer, hear me roar!In the article William explains how to setup Ironclad, gives examples of what works and how to use it, and also discusses the limitations and future directions for the Ironclad project.
Hundreds of tests now pass in both NumPy and SciPy, and I hope that some of you will be inspired to test it against your own requirements. For example, the Gaussian error function has been mentioned a few times on this blog (and, crucially, I have a vague idea of what it actually is), and I can demonstrate that scipy.special.erf works perfectly under Ironclad.
An interesting project occurs to me. It would be fun, but I haven't the bandwidth to attempt it at the moment. It should be easy to write a wrapper around Ironclad so that Python extensions like SciPy and Numpy could be used from IronRuby or even from C# and F#...