January 2009

Here’s a surprising theorem. Suppose you have a convex set in Rⁿ you pass a plane through its center of mass. How much of the volume of the set can be on one side? No more than about 63%. (Precisely, 1 – 1/e.) This holds for any dimension n and for any direction through the center. I don’t have a reference for this theorem except that it is mentioned near the end of lecture 5 in this course.

Update:  See Splitting a convex set through its center for an illustration and a partial proof.

Sergey Fomel left a valuable comment on my post about computing the error function erf(x). He gave some sample code for exposing C functions to Python via SWIG on Linux. I haven’t used SWIG, but I’ve heard good things about it. I believe Google uses SWIG extensively to make C++ code callable from Python.

Here’s Sergey’s code.

bash$ cat erf.i

%module erf
#include
double erf(double);

bash$ swig -o erf_wrap.c -python erf.i
bash$ gcc -o erf_wrap.os -c -fPIC -I/usr/include/python2.4 erf_wrap.c
bash$ gcc -o _erf.so -shared erf_wrap.os
bash$ python
>>> from erf import erf
>>> erf(1)
0.84270079294971489

C. S. Lewis on the value of reading old books:

Every age has its own outlook. It is specially good at seeing certain truths and specially liable to make certain mistakes. We all, therefore, need the books that will correct the characteristic mistakes of our own period. And that means the old books. All contemporary writers share to some extent the contemporary outlook—even those, like myself, who seem most opposed to it. … To be sure, the books of the future would be just as good a corrective as the books of the past, but unfortunately we cannot get at them.

The quote comes from an introduction he wrote for a translation of On the Incarnation. In the same book, Lewis recommended reading one old book for every contemporary book or two. I agree that reading old books provides perspective on the present, and I do read old books from time to time, but I don’t come close his recommended ratio of old books.

What are some old books you’ve enjoyed or intend to read?

Update: Related post, Old math books

I’ve seen several people ask lately how to compute the distribution (CDF) function for a standard normal random variable, often denoted Φ(x). They want to know how to compute it in Java, or Python, or C++, etc. Every language has its own standard libraries, and in general I recommend using standard libraries. However, sometimes you want to minimize dependencies. Or maybe you want more transparency than your library allows. The code given here is in Python, but it is so compact that it could easily be ported to any other language.

I just posted Python code for computing the error function, erf(x). The normal density Φ(x) is a simple transformation of erf(x). Given code for erf(x), here’s code for Φ(x).

def phi(x):
    return 0.5*( 1.0 + erf(x/math.sqrt(2)) )

After deriving the transformations between erf(x) and Φ(x) several times, including their complements and inverses, I wrote them down to save. See the PDF file Relating Φ and erf.

See also stand alone code for computing the inverse of the standard normal CDF.

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 – (a₁t¹ + a₂t² + a₃t³ + a₄t⁴ + a₅t⁵)exp(-x²), which is absolutely correct. But directly evaluating an nth order polynomial takes O(n²) operations, while the factorization used in the code above uses O(n) operations. This technique is known as Horner’s method.

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.

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

Here is an amazing slide show of Bill Robertson’s miniature woodworking.

Here is a stairway from a miniature house featured in the slide show.

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

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

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 web site 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:

My Twitter accounts
Twitter daily tip news

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. Here’s a plot of the absolute value of the gamma function over the complex plane.

The volcano-like spikes are centered at the poles at non-positive integers. The steep ramp on the right shows how quickly the gamma function increases as the real part of the argument increases.

Related post: generalizing binomial coefficients

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.

I just ran across a nice series of PowerShell posts by Joe Pruitt called the “PowerShell ABC’s,” one post for each letter of the alphabet. He’s up to O so far.

Arithmetic operators
Begin
CmdLet
Debugging
Execution policy
Format operator
Generics
Here string
Is
JavaScript
Keywords
Location
Matching
Numbers
Output

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 + p² + p³ – p⁴ + …

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 p². 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 p²/(1+p). This confirms that 1/(1+p) is greater than 1-p for all p larger than -1.