# Prettify math formula in code

I have a function to calculate the normal distribution in Python:

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


But I have to adhere to PEP8 and 80 char margin, hence the ugly \s. How else should I prettify my code?

In mathematical form,

\begin{align*} \textrm{norm_cdf}(z) = 1 - 0.3989423 e^\frac{-z^2}{2} (&1.330274429 t^5 - 1.821255978 t^4 \\ &+ 1.781477937 t^3 - 0.356563782 t^2 + 0.319381530 t ) \end{align*}

where

$$t = \frac{1}{1 + 0.2316419 z}$$

• @200_success likes for the formula image!! – alvas Aug 7 '14 at 16:43
• Jeez comments in Python are ugly... on another note, can you do something like a C-style typedef for your mathematical constants? It would make the code a hell of a lot more readable, imo. Something like typdef 0.319381530 NORMALITY. But you know, pick better names and stuff. For Python, you could just have all your constants at the top of the method (for readability), something like NORMALITY = 0.319381530. That's my preferred method of organizing mathematical functions anyway. – Chris Cirefice Aug 7 '14 at 17:21
• @ChrisCirefice there are not such a thing as typedefs in Python. Saving the coefficients in a list is a more readable approach. – Davidmh Aug 7 '14 at 17:43
• @Davidmh I know that there is no typedef in Python, which is why I suggested just defining a "constant" like NORMALITY = 0.319381530. Personally, I don't like the list approach as suggested in @Davidmh's answer: coeff = [1, 0.319381530, -0.356563782, 1.781477937, -1.821255978, 1.330274429] If there are a lot of constants, using coeff[0], coeff[1] etc. starts to become unreadable in a long function, especially functions that use the constants multiple times. Again, just my opinion. – Chris Cirefice Aug 7 '14 at 17:51
• @ChrisCirefice A list is useful in a loop, the coefficients don't hold any particular meaning by themelves, and they would only be reused in expresions derivated from the polynomial (derivatives, integrals, etc). – Davidmh Aug 7 '14 at 17:58

Let me quote the wonderful book Numerical Recipes in C++ (but also applicable):

We assume that you know enough never to evaluate a polynomial this way:

p=c[0]+c[1]*x+c[2]*x*x+c[3]*x*x*x+c[4]*x*x*x*x;


or (even worse!),

p=c[0]+c[1]*x+c[2]*pow(x,2.0)+c[3]*pow(x,3.0)+c[4]*pow(x,4.0);


Come the (computer) revolution, all persons found guilty of such criminal behavior will be summarily executed, and their programs won't be!

(You can find the page in your edition in the analytical index, under the entry "puns, particularly bad". I love this book.)

There are two reasons not to do that: accuracy and performance. The correct way to evaluate the polynomial is like this:

-t * (0.319381530  +  t * (-0.356563782 + t * (1.781477937 + t * (-1.821255978 + 1.330274429 * t))))


And you can, of course, split at your convenience, as the newlines inside parenthesis are ignored. Remember the PEP: " The preferred place to break around a binary operator is after the operator, not before it."

-t * (0.319381530  +  t * (-0.356563782 +
t * (1.781477937 + t * (-1.821255978 + 1.330274429 * t))))


Another alternative is to save the coefficients in a list:

coeff = [0, 0.319381530, -0.356563782, 1.781477937, -1.821255978, 1.330274429]
poly = coeff[-1]
for c in coeff[-2::-1]:
poly *= x
poly += c


I am doing operations in place to avoid allocating and deallocating memory, but this is only relevant if x is a NumPy array. If you are evaluating on a single number, you can just use instead the nicer expression:

poly = poly * x + coeff[i]


But I would stick with the first one because it is more general.

Of course, the result has to be multiplied by the prefactor:

return 1 - 0.3989423*math.exp(-z*z/2) * poly


Or, if you want to do it in place:

z2 = z * z # Be careful not to modify your input!
z2 *= 0.5  # Multiplication is faster than division.
np.exp(z2, out=z2)

probd = z2 * poly
probd *= -0.3989423
probd += 1
return probd


Bonus track:

If you are applying this function to large arrays (more than a thousand numbers), you may benefit from using the first technique in numexpr:

expr += '1 - 0.3989423* exp(-z * z / 2) * '
expr += '(-t * (0.319381530  +  t * (-0.356563782 +  t * '
expr += '(1.781477937 + t * (-1.821255978 + 1.330274429 * t)))))'
ne.evaluate(expr)


This will compile the expression for you and transparently use as many cores as you have.

• I can see the severe performance penalty of implementing it in any of the discouraged ways, but not the accuracy penalty. Why is it bad for accuracy? I asked it on SO: stackoverflow.com/q/25203873/174365 – Emilio M Bumachar Aug 8 '14 at 12:46
• Maybe I'm missing some context, but it looks like you're making longer, uglier code for the sake of speed, which did not seem like a concern to the OP at all. – raptortech97 Aug 8 '14 at 13:20
• @EmilioMBumachar Depends on the polynomial. But for many polynomials used in practice the naive approach will result in subtracting large similar sized values which increases the absolute error since larger numbers have larger absolute representation errors using floatingpoints. – CodesInChaos Aug 8 '14 at 16:21
• This method of evaluating polynomials is known as "Horner's scheme" (or method depending on where you live). – Emily L. Aug 9 '14 at 17:03
• Instead of using xrange(len(coeff) - 2, -1, -1), and then indexing coeff by i, I'd iterate through coeff[-2::-1] – Francisco Couzo Oct 19 '16 at 6:06

As it turns out, a similar question was asked recently on Math.SE. Rather than , 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 have the full code with a try-except when importing scipy.stats.norm.cdf(), just in case, the users didn't install scipy but the math.erfc would be good. – alvas Aug 7 '14 at 19:08
• but i guess that will weed out the <= py 2.6 users... – alvas Aug 7 '14 at 19:10
• Why not just require or redistribute SciPy if running on Python ≤ 2.6? – 200_success Aug 7 '14 at 19:12
• Or, just steal this public implementation of erf(). – 200_success Aug 7 '14 at 19:16
• @200_success: A & S 26.2.17 explicitly says the valid inputs are non-negative. The code shouldn't be called with negative values, and in an environment where invalid calls are likely, maybe it should include a check. The error bound is listed as 7.5 * 10 ** -8. Is the accuracy specified somehow for math.erfc? I'm not suggesting either algorithm is better or worse, I'm just saying the tradeoffs should be understood in the context of the problem being solved. – GraniteRobert Aug 8 '14 at 13:05

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)

• If this is true then you could define another function that calls the code snippet you described, as some sort of helper method. – Pimgd Aug 7 '14 at 13:37
• I obviously don't have enough knowledge to put this in a function with a meaningful name but I am not even sure this is required. I think it brings the level of abstraction required to make it look like the mathematical objects involved and moving this to a different place might break the mathematical formula apart. – Josay Aug 7 '14 at 13:42
• The only nitpick i have is maybe rename pol to coeff or coefficients to make it more clear, or possibly use the builtin as described below? – C.B. Aug 7 '14 at 15:53

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 books1.

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

• Welcome to StackExchange Code Review! Cool bit of history, thanks for sharing. – Stephen Rauch Mar 20 '17 at 2:44
• I would also recommend improving the formulas here with MathJax. If you don't know how to do that, I can work on it later. – Jamal Mar 20 '17 at 14:32
• hmm I've never worked with Mathjax before, would be interested to see the result Jamal; thanks for the welcome Stephen Rauch – litepresence Mar 20 '17 at 14:41

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.