The best tools let you focus on your work rather than on the limitations of the tool itself. This is a fairly common metric when comparing the suitability of programming languages or libraries—the best ones to use are those that “get out of your way” and allow you to direct your thoughts and attention to the problem domain. You want to be able to think in pure terms, “How do I solve X?”, instead of the tainted “How do I solve X with Y?”

In numerical methods it is often not possible to ignore the limitations of the tools, especially when the tool is something as basic as floating point representation. A great deal of care needs to be taken to ensure that the implementation of a numerical algorithm will be accurate despite the limits of machine precision. If instead you try to think purely in terms of the algorithm, your results may not be what you expect. The constraints of the floating point representation of Real numbers cannot be ignored.

I recently relearned this lesson when adapting an algorithm for calculating the Normal Quantile function in Python. Wichura’s 1988 algorithm AS241 uses parameters from a 7-term polynomial approximation of the inverse cumulative normal distribution to generate normal deviates corresponding to a given lower-tail probability with error on the order of 1e-16. This is essentially the limit of machine precision available with a 64-bit floating point representation. In the original Fortran algorithm for function PPND16, Wichura writes the multiplication of the terms in reverse order so that the only operations are multiplication and addition:

``````    PPND16 = Q * (((((((A7 * R + A6) * R + A5) * R + A4) * R + A3)
*             * R + A2) * R + A1) * R + A0) /
*             (((((((B7 * R + B6) * R + B5) * R + B4) * R + B3)
*             * R + B2) * R + B1) * R + ONE)
RETURN``````

Translated literally into Python with A and B as lists of length 8 containing the parameters, I started with this:

``````result = Q * (((((((A * R + A) * R + A) * R + A) * R + A) *  \
R + A) * R + A) * R + A) /                             \
(((((((B * R + B) * R + B) * R + B) * R + B) *   \
R + B) * R + B) * R + 1.0)
return result``````

I verified that the results from my Python implementation matched exactly what is produced by the `qnorm` function in R, which uses a C version of Wichura’s algorithm and is certainly peer-reviewed enough to be considered authoritative. I could have stopped there, but looking at the above code I really wanted to clean it up. I mean, look at it, it’s ugly! Since I was implementing this algorithm in Python, why not make it more Pythonic? A couple of list comprehensions will make that jumbled mess look a lot nicer:

``````result = Q * sum([a * pow(R, i) for i, a in enumerate(A)]) /
sum([b * pow(R, i) for i, b in enumerate(B)])``````

There is perhaps a more clever way of doing this using `reduce` and `operator.mul`, but even this is a welcome improvement, aesthetically at any rate, over the original Fortran.

Note that these are algebraically equivalent forms of the two summations. In both the numerator and denominator, there are 8 terms (including a constant) and the result is the sum of the ith term of `A` (or `B`) multiplied by `R` raised to the ith power. Unfortunately, the aesthetically-pleasing Pythonic calculation is not numerically equivalent: there are small, systematic error perturbations across the range of possible inputs.

Using the rpy2 module I am validating my custom Z-quantile function against the `qnorm` function of R, at 1000 p-values ranging from 0.001 to 0.999:

``````    import rpy2, rpy2.robjects
P = range(1, 1000)
for p in P:
p /= 1000.0
rpy2.robjects.r.assign('p', p)
rpy2.robjects.r('z <- qnorm(p)')
z = rpy2.robjects.r('z')
print p, "\t", "{:.20f}".format(quantile_z(p) - z)``````

Whereas the first translation of the algorithm—the ugly one that is faithful to the Fortran original—produces results that are binary exact matches to R’s `qnorm`, the Pythonic version introduces small errors that exhibit an interesting pattern. Here is a tabulation of those error values:

``````Value           Count      Pct
=====           =====      ===
-8.8818e-16       3       0.30%
-6.6613e-16       8       0.80%
-4.4409e-16      32       3.20%
-3.3307e-16      12       1.20%
-2.2204e-16      85       8.51%
-1.6653e-16       5       0.50%
-1.1102e-16      85       8.51%
-8.327e-17        2       0.20%
-5.551e-17       45       4.50%
-4.163e-17        1       0.10%
-2.776e-17       20       2.00%
-2.082e-17        1       0.10%
-1.388e-17       14       1.40%
-6.94e-18         8       0.80%
-5.2e-18          2       0.20%
-3.47e-18         3       0.30%
-1.73e-18         2       0.20%
-8.7e-19          1       0.10%
0.0             345      34.53%
8.7e-19           1       0.10%
1.73e-18          2       0.20%
3.47e-18          3       0.30%
5.2e-18           2       0.20%
6.94e-18          8       0.80%
1.388e-17        12       1.20%
2.082e-17         1       0.10%
2.776e-17        20       2.00%
5.551e-17        56       5.61%
8.327e-17         3       0.30%
1.1102e-16       80       8.01%
1.6653e-16        1       0.10%
2.2204e-16       87       8.71%
3.3307e-16       11       1.10%
3.8858e-16        1       0.10%
4.4409e-16       26       2.60%
6.6613e-16        7       0.70%
8.8818e-16        4       0.40%``````

A plot of the error values against p reveals the systematic pattern:

Clearly there was a reason he explicitly used the ugly expansion of the polynomials in the original code: there can be subtle errors in floating-point arithmetic when using alternate formulations even though they are algebraically equivalent. The Pythonic expression uses the `pow()` function to raise each term to its corresponding power in the series and this method produces different round-off characteristics than Wichura’s method wherein the terms are expanded out from right to left. Wichura must have recognized that exponentiating the r terms is a less numerically stable calculation.