Prettify math formula in code

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, particullary 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 i in xrange(len(coeff) - 2, -1, -1):
    poly *= x
    poly += coeff[i]

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.

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)

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

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.

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.

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.