Introduction to Statistical Methods

PH3010 Advanced Skills

RHUL Physics 2019/2020

Abstract

In this module several important areas of statistical data analysis are explored, in-

cluding parameter fitting and model testing with the method of least squares. This is

done through a series of exercises that include fitting a polynomial and using the resulting covariance matrix of estimated parameters for error propagation, as well as analysis of historical data sets from Galileo and Ptolemy. Using measures of goodness-of-fit from the least squares fits, the level of agreement between the data and competing hypotheses is assessed.

1. Introduction

This module will provide an introduction to parameter fitting and model testing with the method of least squares, which are simple but powerful tools in the analysis of experimental data. The basic problem is illustrated in Fig. 1. Suppose we have measurements of a quantity y, each of which is carried out corresponding to a fixed value of some other variable x, called the control variable. Roughly speaking the goal is to find a smooth curve that passes close to the data points, called “fitting” a curve to the data.

10

8

6

y

Figure 1: Typical ingredients for

4

the least-squares problem: measure-

ments y with error bars σ are car-

2

ried out at known values of a con-

trol variable x. The curve has been

fitted to the data.

00

2

4

6

8

10

x

In Sec. 2 we will formulate this problem more precisely and state how the method of least squares is used to find the fitted curve. Then in Secs. 3, 4, 5 and 6 we will examine in greater detail several aspects of this problem: determining the fitted curve itself, quantifying its statistical uncertainty, testing whether the hypothesized functional form of the curve provides an adequate description of the data, and extending the procedure to the case of correlated data. Finally we present a set of exercises including analysis of historical data from Galileo and Ptolemy.

2. The method of least squares

Referring back to Fig. 1 we will suppose here that the x values are known very precisely, but the measured y values are subject to random fluctuations. For any given measurement we do not know exactly how far or in what direction the observed value has departed from its “true”

value, but we will assume we have some measure σ of the typical size of the fluctuations.

Suppose we have N measurements each labeled with an index i, and thus the data consist of the set of values (xi, yi, σi) with i = 1, . . . , N . These are displayed in Fig. 1 with the σi values as error bars.

Looking at the data in Fig. 1, one can believe that if it were possible to carry out the measurements without the random errors, i.e., in the case where σ is zero or negligibly small, then all of the (x, y) points would lie on a smooth curve. That is, there could be some function f (x) that gives the “true” values, and the observed values differ from these as a result of the fluctuations.

More precisely we will treat the yi as independent random variables. If one were to repeat the measurement of yi many times under the same circumstances (i.e., at the same xi), the observed yi values would come out different each time, and would follow a distribution with a certain width. To characterize this width we give the standard deviation σi. Note that in 1

the real problem we only have one value yi for a given xi; here the repeated measurements used to interpret the meaning of the standard deviation are purely hypothetical. The further statement that the set of yi values are independent means that the measured value of any one of them has no influence on the fluctuations of any of the others.

If we only have the set of (x, y, σ) measurements, how can we figure out the curve f (x)?

Given the available data, there is no way to give a definitive answer to this question. But we can estimate the curve using statistical methods. The approach we will take here requires that we hypothesize a formula for f (x) that contains some adjustable parameters, i.e., constants whose values are not yet determined. We will usually use Greek letters for the parameters and may write them as a vector if we have more than one, e.g., f (x; θ) with θ = (θ1, . . . , θm).

For example, we could suppose that f (x) is a straight line, i.e.,

f (x; θ) = θ0 + θ1x ,

(1)

where here θ = (θ0, θ1). The problem of estimating the curve is now reduced to estimating the unknown parameters. In this example, each possible choice for the values will yield a possible line, and the question is how to decide which point in the parameter space is best.

Suppose we pick an arbitrary point (θ0, θ1), which corresponds to a given curve, as shown in Fig. 2. Our initial choice is evidently not very good, as the amount by which the data points miss the curve, yi − f(xi; θ), is quite large for most of the points.

10

8

6

y

Figure 2: The line corresponding to

4

an arbitrarily chosen point in the pa-

resid

si u

d a

u l

a

rameter space (θ0, θ1) together with

2

the data points. The thick vertical

bars show the residuals yi −f(xi; θ).

00

2

4

6

8

10

x

Notice that for some of the points in Fig. 2, the difference or “residual” yi − f(xi; θ) is large but also the error bar σi is large. For these points one expects a bigger fluctuation in the measured value. Some of the measurements have small error bars, and so these should lie closer to the hypothesized curve. So to characterize how far a measured yi is from the predicted value f (xi; θ), we can measure the residual yi − f(xi; θ) in units of σi, i.e., we use yi − f(xi; θ) ,

(2)

σi

sometimes called the normalized residual. The “best” curve should give small values for the normalized residuals, but these can be positive or negative for different data points, so we cannot simply minimize their sum. In the method of Least Squares (LS), we define the best parameter estimates to be those that minimize the sum of squares of the normalized residuals, a quantity called the chi-squared,

2

N (y

χ2(

X

i − f (xi; θ))2

θ) =

.

(3)

σ2

i=1

i

The values of the parameters that give the minimum chi-squared, called the estimators of the parameters, are written with hats, e.g., ˆ

θ0 and ˆ

θ1. Figure 3 shows the example of the fitted

line, i.e., the curve f (x; θ) evaluated with the LS estimators for θ.

10

8

6

y

4

2

data

Figure 3: Least-squares fit of a

fit result

straight line (see text).

00

2

4

6

8

10

x

One can derive the method of least squares from a more general starting point called the method of maximum likelihood in the special case where the data values yi are statistically independent and each follow a Gaussian distribution (see, e.g., Ref. [2]). But it is equally valid to simply use the minimum of the χ2(θ) to define the estimators and it is important to emphasise that it is not possible in general to construct estimators that are “best” in every relevant sense. In practice it turns out that LS estimators have desirable properties in several respects, as we will see in examining the following questions:

1. How do we find the estimators, i.e., how do we minimize χ2(θ)?

2. How do we quantify the statistical uncertainty in the estimated parameters that stems from the random fluctuations in the measurements, and how is this information used

in an analysis problem, e.g., using error propagation?

3. How do we assess whether the hypothesized functional form f (x; θ) adequately describes the data?

3. Finding the minimum of χ2(θ)

Depending on the problem there can be a variety of different methods employed to find the estimators, i.e., to minimize the quantity χ2(θ). We will illustrate using a simple example first in a case where we can write down the estimators in closed form in Sec. 3.1, and then by carrying out the minimization numerically in Sec. 3.2.

3.1. Finding LS estimators in closed form

The curve f (x; θ) that we are fitting to the data does not of course have to be a straight line; it can be any function. It turns out that the method of least squares is particularly easy for 3

problems where f (x; θ) is a linear function of the parameters (it does not matter whether it is a linear function of the control variable x). As an example one could have a polynomial of order M , i.e.,

M

f (x;

X

θ) =

θnxn .

(4)

n=0

For this type of problem it is possible to find the values of the parameters that minimize the chi-squared in closed form, and we will show this here for the special case of a first-order polynomial, i.e., a straight line. Here the the function χ2(θ) takes on the form

N (y

χ2(θ

X

i − θ0 − θ1xi)2

0, θ1) =

.

(5)

σ2

i=1

i

To find its minimum we set the derivatives of χ2 with respect to both parameters equal to zero,

∂χ2

N −2(y

= X

i − θ0 − θ1xi) = 0 ,

(6)

∂θ0

σ2

i=1

i

∂χ2

N −2x

= X

i(yi − θ0 − θ1xi) = 0 .

(7)

∂θ1

σ2

i=1

i

These two equations can be rewritten using matrix notation as



PN

1

PN

xi 



y







PN

i

i=1 σ2

i=1 σ2

θ0

i=1

i

i

σ2i

















= 



,

(8)











x



PN

x

x2

PN

i

PN

i

iyi

i

θ

=1 σ2

i=1 σ2

1

i=1 σ2

i

i

i

which has the general form



a b   θ 





0

e



 

 = 



.

(9)

c d

θ1

f

By comparing Eqs. (8) and (9) we can read off the values of a, b, c, d, e and f . Then using the fact that a 2 × 2 matrix A and its inverse A−1 can be written

a b !

1

d

−b !

A =

,

A−1 =

,

(10)

c d

ad − bc

−c

a

and applying A−1 to both sides of Eq. (9), we can find solutions

ˆ

de − bf

θ0 =

,

(11)

ad − bc

âf − ec

θ1 =

.

(12)

ad − bc

Here we write the solutions with hats because these are the estimators for the parameters θ0

and θ1.

4

We can extract a few lessons from this example. First, we see that the solutions can be obtained by solving a system of linear equations using standard matrix techniques. Fur-thermore, the data values yi enter into the estimators through the quantities e and f, and therefore from Eqs. (11) and (12) we see that the estimators ˆ

θ0 and ˆ

θ1 are linear functions

of the yi. This remains true for an arbitrary number of parameters, as long as the function that we are fitting, f (x; θ), is linear in the parameters (see, e.g., Ref. [2]). Finally we see that although we can write down expressions for the solutions in closed form, they are somewhat involved even for the simple case of two parameters. If the number of parameters is higher then this approach may become impractical, and in any case if f (x; θ) is nonlinear in any of the parameters then we must use other methods to find the minimum of χ2(θ).

3.2. Numerical minimization of χ2(θ)

In most problems of practical interest it is too difficult or impossible to find formulas for the estimators in closed form as done in the previous section. Rather, we must find the minimum of χ2(θ) (called the objective function) numerically. The study of algorithms used to carry out this task constitutes an entire branch of applied mathematics called mathematical optimization.

Here we will not go into the details of these algorithms other than to mention some of their general principles. The basic idea is to evaluate the objective function χ2(θ) at different points in the parameter space and in this way to try to figure out where it has its minimum.

Some algorithms make use of the derivatives of the objective function with respect to the parameters if these are known. One must start at a certain point in the parameter space, say, θstart. Sometimes default values can be used for this; in other problems one may need to have a rough preliminary estimate.

For example, starting with some initial guess for the parameters, hold all but one of them constant, say, θ0, and determine the value of θ0 that minimizes the χ2. This can be done by finding the solution to ∂χ2/∂θ0 = 0 using, e.g., the Newton–Raphson method. Then hold all parameters constant at the newly found point except the next one, say, θ1, and determine its value so as to minimize the χ2. After cycling through all of the parameters one will have moved to a point closer to the minimum. Repeat the procedure as many times as necessary until the change in the parameters from one iteration to the next drops below some threshold, or alternatively, until the change in the χ2 between iterations falls below a specified value.

In fact the most successful algorithms do not minimize by varying separately individual parameters but rather by moving along certain well-chosen directions in the parameter space, e.g., in the direction where the objective function is decreasing most rapidly (called the method of steepest descent) or using so-called conjugate gradient directions. More details on these types of algorithms can be found in Chapter 10 of Ref. [1].

Here we will show a simple example of numerical minimization using routines in Python.

Let us suppose we have the data points shown above in Fig. 1 and we want to fit a straight line. A simple but relatively flexible tool in Python is the routine curve fit, which is part of the package scipy.optimize. To use this we first need to import the relevant packages: import numpy as np

from scipy.optimize import curve_fit

We need to define the function that we will be fitting to the data. For this example suppose it is a straight line:

5

def func(x, *theta):

theta0, theta1 = theta

return theta0 + theta1*x

Then we need to have the data values (xi, yi, σi) in the form of NumPy arrays. In practice we would probably read the data values in from a file, but for this example let us suppose we have the arrays defined by

x

= np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0])

y

= np.array([2.7, 3.9, 5.5, 5.8, 6.5, 6.3, 7.7, 8.5, 8.7])

sig = np.array([0.3, 0.5, 0.7, 0.6, 0.4, 0.3, 0.7, 0.8, 0.5])

To carry out the least-squares fit we usually want to supply initial values for the parameters.

If we have reasonable guesses then these can be used; if they are not supplied then the program will take the initial values equal to 1.0. Looking at the graph we see that the slope and y intercept are both of order unity, so we define

p0 = np.array([1.0, 1.0])

In practice the result should be insensitive to the initial values, but for some special problems one may need to choose this point with some care. To actually fit the curve we use

thetaHat, cov = curve_fit(func, x, y, p0, sig, absolute_sigma=True)

Notice that it is necessary here to include the argument absolute sigma=True in order for the statistical errors of the estimated parameters to have the interpretation that we will require. This results in ˆ

θ0 = 2.26 ± 0.29 and ˆ

θ1 = 0.741 ± 0.057, as represented by the fitted

curve shown earlier in Fig. 3. Here the indicated statistical errors correspond to the standard deviations of the estimators, as described in the next section.

4. Statistical errors of the fitted parameters

In this section we will examine the statistical uncertainties in the estimated parameters, which stem from the random fluctuations in the original measurements. First in Sec. 4.1 we will discuss in some detail how to interpret the statistical errors and how to find them using the Monte Carlo method. In Sec. 4.2 we will show a much simpler way to estimate the errors using derivatives of the function χ2(θ), and in Sec. 4.3 we show how to use the statistical uncertainties with error propagation.

4.1. Meaning of statistical errors and estimation by Monte Carlo method

First let us consider carefully the statement that the estimators ˆ

θ (e.g., ˆ

θ0 and ˆ

θ1 in the

example above) have statistical errors. What this means is that if we were to repeat the entire experiment with a new, statistically independent set of measurements yi with i = 1, . . . , N , then the fluctuations in the data would in general be different from before and therefore the numerical values of the estimates ˆ

θ would also differ. Equivalently we can say that the

estimators ˆ

θ are functions of random variables y1, . . . , yN , and they are therefore random

variables themselves.

Although we may only carry out the real experiment once, we can simulate it many

times with the Monte Carlo method. Figure 4 shows the outcomes of many independent

6

simulations of the same set of measurements of the straight-line fit described above, each of which is characterized by two fitted parameters, ˆ

θ0 and ˆ

θ1.

1̂0

0̂8

̂ θ1

Figure 4: Estimates ˆ

θ0 and ˆ

θ1 from

0̂6

repeated simulations of the straight-

line fit, each based on a statistically

independent set of measurements.

1

2

3

̂θ0

Figures 5 show the projections of the points from the scatter plot onto the ˆ

θ0 and ˆ

θ1 axes,

which give histograms of the corresponding quantities. The widths, or more precisely, the standard deviations of these distributions are used to quantify the statistical errors in the parameters, which we can write as σând σ

(note the hats).

θ

ˆ

0

θ1

(a)

1̂5

8

(b)

6

1̂0

4

0̂5

̂

0̂0

0

1

2

3

0.6

0.8

1.0

̂θ

̂

0

θ1

Figure 5: Projections of the scatter plot from Fig. 4 onto (a) the ˆ

θ0 and (b) the ˆ

θ1 axes.

Returning to the scatter plot in Fig. 4, we see another interesting feature, namely, that the cloud of points is tilted. This shows that if the value of, say, ˆ

θ0 comes out higher than

average, then there is a higher probability for the value of ˆ

θ1 to be lower than its average.

That is, the quantities ˆ

θ0 and ˆ

θ1 are correlated, in this case the correlation being negative.1

Before we quantify the statistical errors in the estimators we need to review briefly how to characterise the typical amount of fluctuation in random variables. Recall that to do this for a single random variable, say, u, one defines the variance V [u] as the mean of its square minus the square of its mean, i.e.,

V [u] = hu2i − hui2 .

(13)

1 Often one hears statements such as “the parameters are correlated”, but this is not strictly correct. Here we are treating the parameters θ0 and θ1 as unknown constants, not as random variables. It is the estimators, ˆ

θ0 and ˆ

θ1, which are functions of the random data, that vary upon repetition of the experiment, and therefore we should say that the estimators of the parameters are correlated.

7

Often we use the notation V [u] = σ2u, and the square root of the variance, σu, is the standard deviation. In a similar way, to characterise the amount of variation in a pair of random variables, say, u and v, we define the covariance as the mean of their product minus the product of the means,

cov[u, v] = huvi − huihvi .

(14)

If u and v have some units such as kg or cm, then the covariance has the units given by their product. It is often convenient to use the dimensionless correlation coefficient ρ, defined as the covariance of the two variables divided by the product of their standard deviations, cov[u, v]

ρ =

.

(15)

σuσv

One can show that the correlation coefficient for any pair of variables lies in the range

−1 ≤ ρ ≤ 1. The extreme values of ρ = ±1 corresponding to a 100% positive (+1) or

negative (−1) correlation (i.e., perfect anticorrelation). The corresponding scatter plot is then infinitesimally thin and tilted up (ρ = +1) or down (ρ = −1). If ρ = 0 then the

variables are uncorrelated, and their scatter plot is not tilted.

If one has more than two random variables, e.g., a set of m estimators ˆ

θ = (ˆ

θ1, . . . , ˆ

θm),

then there is a covariance for each pair, and these numbers can be arranged in a covariance matrix,

Uij = cov[ˆ

θi, ˆ

θj] .

(16)

From the definition of covariance it is easy to check that the matrix is square and symmetric.

The diagonal elements simply give the variances

Uii = cov[ˆ

θi, ˆ

θi] = σ2ˆ ,

(17)

θi

and so by taking their square roots one finds the standard deviations σˆ . The set of corre-

θi

lation coefficients also form a matrix ρij, given by

U

ρ

ij

ij =

.

(18)

σˆ σ

θ

î

θj

It is easy to verify that the diagonal elements are ρii = 1 for all i, which is to say that a random variable is always 100% positively correlated with itself.

From the scatter plot of the simulated estimates ˆ

θ0 and ˆ

θ1 we see that they have a

negative correlation. Using standard formulas or NumPy routines to estimate covariance from a sample of observed values (see, e.g., Refs. [2, 3]), we find

σˆ

= 0.29 ,

(19)

θ0

σˆ

= 0.057 ,

(20)

θ1

cov[ˆ

θ0, ˆ

θ1] = −0.0142 ,

(21)

ρ = −0.86 .

(22)

8

4.2. Statistical errors from derivatives of χ2(θ)

The Monte Carlo technique described in the previous section provides estimates of the covariance of the estimators but it requires a certain amount of work to set up the simulations and work out the required numbers. The covariance matrix can also be determined directly from the single real data set, with no need of simulated data, by using the second derivatives of χ2(θ) with respect to the parameters.

Derivation of the required formulas is nontrivial and goes beyond the scope of this module.

We claim here without proof that the inverse of the covariance matrix for a set of estimators Uij = cov[ˆ

θi, ˆ

θj] can be estimated by

1 ∂2χ2 !

U −1

ij

=

.

(23)

2

∂θi∂θj θ=ˆθ

In general this formula is an approximation but it is exact or at least close to being so in most cases of practical interest. When using the Python routines described above, the curve fit function finds the minimum of χ2(θ) numerically and then by taking small steps about the minimum, it finds the matrix of second derivatives to get U −1. It then inverts this and returns the covariance matrix U .

For the example of the straight line shown above, curve fit returns the covariance matrix



0.08537

−0.01438 

U = 



,

(24)

−0.01438 0.003275

where the first row/column corresponds to the estimator ˆ

θ0 and the second to ˆ

θ1. By taking

the square root of the diagonal elements one finds the standard deviations σˆ = 0.29 and θ0

σˆ = 0.057, and combining with the covariance one finds the correlation coefficient ρ

θ

01 =

1

−0.86. All of these numbers agree to good level of precision with the values found using the Monte Carlo calculation in Sec. 4.1.

4.3. Using the covariance of estimators with error propagation

Suppose we have carried out a least-squares fit and found a set of estimators ˆ

θ = (ˆ

θ1, . . . , ˆ

θm)

and their covariance matrix Uij = cov[ˆ

θi, ˆ

θj]. Very often we then form some function u of the

estimators as the quantity of interest, e.g., u(ˆ

θ) = (ˆ

θ ˆ

1θ2 − 2 ˆ

θ3)2. Because the estimators ˆ

θi

are themselves random variables, a function of them is also a random variable with a certain variance. To characterise the statistical uncertainty in the function, we need to “propagate”

the covariance of the ˆ

θi through to u. An analogous procedure is carried out with more

than one function, e.g., u and v, whereby one finds their variances as well as the covariance cov[u, v].

One way to carry out this procedure would be to simulate the entire experiment many

times, just as in Sec. 4.1. Each time we determine the estimates ˆ

θ, with these evaluate the

function u(ˆ

θ), and record its value, e.g., in a histogram. By repeating a sufficient number of

times one obtains the distribution of u and thus its standard deviation can be found.

Although the Monte Carlo procedure is robust in the sense that it is essentially guaranteed to produce a meaningful result, it may require some effort to set up. So it is useful to have a simpler approximate method, called (linear) error propagation. This is based on a linear 9

approximation to the function in the neighborhood of the estimates actually found. For a pair of functions u(ˆ

θ) and v(ˆ

θ) one finds for the covariance (see, e.g, Ref. [2]),

m

∂u ∂v

cov[u, v] ≈ X

Uij ,

(25)

i,j

∂ ˆ

θ ∂ ˆ

θ

=1

i

j

where Uij = cov[ˆ

θi, ˆ

θj] is the covariance matrix of the estimators of the parameters and the

derivatives are evaluated at the estimated values of the parameters. In the case where there is only one function, say, u, its variance is given by σ2u = cov[u, u], i.e.,

m

∂u ∂u

σ2

X

u ≈

Uij .

(26)

i,j

∂ ˆ

θ ∂ ˆ

θ

=1

i

j

In the special case where the covariance may be diagonal, i.e., Uij = δijσ2 , then one can ˆ

θi

carry out one of the two sums to find

m ∂u !2

σ2

X

u ≈

σ2ˆ .

(27)

θi

i

∂ ˆ

θ

=1

i

It should be emphasised, however, that estimators for parameters are correlated in general, i.e., the covariance matrix usually has nonzero off-diagonal elements, and thus one must use the full covariance in propagating the errors.

The approximation on which Eqs. (25), (26) and (27) are based is that the functions

should be linear in the neighborhood of the estimated values, where the size of the relevant neighborhood is determined by the standard deviations of the corresponding estimators.

Suppose, for example, we have a function of the form u = ˆ

θ1/ˆ

θ2, and that ˆ

θ2 = 10 ± 1. The

function is (exactly) linear in ˆ

θ1, and although it is nonlinear in ˆ

θ2, this is approximately linear

if we only consider the small region around 10 ± 1. So in this case linear error propagation should be a good approximation. If, on the other hand, we have ˆ

θ2 = 10±10, then the function

is very nonlinear if we consider this larger range. In this case linear error propagation can give wildly wrong results. The Monte Carlo approach, however, could be used in such a case to infer the correct distribution of the function u.

5. Goodness-of-fit from the minimum of χ2(θ)

In the example from Sec. 2 we hypothesized that the data shown in Fig. 3 could be described by a straight line. Having carried out the fit we then found the values of the parameters for the slope and y-intercept. But how do we know that this hypothesis is valid? Perhaps we should have used a curve with more adjustable parameters, such as the third-order polynomial fitted to the same data in Fig. 1. Or perhaps the data come out as in Fig. 6 below; in this case the hypothesis of a straight line is clearly not very good. Beyond a visual comparison of the fitted curve with the data, how can we quantify the “goodness of fit”?

Before we go any further we should head off a common misunderstanding. One might

think that if the hypothesized function is bad, then this should be reflected in large values of the statistical errors for the fitted parameters. This is not true. The level of agreement between data and hypothesis and the size of statistical errors are in general very separate questions, and in the particular example shown here they are completely unrelated. In fact, 10

10

8

6

y

4

Figure 6: Fit resulting in a poor

2

level of agreement between the hy-

data

pothesized fit function and the data.

fit result

00

2

4

6

8

10

x

to obtain the data points shown in Fig. 6, the points from the earlier example (Fig. 3) were simply shifted up or down; their error bars were left unchanged. And as a result, the standard deviations of the estimators of the parameters shown in Fig. 6 are exactly the same as before.

The reason for this perhaps surprising behaviour lies in the meaning of the statistical errors. What they reflect is the level with which the parameter estimates would fluctuate upon repetition of the experiment with a new set of measured yi. It could be that each repetition results in poor agreement, as found in Fig. 6, but that the fitted parameters do not vary much for each repetition, so the standard deviations of the estimators may be small.

What we would like to quantify, however, is something different, namely the level to which a straight line provides an adequate description of the data.

One way to address this question is to look at the value of the χ2 at its minimum, i.e., at θ = ˆ

θ,

N (y

χ2

= X

i − f (xi; ˆ

θ))2 .

(28)

min

σ2

i=1

i

From the form of this equation we can immediately see that if the fitted function is in good agreement with the data, then the contribution of each term in the sum is small and therefore χ2

should be small, so this statistic can be used to reflect the goodness-of-fit. But we should min

not be surprised to see some nonzero value even if the fit function is correct, because of the statistical fluctuations in the data. To make a quantitative statement about the goodness-of-fit we need to know how large of a χ2

to expect in the case where the hypothesized fit

min

function is correct.

Recall that the standard deviation σi reflects the typical (more precisely, the root mean square) level of departure of the measurement yi from its mean. That is to say, if the fitted value f (xi; ˆ

θ) were equal to the true mean of yi, then one would expect |yi − f(xi; ˆθ)| to be

similar to σi and therefore each normalized residual (cf. Eq. (2)) should be of order unity, yi − f (xi; ˆ

θ)

∼ 1 .

(29)

σ

i

Since there are N terms in the sum, one would naively expect χ2

to be roughly equal to

min

N . In fact this cannot be quite right, since we have adjusted some number of parameters, so the fitted values f (xi; ˆ

θ) get pulled even closer to yi than would be obtained from the “true”

function. For example, if one had N data points and m = N parameters, then one would be 11

able to adjust their values so that the curve goes exactly through each of the points and one would find χ2

= 0.

min

The minimized value χ2

is of course a function of the data and as such it is itself a

min

random variable. If one were to repeat the experiment many times, its value would follow a certain distribution. One can show that if the hypothesised form of the function f (x; θ) is correct, and if the yi are Gaussian distributed, and if certain additional conditions are satisfied that usually hold in practice, then the statistic χ2

will follow a probability distribution

min

function (pdf) called, not surprisingly, the chi-squared distribution. Calling the variable χ2min here z to simplify the notation, the chi-squared pdf is given by

1

fχ2(z; n) =

zn/2−1e−z/2 .

(30)

2n/2Γ(n/2)

Here the parameter n is called the “number of degrees of freedom” and takes on positive integer values, n = 1, 2, . . .. If there are N measurements y1, . . . , yN and m fitted parameters θ = (θ1, . . . , θm), then the number of degrees of freedom is given by n = N − m.

The distribution of χ2

for the example of the straight-line fit from repeated simulations

min

of the measurement is shown as a histogram in Fig. 7. The fit had N = 9 measured y values and m = 2 fitted parameters, so the relevant number of degrees of freedom is n = 7. Also shown on the plot as a solid curve is the chi-squared pdf for 7 degrees of freedom.

0.15

0.10

Figure 7: Distribution of the value

0.05

of χ2

from Monte Carlo simula-

min

tions of the measurement and fit of a

straight line. The solid curve shows

the chi-squared pdf for n = 7 de-

0.00

grees of freedom.

0

5

10

15

20

25

χ2min

The mean and standard deviation of a chi-squared distributed variable for n degrees of freedom are given by

hzi = n ,

(31)

√

σz =

2n .

(32)

Since the mean is equal to the number of degrees of freedom, we can immediately get a good idea of the goodness-of-fit simply by comparing the χ2

obtained to n = N

min

− m, the number

of measurements minus number of fitted parameters. Sometimes this is given as a ratio, χ2 /n, called the chi-squared per degree of freedom. If this ratio comes out close to unity, min

the fit is “good”; if it is very large, the fit is bad, and if it is much less than one then it means that the data points have come out closer to the hypothesised curve than what one would expect, given the random fluctuations that should normally be present.

12

For the fit shown in Fig. 6, there are 9 measurements and 2 fitted parameters, so if the straight-line hypothesis is true, one would expect a value of χ2

of around 7. In the fit

min

shown, a value χ2

= 20.9 was found, or a chi-squared per degree of freedom of almost 3.

min

This can be contrasted with the fit shown in Fig. 3, which has χ2

= 8.2, for the same

min

number of degrees of freedom.

The fact that the value of χ2

in Fig. 6 came out quite high does not prove that the

min

straight-line hypothesis is wrong. It can happen, as a result of of the random fluctuations in the data, that the value of χ2

comes out sometimes high and sometimes low. After all it

min

is a random variable, as illustrated by its distribution shown in Fig. 7. We can take a step further in quantifying the goodness-of-fit by giving the probability, under assumption of the hypothesis used (here a straight line), to find a value of χ2

as high as the one found or

min

higher. This can be found by integrating the chi-squared pdf from the value of χ2

observed

min

to infinity,

Z ∞

p =

fχ2(z; n) dz ,

(33)

χ2min

where fχ2(z; n) is the chi-squared pdf for n degrees of freedom from Eq. (30). The probability p obtained from Eq. (33) is called the p-value of the hypothesis. If it comes out very small, it means that there is only a very small probability, if the hypothesis is true, to find data at least as incompatible with it or more so. We may therefore be tempted to think of the p-value as something like the probability that the hypothesis is true, but this is not a correct statement. There is an important distinction between these quantities, which we will not consider in detail here (further discussion can be found, e.g., in Ref. [2]). For now we will simply emphasize that the p-value is a measure of compatibility between the data and a given hypothesis, and if it comes out very small then the hypothesis is disfavoured.

For the straight-line fit shown in Fig. 3, we had χ2

= 8.2 and 7 degrees of free-

min

dom. Using these values in Eq. (33) and carrying out the integral with the Python routine scipy.stats.chi2.sf gives a p-value of 0.31. Thus if the straight-line hypothesis is correct, there is a substantial probability to get an even higher value of χ2 , and so we would not min

reject the hypothesis in this case. But in the fit in Fig. 6 we found χ2

= 20.9, corresponding

min

to a p-value of 0.0039. It is still possible that the straight-line hypothesis is correct, but if it is, there is only a 4 in a thousand chance to see a χ2

at least as high as the one we found.

min

If we find a very low p-value (or equivalently, χ2

substantially greater than the number

min

of degrees of freedom), then we may wish to try another hypothesis. Using the same data as in Fig. 6 we can try, for example, a second-order polynomial,

f (x; θ) = θ0 + θ1x + θ2x2 .

(34)

The fit using this function is shown in Fig. 8. It gives χ2

= 3.5 for 6 degrees of freedom,

min

corresponding to a p-value of 0.75, and thus a much better level of agreement between the data and hypothesis.

6. Least squares with correlated measurements

In the method of least squares as described above, the quantity χ2(θ), whose minimum defines the estimators ˆ

θ, is given by

13

10

8

6

y

4

Figure 8: Fit to the same data val-

2

ues as in Fig. 6 but here using a fit

test data

function with three adjustable pa-

fit result

rameters.

00

2

4

6

8

10

x

N (y

χ2(

X

i − f (xi; θ))2

θ) =

.

(35)

σ2

i=1

i

As mentioned in Sec. 2, one can show that this recipe follows from a more general principle, namely, the method of maximum likelihood, for the special case where the measured quantities yi are independent and Gaussian distributed. In practice the measured quantities we deal with very often satisfy these two criteria, and this is an important reason why the method of least squares is so widely used.

It is not uncommon, however, that measured values are Gaussian but they are correlated, and hence not independent. That is, each measurement does not simply have a variance V [yi], but the set of N measurements is described by an N × N covariance matrix Vij = cov[yi, yj].

The elements of this matrix are related to the correlation coefficients through

Vij = σiσjρij .

(36)

Equation (36) includes the special case of i = j, i.e., the diagonal elements have a correlation coeffficient ρii = 1 and the corresponding elements Vii = σ2i are the variances. In such a problem, one can generalise the method of least squares by replacing the sum of squares given above by

N

χ2(

X

θ) =

(yi − f(xi; θ))V −1

ij (yj − f (xj ; θ)) ,

(37)

i,j=1

where V −1 is the inverse of the covariance matrix. As in the uncorrelated case, this version of χ2(θ) can be derived from the method of maximum likelihood where the data follow a

multivariate Gaussian distribution (see, e.g., [2]).

If we use Eq. (37) for the case where the measurements are in fact uncorrelated, then the covariance matrix is diagonal, i.e., Vij = δijσ2i, and thus its inverse is simply V −1

ij

= δij/σ2i.

Substituting this expression into Eq. (37) allows one to carry out one of the two sums, and the formula for χ2(θ) becomes again a sum of squares as in Eq. (35).

14

Exercises

Here we give several exercises that illustrate the methods described above. Exercise 1 is similar to the straight-line fit shown above and requires an extended error analysis. Exercises 2

and 3 use historical data on projectile motion from Galileo and on refraction of light from Ptolemy. The goodness-of-fit is analyzed to choose between competing hypotheses. Exercise 4

explores algorithms that can be used to find the minimum of χ2(θ) numerically.

Exercise 1: Polynomial fit and error analysis

Consider the following set of (x, y, σ) data points:

x

= np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0])

y

= np.array([2.7, 3.9, 5.5, 5.8, 6.5, 6.3, 7.7, 8.5, 8.7])

sig = np.array([0.3, 0.5, 0.7, 0.6, 0.4, 0.3, 0.7, 0.8, 0.5])

1(a): Using these data carry out the least-squares fit of an M th order polynomial, i.e., with M + 1 adjustable parameters, for M = 1, 2, 3.

1(b): For each fit, use the error propagation formula Eq. (26) to find the standard deviation of the fitted function σf as a function of x. Note that to do this you will need to compute the derivatives of f (x; ˆ

θ) with respect to the components of ˆ

θ. Display the fitted curve plus-

or-minus one standard deviation as a shaded band, and extend the x axis to at least 20.

(The shaded band can be made with the function matplotlib.fill between.) The result

for M = 1 is shown in Fig. 9. Note how the size of the error band increases when one goes to x values outside the region where data are available; investigate how this behaviour changes as the order of the polynomial is increased.

16

12

y

8

4

Figure 9: Result of least-squares fit

test data

with the one-sigma error shown as a

fit result

shaded band.

00

5

10

15

20

x

1(c): Consider the fit with M = 3 and define the difference

∆ab(ˆθ) = f(a; ˆθ) − f(b; ˆθ) .

(38)

Using error propagation, find the standard deviation of ∆ab(ˆθ) for a = 5 and b = 6, 10, 20.

Compare these values you find with the standard deviation of f that you plotted for this fit as a shaded band evaluated at both a and b.

15

Exercise 2: Analysis of Galileo’s ball and ramp data

In this exercise we consider an experiment performed by Galileo using a ball and an inclined ramp as shown in Fig. 10. The ball starts at a height h above the edge of the ramp, and its trajectory is forced to be horizontal before it falls over the edge. The horizontal distance d from the edge to the point of impact is measured for different values of h. Five data points obtained by Galileo in 1608 are shown in Table 1 (from Ref. [4]).

h

Figure 10:

The configuration

Galileo’s ball and ramp experiment.

d

Table 1: Galileo’s data on horizontal distance before impact d for five values of the starting height h. The units are punti (points); one punto is slightly less than 1 mm.

h

d

1000

1500

828

1340

800

1328

600

1172

300

800

We will assume the heights h are known with negligible error, and that the horizontal

distances d have uncertainties of σ = 15 punti (points); one punto is slightly less than 1 mm.

It is not actually known what the measurement uncertainties were, but 1–2% is plausible. In addition, we know that if h = 0, then the horizontal distance d will be zero, i.e. if the ball is started at the very edge of the ramp, it will fall straight down to the floor.

We could proceed by applying the known laws of mechanics to the system and deriving

the relationship between d and h. For purposes of this exercise, however, we will pretend we do not yet know Newton’s laws (as was the case for Galileo), and we will simply try different hypotheses and compare their predictions with the data.

Consider the following three hypotheses for the functional relationship between d and h: a linear relation with a single free parameter α,

d = αh ,

(39)

a quadratic relation with two parameters, α and β,

d = αh + βh2,

(40)

16

and something which is a nonlinear function of the parameters,

d = αhβ .

(41)

The goals of this exercise are to investigate these three hypotheses. More specifically, you should write a Python program to do the following:

2(a) Carry out a least-squares fit of the parameters for each hypothesis. Find the estimated parameter values, their standard deviations and if relevant their covariance and correlation coefficient.

2(b) Produce plots of the fitted curves along with the data points and their errors.

2(c) Quantify the level of agreement between the data and each hypothesis with both a

p-value and the chi-squared per degree of freedom. Discuss which of the hypotheses are disfavoured and which are preferred.

2(d) Use Newton’s Laws to derive the relation between d and h and predict the values of the parameters. Discuss how you would compare this prediction to the fit results obtained above.

What parameters of the system can be related using your theoretical prediction to the fitted parameters? What important experimental uncertainties could render such a comparison

difficult?

Bonus question: From the relation between d and h that you derive in 2(d), one can predict the height of the table in punti. Find the height and using error propagation determine its statistical error from those of the fitted parameters. Compare your result to the value that Galileo recorded in his notebook. (See Appendix A and look where it says ”altezza della tavola”.)

Exercise 3: Analysis of refraction data from Ptolemy

The astronomer Claudius Ptolemy performed experiments on the refraction of light using a circular copper disc submerged to its centre in water, as illustrated in Fig. 11. Angles of refraction θr for 8 values of the angle of incidence θi obtained by Ptolemy around 140 a.d.

are shown in Table 2 (from Ref. [5]).

θ i

incident ray

air

water

Figure 11: The apparatus used by

Ptolemy to investigate the refraction

copper disc

of light.

refracted ray

θ r

For purposes of this exercise we will take the angles of incidence to be known with negligible error and treat the angles of reflection as independent Gaussian-distributed measurements 17

Table 2: Angles of incidence and refraction (in degrees).

θi

θr

10

8

20

15 12

30

22 12

40

29

50

35

60

40 12

70

45 12

80

50

with standard deviations of σ = 1 ◦. (This is a reasonable guess given that the angles are 2

reported to the nearest half degree. Note that we can absorb an error in θi into an effective error in θr.)

Until the discovery of the correct law of refraction (see below), a commonly used hypothesis was

θr = αθi,

(42)

although it is reported that Ptolemy preferred the formula

θr = αθi − βθ2.

(43)

i

3(a): Find the least-squares estimates of the parameters for both hypotheses and determine the minimized χ2. Comment on the goodness-of-fit for both hypotheses. Is it plausible that all of the data values are based on actual measurements?2

The law of refraction discovered by the Persian mathematician and physicist Ibn Sahl in the 10th century and rediscovered by others including Snell in 1621 is

sin θ

i

θr = sin−1

,

(44)

r

where r = nr/ni is the ratio of indices of refraction of the two media.

3(b): Determine the least-squares estimate for r and find value of the minimized χ2. Comment on the validity of the Gaussian assumption for θ

◦

r with σ = 1 .

2

For all of the fits, report the parameter estimates, their standard deviations, and where relevant the covariance and correlation coefficient. Make plots of the fitted curves and the data.

Exercise 4: Algorithms for numerical minimization

In most applications we use a specialised software package such as scipy.optimize to minimize an objective function like χ2(θ). It is nevertheless an interesting exercise to create such 2 See R. Feynman, R. Leighton and M. Sands, The Feynman Lectures on Physics, Vol. I, Addison-Wesley, Menlo Park, 1963, Section 26-2.

18

a routine yourself. The goal here is to implement a minimization algorithm in Python and to compare its results to those obtained from curve fit.

Because the mathematics of function minimization can be treated separately from the

method of least squares, here we will change notation and suppose that we have an objective function f (x) of an n-dimensional argument x = (x1, . . . , xn). The goal is to find the point in x-space that gives the minimum of f (x). We will assume in doing so that the analyst supplies a starting point xs, which can be based on an initial guess.

All methods that achieve this goal require that one be able to evaluate the objective

function f at an arbitrary point x. If, in addition, one is able to evaluate the derivatives of f with respect to the components of x then algorithms can be used that will find the minimum much faster. In some cases one may know the derivatives in closed form; in other situations one may have to estimate them numerically. For the functions we have considered in the examples above it is relatively simple to calculate the first and second derivatives of χ2(θ) with respect to the parameters, and so for the method we implement here we will suppose that these are available.

One simple strategy to minimize f (x) is to consider all of the components of x as fixed except one. Starting from the point x0, we vary only this chosen component to minimize f .

We then repeat this procedure for all of the components of x in turn.

Consider first the function f as depending only on one non-constant component, which

here for simplicity we will call x, without a subscript. Suppose that at a particular step of the algorithm the value of x is x0, i.e., here the index refers to a step of the algorithm, not to the component. We can expand f (x) in a Taylor series to second order about its initial value x0, i.e.,

1

f (x) ≈ f(x0) + f′(x0)(x − x0) +

f ′′(x0)(x − x0)2 .

(45)

2!

To minimize this function with respect to x, we set its derivative equal to zero,

f ′(x) ≈ f′(x0) + f′′(x0)(x − x0) = 0 .

(46)

Solving for x gives

f ′(x0)

x = x0 −

.

(47)

f ′′(x0)

This x is used to updated the corresponding component of x, which in turn gives an updated estimate of the position of the minimum of f (x). Having updated one component, we then cycle through all of the remaining ones and minimize with respect to each while holding the rest constant.

After carrying out minimizations with respect to all of the components we can evaluate the objective function and see how much it changed from its value at the previous full step.

The procedure is repeated until the change in f (x) from one iteration to the next falls below a desired threshold, which determines the final numerical accuracy.

One can easily generalise this algorithm to minimize not with respect to the individual coordinates, but rather along an arbitrary direction in x-space. If we have available the derivatives of f , i.e,

19

∂f

∂f

∇f =

, . . . ,

,

(48)

∂x1

∂xn

then we can generalise the procedure described above to minimize f along the direction of

−∇f, i.e., the direction in which the function is decreasing most rapidly. This is called the method of steepest descent.

Although the method of steepest descent should eventually find the minimum of the

objective function, in many cases it turns out to be faster to choose what are called conjugate gradient directions; more details can be found in Chapter 10 of Ref. [1].

4(a): The goal of this exercise is to create a Python program that will find the minimum of the objective function χ2(θ) numerically. Implement the algorithm with minimization along the coordinate directions as described above and test it by fitting a polynomial to the data shown in Sec. 6. Use different orders for the polynomial and thus different numbers of fitted parameters. Compare the results you obtain with those from the routine curve_fit and

comment on any differences.

4(b): Extend exercise 4(a) by using the method of steepest descent.

References

[1] S. Brandt, Statistical and Computational Methods in Data Analysis, Springer, New

York (1997).

[2] G. Cowan, Statistical Data Analysis, Oxford University Press, 1998.

[3] RHUL Physics Department, Mathematical Formula Sheet for Physics (2017); available

on moodle.

[4] Stillman Drake and James Maclachlan, Galileo’s discovery of the parabolic trajectory, Scientific American 232 (March 1975) 102; Stillman Drake, Galileo at Work, University

of Chicago Press, Chicago (1978).

[5] Olaf Pedersen and Mogens Pihl, Early Physics and Astronomy: A Historical Introduc-

tion, MacDonald and Janes, London, 1974.

[6] Galileo Galilei,

Manuscript f.116,

Biblioteca Nazionale Centrale di Firenze,

www.bncf.firenze.sbn.it.

20

A. Galileo’s notes on the ball and ramp experiment

Manuscript f.116 from Galileo’s notebooks showing the ball and ramp experiment (from [6]).

21

B. Python code

Program simpleFit.py for polynomial fits.

1

#

simpleFit.py

2

#

G. Cowan / RHUL Physics / October 2017

3

#

Simple program to illustrate least-squares fitting with curve_fit

4

5

import matplotlib

6

import matplotlib.pyplot as plt

7

import numpy as np

8

from scipy.optimize import curve_fit

9

10

# define fit function

11

def func(x, *theta):

12

theta0, theta1 = theta

13

return theta0 + theta1*x

14

15

# set data values

16

x

= np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0])

17

y

= np.array([2.7, 3.9, 5.5, 5.8, 6.5, 6.3, 7.7, 8.5, 8.7])

18

sig = np.array([0.3, 0.5, 0.7, 0.6, 0.4, 0.3, 0.7, 0.8, 0.5])

19

20

# set default parameter values and do the fit

21

p0 = np.array([1.0, 1.0])

22

thetaHat, cov = curve_fit(func, x, y, p0, sig, absolute_sigma=True)

23

24

# Retrieve minimized chi-squared, etc.

25

numPoints = len(x)

26

numPar = len(p0)

27

ndof = numPoints - numPar

28

chisq = sum(((y - func(x, *thetaHat))/sig)**2)

29

print ("chisq = ", chisq, ",

ndof = ", ndof)

30

31

# Print fit parameters and covariance matrix

32

print ("\n", "Fitted parameters and standard deviations:")

33

sigThetaHat = np.sqrt(np.diag(cov))

34

for i in range(len(thetaHat)):

35

print ("thetaHat[", i, "] = ", thetaHat[i], "

+-

", sigThetaHat[i])

36

37

print ("\n", "i, j, cov[i,j], rho[i,j]:")

38

for i in range(len(thetaHat)):

39

for j in range(len(thetaHat)):

40

rho = cov[i][j] / (sigThetaHat[i]*sigThetaHat[j])

41

print (i, "

", j, "

", cov[i][j], "

", rho)

42

43

# Set up plot

44

matplotlib.rcParams.update({'font.size':18})

# set all font sizes

45

plt.clf()

46

fig, ax = plt.subplots(1,1)

47

plt.gcf().subplots_adjust(bottom=0.15)

48

plt.gcf().subplots_adjust(left=0.15)

49

plt.errorbar(x, y, yerr=sig, xerr=0, color='black', fmt='o', label='data')

50

plt.xlabel(r'$x$')

51

plt.ylabel(r'$y$', labelpad=10)

52

xMin = 0

53

xMax = 10

54

yMin = 0

55

yMax = 10

56

plt.xlim(xMin, xMax)

57

plt.ylim(yMin, yMax)

58

xPlot = np.linspace(xMin, xMax, 100)

# enough points for a smooth curve

59

fit = func(xPlot, *thetaHat)

60

plt.plot(xPlot, fit, 'red', linewidth=2, label='fit result')

61

62

# Tweak legend

63

handles, labels = ax.get_legend_handles_labels()

64

handles = [handles[1], handles[0]]

65

labels = [labels[1], labels[0]]

66

handles = [handles[0][0], handles[1]]

# turn off error bar for data in legend

22

67

plt.legend(handles, labels, loc='lower right', fontsize=14, frameon=False)

68

69

# Make and store plot

70

plt.show()

71

plt.savefig("simpleFit.pdf", format='pdf')

G. Cowan

RHUL Physics

Version 3.0 / September 2019

23