Prettify math formula in code

Question

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}$$

Answer 1

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.

Answer 2

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()!

Answer 3

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.

Answer 4

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.

Answer 5

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)

Answer 6

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."}

Answer 7

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

Answer 8

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.