the second degree. If knowledge is available from other sources about the

functional form that is to be expected, a direct least-squares ¬t to such a

functional form is to be preferred above a ¬tting spline. The latter will

always choose the smoothest curve, statistically consistent with the data.

536 Splines for everything

7

original function

curv = 237, χ2/n = 0.001

6

curv = 0.70, χ2/n = 1.02

curv = 0.001, χ2/n = 3.06

5

4

3

2

1

0

’1

’2

’6 ’4 ’2 0 2 4 6

Figure 19.4 Fitting splines approximating noisy data points. Three curves were

generated with di¬erent weighting factors (1000; 0.2; 0.001). The results range

from a spline through the points to an almost straight line. The “best” ¬t has an

average chi-square deviation close to 1.

What we are looking for is the smoothest curve consistent not only with the

data, but with all the knowledge we may possess about the expected curve.

19.4 Fitting distribution functions

Suppose we wish to generate a probability distribution function ρ(x) from

a number of statistical samples from that distribution. For example, we

wish to generate a radial distribution function between two speci¬ed species

of atoms, or an angular distribution of relative orientations of molecular

dipoles. If there are su¬cient independent data, this task simply reduces

to constructing a histogram by sorting the data into bins x0 , . . . , xn . The

number of samples ni in each bin follows multinomial statistics: when the

probability to ¬nd a sample in the i-th bin is pi , then the joint probability

to ¬nd the distribution n1 , . . . , nn , given the total number of samples N , is

given by

ni

n pi

P (n1 , . . . , nn ) = N ! Πi=1 . (19.40)

ni !

19.4 Fitting distribution functions 537

The (marginal) probability to ¬nd exactly ni samples in the i-th bin, irre-

spective of the occupancy of other bins, follows binomial statistics:

N!

pni (1 ’ pi )N ’ni .

P (ni |N ) = (19.41)

ni !(N ’ ni )! i

For a su¬cient number of bins, the probability to ¬nd a sample in one

particular bin becomes small. In the limit of pi ’ 0 (but keeping the

expected number in the bin μi = N pi constant), the binomial distribution

becomes a Poisson distribution:

μni e’μ

P (ni |μi ) = . (19.42)

ni !

The expected value μi in bin (xi , xi+1 ) equals

xi+1

μi = N ρ(x) dx. (19.43)

xi

The observed number can be written as

ni = μi + µi , (19.44)

with µi a random variable with average zero and variance μi . Thus the

observed number ni has a standard deviation approximately equal to the

square root of ni itself. When the number of samples in each bin is large,

the histogram itself is a good representation of the required distribution

function.

When the available number of samples is small, it is not possible to con-

struct a histogram that has both a high resolution and a high precision.

Can we construct a smooth distribution function that is compatible with

the sampled data? There is one di¬culty: the number in each bin is not a

sample of the distribution function at one point (say, the center of the bin),

but rather a sample of the integral of the distribution function over the bin.

This di¬culty is resolved if we start with the cumulative distribution Ni :

i

Ni = nj , (19.45)

j=1

which is a sample of (N —) the cumulative distribution function F (xi ):

x

F (x) = ρ(x) dx. (19.46)

’∞

So we ¬t a spline to the observed numbers Ni /N , knowing that the function

must start at 0 and end at 1. We also know that the derivative of the

538 Splines for everything

%

100

80

60

40

20

-3 -2 -1 0 1 2 3

Figure 19.5 Cumulative distribution of 50 samples taken from a normal distribu-

tion.

function, being a density, must always be non-negative. But we must also

know the standard deviations of Ni /N .

Since

i

Ni = E[Ni ] + µj , (19.47)

j=1

we can approximate by the ¬rst term in a Taylor expansion

i n

Ni 1 1

≈ F (xi )(1 + µj ’ µj = F (xi ) + ·i . (19.48)

N E[Ni ] E[N ]

j=1 j=1

The variance of the relative error ·i can be easily computed, using the facts

that the variance of µj equals μj (expectation of nj ) and that µi and µj are

uncorrelated for i = j. We ¬nd

1 Ni Ni

Fi (1 ’ Fi ) ≈ 1’

2

σi = . (19.49)

N N N

The standard deviations are zero at both ends, which √ therefore forced

are