Take the 2-minute tour ×
Code Review Stack Exchange is a question and answer site for peer programmer code reviews. It's 100% free, no registration required.

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

share|improve this question
1  
@200_success likes for the formula image!! –  alvas Aug 7 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 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 at 17:43
1  
@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 at 17:51
1  
@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 at 17:58

7 Answers 7

up vote 29 down vote accepted
+50

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

share|improve this answer
1  
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 at 12:46
1  
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 at 13:20
    
@raptortech97 I am giving both options. Anyway, the important part is that the evaluation of the polynomial has to be done correctly to be accurate. –  Davidmh Aug 8 at 13:26
3  
@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 at 16:21
1  
This method of evaluating polynomials is known as "Horner's scheme" (or method depending on where you live). –  Emily L. Aug 9 at 17:03

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

share|improve this answer
    
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 at 19:08
    
but i guess that will weed out the <= py 2.6 users... –  alvas Aug 7 at 19:10
    
Why not just require or redistribute SciPy if running on Python ≤ 2.6? –  200_success Aug 7 at 19:12
2  
Or, just steal this public implementation of erf(). –  200_success Aug 7 at 19:16
1  
@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 at 13:05

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)
share|improve this answer
    
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 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 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 at 15:53

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.

share|improve this answer

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.

share|improve this answer

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

share|improve this answer

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.