Stand-alone error function erf(x)

The question came up on StackOverflow this morning how to compute the error function erf(x) in Python. The standard answer for how to compute anything numerical in Python is “Look in SciPy.” However, this person didn’t want to take on the dependence on SciPy. I’ve seen variations on this question come up in several different contexts lately, including questions about computing the normal distribution function, so I thought I’d write up a solution.

Here’s a Python implementation of erf(x) based on formula 7.1.26 from A&S. The maximum error is below 1.5 × 10-7.

import math

def erf(x):
    # constants
    a1 =  0.254829592
    a2 = -0.284496736
    a3 =  1.421413741
    a4 = -1.453152027
    a5 =  1.061405429
    p  =  0.3275911

    # Save the sign of x
    sign = 1
    if x < 0:
        sign = -1
    x = abs(x)

    # A & S 7.1.26
    t = 1.0/(1.0 + p*x)
    y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*math.exp(-x*x)

    return sign*y

This problem is typical in two ways: A&S has a solution, and you’ve got to know a little background before you can use it.

The formula given in A&S is only good for x ≥ 0. That’s no problem if you know that the error function is an odd function, i.e. erf(-x) = -erf(x). But if you’re an engineer who has never heard of the error function but needs to use it, it may take a while to figure out how to handle negative inputs.

One other thing that someone just picking up A&S might not know is the best way to evaluate polynomials. The formula appears as 1 – (a1t1 + a2t2 + a3t3 + a4t4 + a5t5)exp(-x2), which is absolutely correct. But directly evaluating an nth order polynomial takes O(n2) operations, while the factorization used in the code above uses O(n) operations. This technique is known as Horner’s method.

Draw a bigger picture

When I was in college, I had the pleasure of taking a couple math classes taught by Mike Starbird. One of the things he told us about problem solving was this: when you’re stuck, draw a picture. Good advice, though hardly original. But then he said something else: If you’re still stuck, draw a bigger picture.

When you draw small, you think small. It’s surprising how much your thinking opens up when you draw a bigger picture. You can draw some little diagram and think “that didn’t help.” But when you take that same unhelpful diagram and make it much larger, it may be just what you need.

Dating medieval manuscripts via DNA

A couple days ago the 60 Second Science podcast had an interesting story on dating medieval manuscripts. It turns out you can take DNA samples from the animal products in the pages. The hope is that the pages in undated manuscripts will have the same DNA signature as extant dated manuscripts, suggesting that the books were written around the same time.

Transcript

How to display side-by-side figures in LaTeX

The other day I tried putting two graphs side-by-side by setting R’s mfrow parameter. That made the graphs look odd because not everything was scaled proportionately. Then someone told me I could use a LaTeX command to place my original graphs side-by-side rather than creating a new image.

The trick is to use the subfigure package. Include the directive usepackage{subfigure} at the top of your file, then use code something like the following.

\begin{figure}
\centering
\mbox{\subfigure{\includegraphics[width=3in]{fig1.pdf}}\quad
\subfigure{\includegraphics[width=3in]{fig2.pdf} }}
\caption{Text pertaining to both graphs ...} \label{fig12}
end{figure}

If you put text in square brackets immediately after the \subfigure command, that text will be a caption for the corresponding sub-figure. For example, the LaTeX code above might be changed to include the following.

...\subfigure[Description of left graph]{\includegraphics...

How to insert graphics in Twitter messages

I saw a couple interesting messages from John Udell on Twitter yesterday.

screen shot of Twitter posts

Apparently Bill Zeller had the clever idea of using Unicode characters to put sparklines in Twitter messages. It didn’t occur to me that Twitter would accept Unicode. I’m sure the intent of Unicode support is multi-lingual text, not clever hacks such as creating sparklines, but this is fun to play around with.

The sparkline above was created using block element symbols. I tried a few other Unicode characters. Some worked, but most didn’t.  The Twitter protocol supports Unicode, but particular Twitter clients may not have fonts for displaying some characters. For example, I found some characters would display correctly when I went to twitter.com that didn’t display from the Twhirl client.

To enter a Unicode character, you can find out the numeric value here. Once you know the value, here are tips for how to insert Unicode characters in Windows and Linux applications. You might, for example, create the symbols you want using Microsoft Word and then paste them into your Twitter client.

Update (8 January 2010):

Here are some examples of how the differently same tweet may appear on the same computer using some screen shots I took this morning.

First, here’s the view from Twitter’s website using Firefox 3.5.7 on Windows:

Here’s Firefox 3.5.3 on Ubuntu 9.10, nearly the same but a couple characters are missing.

Here’s the view from IE 8 on Windows:

And now Safari on Windows:

Here’s the view from TweetDeck on Windows with its default font:

And now TweetDeck on Windows with its international font:

I imagine TweetDeck would look the same on other operating systems since Adobe Air is largely self-contained.

Related links:

The gamma function

The gamma function Γ(x) is the most important function not on a calculator. It comes up constantly in math. In some areas, such as probability and statistics, you will see the gamma function more often than other functions that are on a typical calculator, such as trig functions.

The gamma function extends the factorial function to real numbers. Since factorial is only defined on non-negative integers, there are many ways you could define factorial that would agree on the integers and disagree elsewhere. But everyone agrees that the gamma function is “the” way to extend factorial. Actually, the gamma function Γ(x) does not extend factorial, but Γ(x+1) does. Shifting the definition over by one makes some equations simpler, and that’s the definition that has caught on.

In a sense, Γ(x+1) is the unique way to generalize factorial. Harald Bohr and Johannes Mollerup proved that it is the only log-convex function that agrees with factorial on the non-negative integers. That’s somewhat satisfying, except why should we look for log-convex functions? Log-convexity is very useful property to have, and a natural one for a function generalizing the factorial.

Here’s a plot of the logarithms of the first few factorial values.

The points nearly fall on a straight line, but if you look closely, the points bend upward slightly. If you have trouble seeing this, imagine a line connecting the first and last dot. This line would lie above the dots in between. This suggests that a function whose graph passes through all the dots should be convex.

Here’s what a graph showing what the gamma function looks like over the real line.

The gamma function is finite except for non-positive integers. It goes to +∞ at zero and negative even integers and to -∞ at negative odd integers. The gamma function can also be uniquely extended to an analytic function on the complex plane. The only singularities are the poles on the real axis.

Related posts

Publish or perish

From @divbyzero:

“If ‘publish or perish’ were really true, Leonhard Euler would still be alive.”—Eric Bach

(Leonhard Euler (1707-1783) published more papers than any mathematician in history.)

Correction: Apparently Euler wrote the most pages in math journals, but Paul Erdős wrote more individual papers.

Sales tax included

Sometimes the wrong answer is more interesting than the right answer. If the wrong approach almost works, it may be fun to understand why.

Suppose you want to figure a price so that the final price including tax has a given value. For example, say you want a T-shirt to sell for $10 after adding 8% sales tax. A common mistake would be to subtract 8% from $10 and sell the shirt for $9.20 plus tax. But this makes the price with tax $9.94. Observe two things: (1) the result was wrong, but (2) it wasn’t far off.

The correct solution would be to divide $10 by 1.08. Then when we add 8%, i.e. multiply by 1.08, we get $10. That says we should price the shirt at $9.26 before tax. That explains (1). But what about (2), that is, why was the result nearly correct? Subtracting a percentage p amounts to multiplying by (1-p). But we should have multiplied by 1/(1+p). We’ve stumbled on the fact that for small p, 1/(1+p) approximately equals 1 – p.

With a little experimentation we might discover a couple more things. First, we notice that multiplying by 1-p always gives a smaller result than multiplying by 1/(1+p), i.e. the common mistake always gets the price too low. Also, with a little experimentation we might notice that the difference between 1/(1+p) and (1-p) gets smaller as p gets smaller. Not only that, making p a little smaller makes the difference between 1/(1+p) and (1-p) a lot smaller. In summary, we noticed that 1/(1+p) – (1-p) is

  1. small,
  2. positive, and
  3. gets small faster than p gets small.

What accounts for these observations? The Taylor series for 1/(1+p) says that for |p| < 1,

1/(1+p) = 1 – p + p2 + p3p4 + …

The error in truncating a Taylor series is roughly equal to the first term that was left out, so the error in approximating 1/(1+p) by (1-p) is roughly p2. That explains why the difference between 1/(1+p) and (1-p) is small, positive, and decreases with p. For example, we should expect that cutting p in half reduces the difference by a factor of four.

The Taylor series argument only necessarily holds for p sufficiently small. If we go back and calculate 1/(1+p) – (1-p) directly we find it’s p2/(1+p). This confirms that 1/(1+p) is greater than 1-p for all p larger than -1.