Prettify math formula in code

As it turns out, a similar question was asked recently on Math.SE. Rather than reinventing-the-wheel, take advantage of built-in functionality in Python.

Your norm_cdf(z) is just a numerical approximation for

$$P(z) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{z} e^{-t^2/2}\ dt = \int_{-\infty}^{z} Z(t)\ dt = \frac{1}{2} \left( 1 + \mathrm{erf}\left( \frac{z}{\sqrt{2}}\right) \right) = \frac{1}{2} \mathrm{erfc}\left( -\,\frac{z}{\sqrt{2}}\right)$$

Therefore, you could just use math.erfc() (available since Python 2.7) and get a more accurate result (especially for very negative values of \$z\$).

import math def norm_cdf(z): return 0.5 * math.erfc(-x / math.sqrt(2)) 

Better yet, just use scipy.stats.norm.cdf()!

I'll only address the so-called "magic numbers" that several reviewers have mentioned.

Sometimes, when you're working in pure mathematics, what seems at first glance to be a "magic number" really isn't. It may be that the numbers themselves are just part of the problem statement. I think the question boils down to this: can you come up with a name that is more descriptive than the number? If there's a good name, you should probably use it.

At first glance, I thought that your numbers were an inherent part of the problem. But when I looked at Abramowitz and Stegun, I saw that the referenced formula has already named your ugly-looking constants. The names are p (which you mentioned in a comment), and b1 through b5. You should use those names in the code, because they create a very clear link to the original formula definition.

When you decided it was a good idea to add the commentp=0.2316419, it was very strong evidence that the number should be named. (And once the code says p=0.2316419, the comment should be removed.)

By the way, it was VERY good of you to include the exact Abramowitz and Stegun reference in the comment.

Instead of math.pow, use the builtin ** operator. You don't need the \s at EOL, because the parentheses surrounding the expression allow it to implicitly span multiple lines. So after making both of those changes, I wind up with:

def norm_cdf(z): """ Use the norm distribution functions as of Gale-Church (1993) srcfile. """ # Equation 26.2.17 from Abramowitz and Stegun (1964:p.932) t = 1.0 / (1 + 0.2316419*z) # t = 1/(1+pz) , p=0.2316419 probdist = 1.0 - ( (0.319381530 * t) + (-0.356563782 * t**2) + (1.781477937 * t**3) + (-1.821255978 * t**4) + (1.330274429 * t**5)) * 0.3989423*math.exp(-z*z/2) return probdist 

I also rearranged one of the multiplications to make the precedence a bit more obvious and readable.

After those changes, the only issue I still see is all of the magic numbers. I don't know how those constants are arrived at, but it may help readability if the constants could be given meaningful names. Sometimes with formulas there isn't really much of a meaningful name to give, though.

Not sure if this helps but you could easily define a function to evaluate the value of a polynom at a given position

def evaluate_polynom(pol, x): return sum(a * math.pow(x, i) for i, a in enumerate(pol)) 

Then

(0.319381530 * t) + (-0.356563782* math.pow(t,2)) + (1.781477937 * math.pow(t,3)) + (-1.821255978* math.pow(t,4)) + (1.330274429 * math.pow(t,5)) 

becomes :

evaluate_polynom([0, 0.319381530, -0.356563782, 1.781477937, -1.821255978, 1.330274429], t) 

At the risk of (further) angering the moderators, I'll advise against following any advice about numerics that originates from any of the Numeric Recipes books¹.

While Horner's method can help in some cases, it's a long ways from a panacea. Rather than fixing the problem of roundoff when summing values, it attempts to avoid that problem. Unfortunately, it's only partially successful in doing so. For most polynomials there will still be inputs for which the results are relatively poor, even at best.

If the values produced by the terms of the polynomial may lead to numeric instability when summing, you might consider generating each individually, then using something like Kahan summation to sum those terms. If you care, this can also give you an error margin along with the sum itself.

Probably better still is to use the Langlois, et al compensated Horner's scheme. At least the last time I looked carefully, this seemed to be pretty much the state of the art in evaluating polynomials. It maintains roughly the same accuracy of result as you'd get from using Horner's scheme using floating point numbers with double the precision (e.g., using a 64-bit double, it gives roughly the same accuracy as Horner's scheme with 128-bit quad-precision floating point, but without nearly the speed penalty that would normally carry). Like Kahan summation, this supports computation of an error bound on the result.

---

_{1. For reference see critiques such as:
http://www.stat.uchicago.edu/~lekheng/courses/302/wnnr/nr.html
http://www.lysator.liu.se/c/num-recipes-in-c.html
I think a more accurate summary than "wonderful" is: "I've found that Numerical Recipes provide just enough information for a person to get himself into trouble, because after reading NR, one thinks that one understands what's going on."}

whereas:

a*x**3 + b*x**2 + c*x = ((a*x + b)*x + c)*x 

h/t @ Emily L. for "Horner" reference:

https://en.wikipedia.org/wiki/Horner%27s_method

and h/t @ Davidmh for noting improvements in computational speed / precision of this method

gale-church cited it like this in 1990:

import math def pnorm(z): t = 1 / (1 + 0.2316419 * z) pd = (1 - 0.3989423 * math.exp(-z * z / 2) * ((((1.330274429 * t - 1.821255978) * t + 1.781477937) * t - 0.356563782) * t + 0.319381530) * t) return pd 

This method conveniently avoids the t^n issue.

citation:

source:

http://www.aclweb.org/anthology/J93-1004

page 21 of 28 in the pdf

page 95 of the journal Computational Linguistics Volume 19, Number 1

I might "prettify" to:

def pnorm(z): t = 1 / (1 + 0.2316419 * z) pd = (1 - 0.3989423 * math.exp(-z * z / 2) * ((((1.330274429 * t - 1.821255978) * t + 1.781477937) * t - 0.356563782) * t + 0.319381530) * t ) return pd 

if you check the

Abromowitz and Stegun, Handbook of Mathematical Functions

page 932 equation 26.2.17

citation:

http://people.math.sfu.ca/~cbm/aands/page_932.htm

you'll find this:

from which we can create a table giving us:

def pnorm(z): p = 0.2316419 b1 = 0.319381530 b2 = -0.356563782 b3 = 1.781477937 b4 = -1.821255978 b5 = 1.330274429 t = 1 / (1 + p * z) pd = (1 - 0.3989423 * math.exp(-z * z / 2) * ((((b5 * t + b4) * t + b3) * t + b2) * t + b1) * t ) return pd 

Then from the previous page; 931 you will find:

Zx = (1/√(2* π))*e(-z*z/2) 

in python:

Zx = (1/math.sqrt(2* math.pi))*math.exp(-z*z/2) 

and we find that (1/√(2* π)) = 0.3989423

also, I think I like this:

t * (b1 + t * (b2 + t * (b3 + t * (b4 + t * b5)))) 

better than:

(((b5 * t + b4) * t + b3) * t + b2) * t + b1) * t 

so then, finally:

import math def pnorm(z): p = 0.2316419 b1 = 0.319381530 b2 = -0.356563782 b3 = 1.781477937 b4 = -1.821255978 b5 = 1.330274429 t = 1 / (1 + p * z) Zx = (1 / math.sqrt(2 * math.pi)) * math.exp(-z * z / 2) pd = Zx * t * (b1 + t * (b2 - t * (b3 + t * (b4 - t * b5)))) return (1 - pd) 

checking my work against the op's

import matplotlib.pyplot as plt import numpy as np import math def norm_cdf(z): """ Use the norm distribution functions as of Gale-Church (1993) srcfile. """ # Equation 26.2.17 from Abramowitz and Stegun (1964:p.932) t = 1.0 / (1+0.2316419*z) # t = 1/(1+pz) , p=0.2316419 probdist = 1 - 0.3989423*math.exp(-z*z/2) * ((0.319381530 * t)+ \ (-0.356563782* math.pow(t,2))+ \ (1.781477937 * math.pow(t,3)) + \ (-1.821255978* math.pow(t,4)) + \ (1.330274429 * math.pow(t,5))) return probdist for z in np.arange (-3,3,0.01): zf = pnorm(z) plt.plot(z,zf, c='red', marker = '.', ms=1) for z in np.arange (-3,3,0.01): zf = norm_cdf(z)+0.1 #offset 0.1 plt.plot(z,zf, c='blue', marker = '.', ms=1) plt.show() plt.pause(0.1) 

I was expecting the Horner method to be faster, so I ran a time test, substituting:

#Zx = (1.0 / math.sqrt(2.0 * math.pi)) * math.exp(-z * z / 2.0) Zx = 0.3989423* math.exp(-z * z / 2.0) 

to make it fair and upping the np.arrange resolution to 0.0001:

t0 = time.time() for z in np.arange (-3,3,0.0001): zf = pnorm(z) t1 = time.time() for z in np.arange (-3,3,0.0001): zf = norm_cdf(z) t2 = time.time() print ('pnorm time : %s' % (t1-t0)) print ('norm_cdf time : %s' % (t2-t1)) 

and the results, spinning my quad core AMD 7950 FM2+ w/ 16G ram pretty hard (albeit with several other apps running)... defied my expectations:

>>> pnorm time : 81.4725670815 norm_cdf time : 80.7865998745 

The Horner method was not faster

I'm gonna be honest here: I totally suck at python.

That said, is it possible to declare some of those numeric literals as constants? That would clean your code up and clarify the code itself a bit more as well.