LEAST- AND CHI-SQUARES FOR THE BUDDING AFICIONADO: ART AND PRACTICE

\begin{matrix}c \\ ◯ \end{matrix}Carl Heiles February 25, 2004

In our never-ending attempt to make your life easier, we present you with this highly instructive, time-saving, and labor-saving informative document! Here we give heuristic derivations, discussions, examples, and the prescription for doing least-squares the easy way using matrix techniques generally, and specifically in IDL. This prescription is given as an example in §5, and the power-user can skip the details and go directly there.

This document is an update, correction, clarification, and elaboration of a previous one made exclusively for the undergraduate lab class. Here we extend the discussion considerably to cover most of what anyone will need in future professional life. This makes the document longer, but the first parts (§1 to 7) are still accessible at the introductory level because they haven’t changed much. We occasionally refer to the books Bevington and Robinson (1992); BR), Cowan (1998), Press et al (2001; Numerical Recipes, NR) and Taylor (1997; T97), and we update the notation to partially conform with NR. We owe sections 11 to 13 to the fantastically excellent website of Stetson, http://nedwww.ipac.caltech.edu/level5/Stetson/Stetson4.html.

We begin with least squares in the classic sense, meaning we minimize the sum of squares instead of minimizing χ^{2}. In astronomy, more often than not you don’t have an independent assessment of the intrinsic uncertainty in the data, which means you cannot evaluate χ^{2}, and the least squares approach is the only option. However, often in astronomy you do want to weight observations differently, e.g. because of integration time, and this requires an approach similar to the χ^{2} one. In later sections we generalize to the χ^{2} and this other weighted-observations case.

Contents

1 LEAST-SQUARE FITTING FOR TWO PARAMETERS, AS WITH A STRAIGHT LINE. 4

1.1 The closed-form expressions for a straight-line fit . . . . . . . . . . . . . . . . . . . . 4

1.2 Better is the following generalized notation. . . . . . . . . . . . . . . . . . . . . . . . 4

2 LEAST-SQUARE FITTING FOR MANY PARAMETERS, AS WITH A CU-BIC 5

3 FAR, FAR BEST AND EASIEST: MATRIX ALGEBRA. 6

4 UNCERTAINTIES IN THE DERIVED COEFFICIENTS. 8 5 A NUMERICAL EXAMPLE AND ITS SOLUTION IN IDL. 10

5.1 Generation of the numerical example . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 5.2 Solution of the Numerical Example in IDL . . . . . . . . . . . . . . . . . . . . . . . . 11 5.3 Discussion of the numerical example . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

6 THE COVARIANCE MATRIX AND ITS NORMALIZED COUNTERPART 13 6.1 Definition of the normalized covariance matrix . . . . . . . . . . . . . . . . . . . . . 13 6.2 The covariance in our numerical example . . . . . . . . . . . . . . . . . . . . . . . . 15

7 REJECTING BAD DATAPOINTS I.: CHAUVENET’S CRITERION 17

8 NONLINEAR LEAST SQUARES 19

9 CHI-SQUARE FITTING, AND WEIGHTED FITTING: DISCUSSION IG-NORING COVARIANCE 22 9.1 The weighted mean: the simplest chi-squared fit . . . . . . . . . . . . . . . . . . . . 23 9.2 Another way to find s_{mean}^{2} . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 9.3 The multivariate Chi-square fit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 9.4 Persnickety Diatribe on Choosing σ_{meas} . . . . . . . . . . . . . . . . . . . . . . . . 27

9.4.1 Choosing and correcting σ_{data,m} . . . . . . . . . . . . . . . . . . . . . . . . . 27 9.4.2 Always use this procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 9.4.3 Think about your results! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 9.4.4 When your measurements are correlated. . . . . . . . . . . . . . . . . . . . . . 28

9.5 The case in which datapoints have different dispersions, or different weights: like a

weighted average . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 9.6 The explicit IDL for the weighted and chi square cases . . . . . . . . . . . . . . . . . 30

10 CHI-SQUARE FITTING, AND WEIGHTED FITTING: DISCUSSION IN-CLUDING COVARIANCE 33 10.1 Phenomenological description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 10.2 Calculating the uncertainties of a single parameter—gedanken experiment . . . . . 35 10.3 Calculating the uncertainties of two parametersr—gedanken experiment . . . . . . . 36 10.4 Calculating the uncertainties of three parametersr—gedanken experiment . . . . . . 36 10.5 Doing these calculations the non-gedanken easy way . . . . . . . . . . . . . . . . . . 36 10.6 Important comments about uncertainties . . . . . . . . . . . . . . . . . . . . . . . . 38

11 REJECTING BAD DATAPOINTS II.: STETSON’S METHOD PLUS CHAU-VENET’S CRITERION 39 11.1 Stetson’s sliding weight . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 11.2 Implementation of the weight in our matrix equations . . . . . . . . . . . . . . . . . 41

12 MEDIAN, INSTEAD OF LEAST SQUARE, FITTING 41 12.1 The median and the double-sided exponential pdf . . . . . . . . . . . . . . . . . . . 42 12.2 Doing a “chi square” fit for the double-sided exponential . . . . . . . . . . . . . . . 43 12.3 IDL’s resources for median fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

13 FITTING WHEN ALL VARIABLES HAVE UNCERTAINTIES 44 13.1 A preliminary: Why the slope is systematically small . . . . . . . . . . . . . . . . . . 45 13.2 Stetson’s method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

13.2.1 Philosophy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 13.2.2 Direct solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

13.3 Commentary on Stetson’s solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 13.3.1 This solution is general and makes sense . . . . . . . . . . . . . . . . . . . . . 51 13.3.2 Yet another point I do not understand . . . . . . . . . . . . . . . . . . . . . 51

13.4 Generalization to Multivariate Case . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

14 USING SINGULAR VALUE DECOMPOSITION (SVD) 54 14.1 Phenomenonological description of SVD . . . . . . . . . . . . . . . . . . . . . . . . . 55 14.2 Important Conclusion for Least Squares!!! . . . . . . . . . . . . . . . . . . . . . . . . 57 14.3 Obtaining [α]^{−1} in IDL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 15 BRUTE FORCE CHI SQUARED AND THE CURVATURE MATRIX 58

15.1 Parameter Uncertainties in Brute Force Chi square fitting . . . . . . . . . . . . . . . 58

16 NOTATION COMPARISON WITH NUMERICAL RECIPES 58

1. LEAST-SQUARE FITTING FOR TWO PARAMETERS, AS WITH A STRAIGHT LINE.

1.1. The closed-form expressions for a straight-line fit

First consider the least squares fit to a straight line. Let y_{m} be the m^{th} measurement of the observed quantity (in this example, y_{m} is zenith distance; t_{m} be the time of the m^{th} measurement; M = the total number of observations, i.e. m runs from 0 to M − 1. Remember that in the leastsquares technique, quantities such as t_{m} are regarded to be known with high accuracy while the quantity y_{m} has uncertainties in its measurement.

We expect the zenith distance y_{m} to change linearly with time as follows:

A + Bt_{m} = y_{m} .  (1)

Given this, one does the maximum likelihood (ML) estimate assuming Gaussian statistics. When all measurements have the same intrinsic uncertainty, this is the same as looking for the solution that minimizes the sum of the squares of the residuals (which we will define later). This leads to the pair of equations (Taylor 8.8, 8.9), called the normal equations

AN + B ∑ t_{m} = ∑ y_{m}  (2a)
A ∑ t_{m} + B ∑ t_{m}^{2} = ∑ t_{m}y_{m}  (2b)

Two equations and two unknowns—easy to solve! The closed-form equations for (A, B) are Taylor’s equations 8.10 to 8.12.

1.2. Better is the following generalized notation.

We want a way to generalize this approach to include any functional dependence on t and even other variables, and to have an arbitrarily large number of unknown coefficients instead of just the two (A, B). This is very easy using matrix math. We will ease into this matrix technique gently, by first carrying through an intermediate stage of notation.

First generalize the straight-line fit slightly by having two functional dependences instead of one. We have something other than the time t_{m}; call it s_{m}. For example, we could have s_{m} = cos(t_{m}) or s_{m} = t_{m}^{2}; or we could have s_{m} = x_{m}, where x_{m} is the position from which the observation was taken. To correspond to equation 1, s_{m} = 1. Then we rewrite equation 1 to include this extra dependence

As_{m} + Bt_{m} = y_{m} .  (3)

There are still only two unknown parameters, so this is an almost trivial generalization; later we’ll generalize to more parameters.

We have M equations like equation 3, one for each measurement. They are known as the equations of condition because they are the equations that specify the theoretical model to which we are fitting the data. There are M equations of condition and only two unknowns (A and B). This is too many equations! We have to end up with a system in which the number of equations is equal to the number of unknowns.

To accomplish this, from equation 3 we form the normal equations. The number of normal equations is equal to the number of unknowns, so in this case we will have two. We could carry through the same ML derivation to derive equations equivalent to equation 2; the result is

A ∑ s_{m}^{2} + B ∑ s_{m}t_{m} = ∑ s_{m}y_{m}  (4a)
A ∑ s_{m}t_{m} + B ∑ t_{m}^{2} = ∑ t_{m}y_{m} .  (4b)

We can rewrite these equations using the notation [st] = ∑ s_{m}t_{m}, etc.:

A[s^{2}] + B[st] = [sy]  (5a)
A[st] + B[t^{2}] = [ty] .  (5b)

This is, of course, precisely analogous to equation 2. And now it’s clear how to generalize to more parameters!

2. LEAST-SQUARE FITTING FOR MANY PARAMETERS, AS WITH A CUBIC

With this notation it’s easy to generalize to more (N) unknowns: the method is obvious because in each equation of condition (like equation 3) we simply add equivalent additional terms such as Cu_{m}, Dv_{m}, etc; and in the normal equations (equation 5) we have more products and also more normal equations.

Let’s take an example with four unknowns (N = 4), which we will denote by A, B, C, D; this would be like fitting a cubic. With N = 4 we need at least five data points (M = 5), so there must be at least five equations of condition. The generalization of equation 4 is the M equations

As_{m} + Bt_{m} + Cu_{m} + Dv_{m} = y_{m} ,  (6)

with m = 0 → (M − 1). Again, the least squares fitting process assumes that the s_{m}, t_{m}, u_{m}, v_{m} are known with zero uncertainty; all of the uncertainties are in the measurements of y_{m}. We then form the four normal equations; the generalization of equation 5 written in matrix format is:

[ \begin{matrix}[ss] & [st] & [su] & [sv] \\ [ts] & [tt] & [tu] & [tv] \\ [us] & [ut] & [uu] & [uv] \\ [vs] & [vt] & [vu] & [vv] \end{matrix} ] [ \begin{matrix}A \\ B \\ C \\ D \end{matrix} ] = [ \begin{matrix}[sy] \\ [ty] \\ [uy] \\ [vy] \end{matrix} ]  (7)

The N × N matrix on the left is symmetric. With N equations and N unknowns, you can actually solve for A, B, C, D!

3. FAR, FAR BEST AND EASIEST: MATRIX ALGEBRA.

The above equations are terribly cumbersome to solve in a computer code because they require lots of loops. However, it becomes trivial if we use matrices. Here we designate a matrix by boldface type.

We illustrate the matrix method by carrying through the above N = 4 example, and we assume that there are 5 independent measurements (M = 5). We first define the matrices

X_{σ} = [ \begin{matrix}s_{0} & t_{0} & u_{0} & v_{0} \\ s_{1} & t_{1} & u_{1} & v_{1} \\ s_{2} & t_{2} & u_{2} & v_{2} \\ s_{3} & t_{3} & u_{3} & v_{3} \\ s_{4} & t_{4} & u_{4} & v_{4} \end{matrix} ]  (8a)

a = [ \begin{matrix}A \\ B \\ C \\ D \end{matrix} ]  (8b)

Y_{σ} = [ \begin{matrix}y_{0} \\ y_{1} \\ y_{2} \\ y_{3} \\ y_{4} \end{matrix} ]  (8c)

so, in matrix form, the equations of condition (equation 6) reduce to the single matrix equation

X_{σ} · a = Y_{σ} .  (9)

The notation for these equations corresponds to NR’s. We write them with subscripts σ to emphasize that they are calculated without dividing by σ_{meas}, i.e. that we are doing least squares instead of chi square fitting. For chi square fitting, see §9 and 10.

Our matrix X_{σ} corresponds to NR’s “design matrix” A of Figure 15.4.1, except that our elements are not divided by σ_{meas,m}, and the matrix equation of condition (equation 9) is identical to the expression inside the square brackets of NR’s equation 15.4.6. The differences arise because here we are discussing least squares fitting instead of chi square fitting, i.e. we have omitted the factors involving σ_{meas,m}, the intrinsic measurement uncertainties (§9).

Again, there are more equations than unknowns so we can’t solve this matrix equation directly. So next we form the normal equations from these matrices. In matrix form, the normal equations (equation 7) reduce to the single equation

[α_{σ}] · a = [β_{σ}] ,  (10)

(NR equation 15.4.10), where

[α_{σ}] = X_{σ}^{T} · X_{σ}  (11a)
[β_{σ}] = X_{σ}^{T} · Y_{σ}  (11b)

The matrix [α_{σ}] is known as the curvature matrix because each element is twice the curvature of σ^{2} (or χ^{2}) plotted against the corresponding product of variables.

The number of equations is equal to the number of unknowns, so the solution of the matrix equation is easy—just rewrite it by multiplying both sides by the inverse of [α_{σ}] (that is, by [α_{σ}]^{−1}), which gives

a = [α_{σ}]^{−1}·[β_{σ}] .  (12)

All of these matrix operations are trivially easy in IDL (§5).

4. UNCERTAINTIES IN THE DERIVED COEFFICIENTS.

How about the uncertainties in the derived quantities contained in the matrix a?

The first thing to do is derive the sample variance s^{2} (the square of standard deviation, or mean error, or dispersion, etc) of the individual data points using the generalization of the usual definition for a straight average of x, s^{2} = [∑_{0}^{M−1}(x_{m} − \overline{x_{m}})^{2}/(M − 1)]. The generalization is, simply, to replace the M − 1 in the denominator by ν = M − N. In the straight-average case, N = 1 so this fits. Here ν is known as the number of degrees of freedom and N, the number of unknown coefficients, is known as the number of constraints. So we have

s^{2} = \frac{1}{M − N} ∑_{m=0}^{M−1}(y_{m} − \overline{y_{m}})^{2} ,  (13)

where \overline{y_{m}} are the values for y_{m} predicted by the derived quantities a. Note the difference: y_{m} are the observed values, while \overline{y_{m}} are the values predicted by the least squares fit. The predicted values are those that are computed from the derived coefficients A, B, C… The M quantities

δy_{m} = y_{m} − \overline{y_{m}}  (14)

are called the residuals or deviations from the fit.

It’s worth reiterating some essentials about s^{2}, and in particular the denominator (M − N). First consider the case of a single-parameter fit, e.g. N = 1. Then we cannot possibly derive a sample variance from only one measurement M = 1; but we can from two M = 2. So the denominator makes sense from that standpoint. The same goes for N > 1.

Next consider the effect of using (M − N) in the denominator: it increases s^{2} by the ratio \frac{M}{M−N} over what you’d get if you just took a straight average and used M. This compensates for the fact that you are subtracting \overline{y_{m}}, which is derived from the data, instead of the truly correct value y^{∗}. (In formal statistical language, y^{∗} is the mean of the parent population from which your set of measurements is drawn). If you used the truly correct value y^{∗}, then the sum would be larger than when using \overline{y_{m}}. The use of M − N in the denominator compensates for this larger value in exactly the right way: the expectation value E_{(s^{2})} for a large number of experiments is precisely equal to the normal variance σ^{2}, which you’d get by using [y^{∗} and M] instead of [\overline{y_{m}} and (M − N)] in equation 14; see Cowan equations 5.9 and 5.10. So s^{2} is, in fact, exactly the number we want: an unbiased estimate of the true variance of our sample. Why not use [y^{∗} and M] in equation 14? The reason is obvious: we don’t know y^{∗}! (If we did, we wouldn’t be doing this analysis!)

It’s easy to calculate the \overline{y_{m}} with matrices. First define the matrix \overline{Y_{σ}} that contains these values:

\overline{Y_{σ}} = [ \overline{y_{0}} ; \overline{y_{1}} ; \overline{y_{2}} ; \overline{y_{3}} ; \overline{y_{4}} ]  (15)

Calculating \overline{Y_{σ}} is simple:

\overline{Y_{σ}} = X_{σ} · a .  (16)

Note that X_{σ} is already defined (equation 8) and a was solved for in equation 12. It’s convenient to define the residual matrix

δY_{σ} = Y_{σ} − \overline{Y_{σ}}  (17)

so we can write

s^{2} = \frac{1}{M − N} δY_{σ}^{T} · δY_{σ} .  (18)

This is the sample variance of the datapoints, not the variances in the derived coefficients. We can obtain these as before, by generalizing the results from the two-parameter case like the straight-line fit discussed in §1. We won’t go through the derivation here; you can find it in Taylor §8.4 and equation 8.16, 8.17. The result is

s_{a}^{2} = s^{2}diag{[α_{σ}]^{−1}}  (19)

Or, to put it simply in words: to get the variance of coefficient n in the matrix a, multiply s^{2} by the n^{th} diagonal element of [α_{σ}]^{−1}.

Although the above equation for s_{a}^{2} is correct, there is more to the story because of covariances, which are the off-diagonal elements. We return to this topic in §6 and §10.

5. A NUMERICAL EXAMPLE AND ITS SOLUTION IN IDL.

If the following sounds like Greek to you, take a look at §3 and 4.

5.1. Generation of the numerical example

Suppose that we make four measurements of the angle y and we want to fit to a parabolic function in time t. In the notation of equation 6, s would be unity, t the time, and u the time squared, so the number of unknowns is three (N = 3). Because there are four independent measurements (M = 4) the subscripts run from m = 0 → 3. Suppose that the four values of time are 5, 7, 9, 11.

First we create the matrix X in IDL

X = fltarr(N, M) = fltarr(3, 4)  (20)

and then we populate it with numbers. In your own work, you would normally do this by reading a data file and transferring the numbers to the matrix using IDL commands; to work through this example, you might manually type them in. After populating the matrix, in direct correspondence with equation 8a we have s_{m} = 1, t_{m} = time_{m}, u_{m} = time_{m}^{2}:

X = [ \begin{matrix}1 & 5 & 25 \\ 1 & 7 & 49 \\ 1 & 9 & 81 \\ 1 & 11 & 121 \end{matrix} ]  (21a)

Suppose that the four measured values of y are (equation 8c)

Y = [ \begin{matrix}142 \\ 168 \\ 211 \\ 251 \end{matrix} ] .  (22a)

Fig. 1.— Our numerical example. Stars are the four datapoints; the solid line is the fit. We perform two fits: one uses the original definition of time; the other uses (time − 8), in effect moving the y-axis to the dashed line. The two fits give the same line but the coefficients and their errors differ greatly.

Figure 1 shows the datapoints, together with the fitted curve.

One word of caution here: in IDL, to get these into a column matrix, which is how we’ve treated Y above, you have to define Y as a two-dimensional array because the second dimension represents the column. When working in IDL it’s more convenient to define a row vector, which has only one dimension; in IDL you do this by defining Y = [142, 168, 211, 251]; you can make it into the necessary column vector by taking its transpose, i.e. Y = transpose(Y).

5.2. Solution of the Numerical Example in IDL

In IDL we calculate the normal equation matrices and denote the [α] in equation 11a by XX:

XX = transpose(X)##X ,  (23a)

and we denote the [β] in equation 11b by XY:

XY = transpose(X)##Y .  (23b)

In IDL we take the inverse of [α] (same as XX) by

XXI = invert(XX) .  (24)

The least-square fitted quantities are in the matrix a (equation 12), which we obtain in IDL with

a = XXI ## XY .  (25)

In IDL we denote the matrix of predicted values \overline{y_{m}} by YBAR, which is

YBAR = X ## a .  (26)

and we can also define the residuals in Y as

DELY = Y − YBAR  (27)

In IDL we denote s^{2} in equations 13 and 18 by s sq and write

s sq = transpose(DELY)##DELY/(M − N) .  (28a)

or

s sq = total(DELY ∧ 2)/(M − N) .  (28b)

It is always a good idea to plot all three quantities (the measured values Y, the fitted values YBAR, and the residuals DELY) to make sure your fit looks reasonable and to check for bad datapoints.

To get the error in the derived coefficients we need a way to select the diagonal elements of a matrix. Obviously, any N × N matrix has N diagonal elements; a convenient way to get them is

diag elements of XXI = XXI[(N + 1) ∗ indgen(N)]  (29)

In IDL, we define the variances of the N derived coefficients by vardc (think of “variances of derived coefficients”). You can get this as in equation 19 from^{1}

vardc = s sq ∗ XXI[(N + 1) ∗ indgen(N)] .  (30)

5.3. Discussion of the numerical example

For this numerical example, the solution (the array of derived coefficients) is

a = [ \begin{matrix}96.6250 \\ 4.5000 \\ 0.8750 \end{matrix} ]  (31a)

and the errors in the derived coefficients [the square root of the σ^{2}’s of the derived coefficients, i.e. [σ_{n}^{2}]^{1/2} or, in IDL, sqrt(vardc) in equations 30] are:

σ_{A} = [ \begin{matrix}34.012 \\ 9.000 \\ 0.5590 \end{matrix} ] .  (31b)

These results look horrible: the uncertainties are large fractions of the derived coefficients,

The reason: we have specifically chosen an example with terrible covariance. And the great thing is this can be fixed easily (at least in this case—certainly not always), without taking more data!

6. THE COVARIANCE MATRIX AND ITS NORMALIZED COUNTERPART

First we provide a general discussion, then we apply it to the above numerical example.

6.1. Definition of the normalized covariance matrix

The variances in the derived coefficients are obtained from the diagonal elements of XXI. The off-diagonal elements represent the covariances between the derived coefficients. Covariance means, specifically, the degree to which the uncertainty in one derived coefficient affects the uncertainty in another derived coefficient.

1If you used equation 28a instead of 28b, then IDL considers s sq an array and you need to use a # instead of a ∗ in this equation.

Because the covariance matrix elements relate pairwise the various coefficients, the units of the matrix elements are all different. This makes it convenient to reduce all the matrix elements to a standard set of units—namely, no units at all. So before discussing the covariance matrix per se, we first discuss its normalized counterpart.

The normalized covariance matrix^{2} ncov is derived from XXI by dividing each element by the square root of the product of the corresponding diagonal elements. Let ncov be the normalized covariance matrix; then

ncov_{ik} = \frac{XXI_{ik}}{\sqrt{XXI_{ii} XXI_{kk}}}  (32)

This is the same normalization that one does with the Pearson linear correlation coefficient of two variables. In fact, the elements of the normalized covariance matrix are the correlation coefficients. In IDL, you do the following:

dc = XXI[(N + 1) ∗ indgen(N)]  (33a)
ncov = XXI/sqrt(dc##dc)  (33b)

In the above, dc##dc is an N × N matrix consisting of products of the diagonals of XXI, so dividing XXI by sqrt(dc##dc) generates the normalized version.

Because ncov is a normalized covariance matrix, you might think that it’s non-normalized parent is XXI—and you’d be almost right. For the least squares case we are discussing, the true covariance matrix C is^{3}

C = s^{2} XXI .  (34)

In ncov, the diagonal elements are all unity and the off-diagonal elements reflect the interdependence of the derived coefficients on each other. The off-diagonal elements can range from −1 → 1. Each matrix element is the correlation coefficient between the uncertainties of its two parameters. In particular, suppose that the data happens to produce a coefficient that differs from its true value by some positive number. If the normalized matrix element is negative, then the other coefficient will tend to differ from its true value by a negative number.

2It is a pleasure to thank Doug Finkbiener for introducing me to this concept.

3For chi square, you use σ_{meas}^{2} instead of s^{2}; see §9.

Here’s a more detailed discussion of what the covariance means. Suppose you are least-squares fitting for two derived coefficients (A_{0} and A_{1}). When you do a least-squares fit to a set of data, you are fitting one set of data out of a possible infinity of possible sets that you’d get by repeating the experiment, and your particular set of data happens to produce specific values of \overline{A_{0}} and \overline{A_{1}}, which differ from the true values (A_{0}^{∗}, A_{1}^{∗}) by amounts δA_{0} and δA_{1}. If their covariance is zero, then in the infinity of data sets you’d find that δA_{0} is uncorrelated with δA_{1}. But if it is nonzero, these two quantities would be correlated.

A high covariance is bad because the derived variables depend on each other. For one, this means that with noisy data power can be shared or passed from one parameter to/from its covariant counterpart(s). As we shall see §10, it also significantly influences the uncertainties in derived coefficients. Often a high covariance results from a poor choice of functions that you are fitting or even a bad choice of the zero point of the independent variable—as in our numerical example (see the next subsection). And, as in that example, you can sometimes eliminate the bad covariance by reformulating the problem—you don’t even need to take more data! The best reformulation involves using a set of orthonormal functions. However, sometimes your interest is a specific set of functions that are not orthogonal, and in such cases it makes no sense to convert to orthogonal functions—because you just have to convert back again and do the error propagation after-the-fact instead of letting the least squares process do it for you.

6.2. The covariance in our numerical example

Apply equation 33 to obtain the covariance matrix for our numerical example:

ncov = [ \begin{matrix}1 & −.989848 & .969717 \\ −.989848 & 1 & −.993808 \\ .969717 & −.993808 & 1 \end{matrix} ] .  (35)

The off-diagonal elements are huge. This is the reason why our solution is so bad.

In this seemingly innocuous example we have an excellent case of a poor choice of zero point for the independent variable (the time). The reason is clear upon a bit of reflection. We are fitting for y = A_{0} + A_{1}t + A_{2}t^{2}. The coefficient A_{0} is the y-intercept and A_{1} is the slope. Inspection of Figure 1 makes it very clear that an error in the slope has a big effect on the y-intercept.

Now we retry the example, making the zero point of the time equal to the mean of all the times, that is we set (time_{m} = time_{m} − 8). We get the same fitted line, but the derived coefficients are completely different—and amazingly better! We get

A = [ \begin{matrix}188.625 \\ 18.500 \\ 0.87500 \end{matrix} ]  (36a)

σ_{A} = [ \begin{matrix}3.58 \\ 1.00 \\ 0.559 \end{matrix} ] .  (36b)

In redefining the origin of the independent variable, we have made the zero intercept completely independent of the slope: changing the slope has no affect at all on the intercept. You can see this from the normalized covariance matrix, which has become

ncov = [ \begin{matrix}1 & 0 & −0.78086881 \\ 0 & 1 & 0 \\ −0.78086881 & 0 & 1 \end{matrix} ] .  (37)

which is nice, but not perfect: Our step is partial because the second-order coefficient A_{2} affects A_{0} because, over the range of [(time − 8) = −3 → +3], the quantity [A_{2} Σ(time_{m} − 8)^{2}] is always positive and is thereby correlated with A_{0}.

We could complete the process of orthogonalization by following the prescription in BR chapter 7.3, which discusses the general technique of orthogonalizing the functions in least square fitting. The general case is a royal pain, analytically speaking, so much so that we won’t even carry it through for our example. But for numerical work you accomplish the orthogonalization using Singular Value Decomposition (SVD), which is of course trivial in IDL (§14).

For some particular cases, standard pre-defined functions are orthogonal. For example, if t_{m} is a set of uniformly spaced points between (−1 → 1) and you are fitting a polynomial, then the appropriate orthogonal set is Legendre polynomials. This is good if your only goal is to represent a bunch of points by a polynomial function, because the coefficients of low-order polynomials are independent of the higher ones. However, it’s more work and, moreover, often you are interested in the coefficients for specific functions that don’t happen to be orthogonal; in such cases, you should just forge ahead.

But always look at the normalized covariance matrix. Suppose one pair of off-diagonal elements departs significantly from zero. Then their corresponding functions are far from being orthogonal and the variances of the derived coefficients will suffer as a result. You might be able to eliminate one of the parameters to make the fit more robust. For example, suppose one function is t cos(t) and the other is sin(t) cos(t). If the range of t is small, these two functions are indistinguishable and have a large covariance; you should eliminate one from the fit. If the range of t is large, there is no problem.

For further discussion of covariance, see §10. Also, you might also want to try out another example in Taylor’s §8.5.

7. REJECTING BAD DATAPOINTS I.: CHAUVENET’S CRITERION

Least squares fitting is derived from the maximum likelihood argument assuming the datapoint residuals δy_{m} have a Gaussian pdf. This means that the errors are distributed as

p_{(δy;σ)} = \frac{1}{\sqrt{2π}σ} e^{−( \frac{δy^{2}}{2σ^{2}} )}  (38)

where σ^{2} is the true variance of the datapoints, i.e. s^{2} in equation 13 (to be precise, s^{2} needs to be averaged over many experiments).

More importantly, the probability of finding datapoints inside the limits ±Δy is

P_{(|δy|<Δy)} = ∫_{−Δy}^{+Δy} p_{(δy;σ)}d(δy) = erf (\frac{Δy}{\sqrt{2}σ} )  (39)

where we use the commonly-defined error function erf_{(X)} = \frac{1}{\sqrt{π}} ∫_{−X}^{+X} e^{−x^{2}}dx. A particularly important value is for Δy = σ, for which

P_{(|δy|<σ)} = 0.683  (40)

If we have an experiment with M datapoints, then the number of datapoints we expect to lie outside the interval ±Δy is

M_{(outside Δy)} = M [1 − erf (\frac{Δy}{\sqrt{2}σ} )]  (41)

Chauvenet’s criterion simply says:

1. Find Δy such that M_{(outside Δy)} = 0.5. This is given by

\frac{Δy}{σ} = erf^{−1} (1 − \frac{1}{2M} )  (42)

2. Discard all datapoints outside this range.

This criterion makes very good sense. It leads to the following rejection criteria

Chauvenet’s criterion versus M

M  \frac{Δy}{σ}
100  2.81
1000  3.48
10^{4}  4.06
10^{5}  4.56

which is a moderately interesting set of numbers.

We offer the following important Comments:

• This assumes data are Gaussian distributed. In real life this doesn’t often happen because of “glitches”. Examples of glitches can be interference in radio astronomy, meteors in optical astronomy, and cosmic rays on CCD chips. These glitches produce bad points that depart from Gaussian statistics. They are often called outliers. It is very important to get rid of the outliers because the least squares process minimizes the squares of the residuals. Outliers, being the points with the largest residuals, have a disproportionately evil effect on the result. On the other hand, if your data don’t follow Gaussian statistics as their intrinsic pdf, then you should think twice before using least squares! (Like, maybe you should try the median fitting discussed in §12).

• You may wish to relax Chauvenet’s criterion by increasing the Δx beyond which you discard points. This is being conservative and, in the presence of some nonGaussian statistics, not a bad idea. But think about why you are doing this before you do it. Maybe the intrinsic statistics aren’t Gaussian? You should never make Chauvenet’s criterion more stringent by decreasing the Δx beyond which you discard points. This rule hardly needs elaboration: it means you are discarding datapoints that follow the assumed pdf!

• Most statistics books (e.g. Taylor, BR) harp on the purity aspect. One extreme: don’t throw out any datum without examining it from all aspects to see if discarding it is justified. The other extreme: apply Chauvenet’s criterion, but do it only once and certainly not repeatedly. Being a real-life astronomer, our approach is different. There do exist outliers. They increase the calculated value of σ. When you discard them, you are left with a more nearly perfect approximation to Gaussian statistics and the new σ calculated therefrom will be smaller than when including the outliers. Because the original σ was too large, there may be points that should have been discarded that weren’t. So our approach is: repeatedly apply Chauvenet’s criterion until it converges. If it doesn’t converge, or if it discards an inordinately large number of datapoints, you’ve got real problems and need to look at the situation from a global perspective.

• Many observers use the 3σ criterion: discard any points with residuals exceeding 3σ. This is definitely not a good idea: the limit 3σ is Chauvenet’s criterion for M = 185 datapoints. Very often M exceeds this, often by a lot.

• To apply Chauvenet’s criterion it’s most convenient to calculate the inverse error function. For this, you have two choices. One (for sissies like myself), you can use inverf.pro from my area ˜ heiles/idl/gen . But the real he-man will want to learn about using a root-finding algorithm such Newton’s method (NR §9.4 and 9.6) together with the error function; both procedures exist in IDL as newton and errorf. You at least ought to skim lightly some of NR’s chapter 9 about root finding, because some day you’ll need it.

8. NONLINEAR LEAST SQUARES

The least squares formulation requires that the data values depend linearly on the unknown coefficients. For example, in equation 1, the unknown coefficients A and B enter linearly.

Suppose you have a nonlinear dependence, such as for example wanting to solve for A and B with equations of condition that look like

sin(At_{m}) + Bt_{m} = y_{m} .  (43)

What do you do here? You linearize the process, using the following procedure.

First, assume a trial values for A and B; call these A_{0} and B_{0}. You should pick values that are close to the correct ones. In our example you wouldn’t need to do this for B, but it’s easier to treat all coefficients identically. These trial values produce predicted values y_{0,m}:

sin(A_{0}t_{m}) + B_{0}t_{m} = y_{0,m} .  (44)

Subtract equation 44 from 43, and express the differences as derivatives. Letting δA = A − A_{0} and δB = B − B_{0}, this gives

δA[t_{m} cos(A_{0}t_{m})] + δBt_{m} = y_{m} − y_{0,m} .  (45)

This is linear in (δA, δB) so you can solve for them using standard least squares. Increment the original guessed values to calculate A_{0,new} = A_{0} + δA and B_{0,new} = B_{0} + δB, These won’t be exact because higher derivatives (including cross derivatives) come into play, so you need to use these new values to repeat the process. This is an iterative procedure and you keep going until the changes become “small”. The generalization to an arbitrarily large number of unknown coefficients is obvious.

We now offer some cautionary and practical remarks.

(0) In linear least squares, the curvature and covariance matrices are set by the values of the independent variable, which here is denoted by t, and are independent of the datapoint values. Here, the matrix elements change from one iteration to the next because they depend on the guessed parameters, and sometimes they even depend on the datapoint values.

(1): Multiple minima: Nonlinear problems often have multiple minima in σ^{2}. A classical case is fitting multiple Gaussians to a spectral line profile. Gaussians are most definitely not orthogonal functions and in some cases several solutions may give almost comparably good values of σ^{2}, each one being a local minimum. For example, for the case of two blended Gaussians, one can often fit two narrow Gaussians or the combination of a wide and narrow Gaussian, the two fits giving almost equal σ^{2}. The lower of these is the real minimum but, given the existence of systematic errors and such, not necessarily the best solution. The best solution is often determined by physical considerations; in this case, for example, you might have physical reasons to fit a broad plus narrow Gaussian, so you’d choose this one even if it’s σ^{2} weren’t the true minimum.

(2): The Initial Guess: When there are multiple minima, the one to which the solution converges is influenced by your initial guess. To fully understand the range of possible solutions, you should try different initial guesses and see what happens. If the solutions always converge to the same answer, then you can have some confidence (but not full confidence) that the solution is unique.

(3): Iterative stability: If your initial guess is too far from the true solution, then the existence of higher derivatives means that the computed corrections can be too large and drive the iterative solution into instability. It is often a good idea to multiply the derived correction factors (δA and δB above) by a factor F less than unity, for example F = 0.5 or 0.75. This increases the number of iterations required for convergence but often allows convergence instead of producing instability.

(4): Convergence criteria: How do you know when the solution has converged? One way: for each iteration, calculate the uncertainties in the derived coefficients. If the uncertainty exceeds the correction, then you are getting close. An alternate way, which I usually use: if the fractional correction (e.g. \frac{δA}{A_{0}}) decreases below some threshold, say 1%, you’re close (some parameters, such as angles, need a threshold that is absolute instead of fractional). At this point, if you are using F ≠ 1, set F = 1, do a few more iterations, and you’re done.

(5): Numerical derivatives: Sometimes the equations of condition are so complicated that taking the derivatives, as in obtaining equation 45, is a huge job and subject to mistakes. So you can take numerical derivatives instead of analytic ones. Be careful, though; it’s safest to use double precision and think a bit about numerical accuracy; take a look at NR’s section 5.7 on evaluating numerical derivatives.

(6): Canned nonlinear least squares: Packages like IDL offer canned nonlinear least squares routines. They are designed to work well for a wide range of different problems. However, often you can do better by tailoring things (such as the factor F and convergence criteria above) for the specific problem at hand. A good example is Gaussian fitting: IDL’s fitting program doesn’t converge for some input data, while for many of these cases the program that I wrote myself works fine.

When convergence is slow or doesn’t occur because your functions are complicated, you might wish to try the Levenberg-Marquardt method (NR §15.5); IDL function LMFIT. This involves increasing the diagonal elements of the curvature matrix by a set of suitably chosen factors; when you get close to the minimum, you reset these factors to unity. I have never had to resort to such tactics. In my opinion, if you are in such difficulty that you need such a technique you should look at the problem with an eye to gaining an appreciation of the approximate locations of the true minima and beginning with your initial guess reasonably close to it. Of course, this might be difficult with a large number of parameters.

(7): Be careful and LOOK at the solution before accepting it! These nonlinear problems can produce surprising results, sometimes completely meaningless results. Don’t rely on them to be automatic or foolproof!

(8): Reformulate! Sometimes you can avoid all this by reformulating the problem. There are two cases: the harmless case and the not-so-harmless case.

An example of the harmless case is fitting for the phase ϕ in the function y = cos(θ + ϕ). This is definitely a nonlinear fit! But its easy to reformulate it in a linear fit using the usual trig identities to write y = A cos θ−B sin θ, where \frac{B}{A} = tan ϕ. Solve for (A, B) using linear least squares, calculate ϕ, and propagate the uncertainties.

An example of the not-so-harmless case is in NR’s §15.4 example: fit for (A, B) with equations of condition y_{m} = Ae^{−Bx_{m}}. They suggest linearizing by rewriting as log(y_{m}) = C − Bx_{m}, solving for (B, C), and deriving A after-the-fact. This is not-so-harmless because you are applying a nonlinear function to the observed values y_{m}; thus the associated errors σ_{meas,m} are also affected. This means you have to do weighted fitting, which is discussed in §9 below. Suppose that A = 1, your datapoints all have σ_{meas,m} = 0.05, and the observed y_{m} ranges from 0.05 to 1. The datapoint with y_{m} = 0.05 has a manageable σ_{meas_{m}}, but what is the corresponding value of σ_{meas,m} for log y_{m} = log 0.05? It’s ill-defined and asymmetric about the central value. Or even, God forbid, you have an observed y_{m} that’s negative??? Even for y_{m} not near zero, you need to calculate new σ_{meas,m} by error propagation; in this case, you need to reassign σ(log y) = \frac{d log y}{dy} σ(y) = \frac{σ(y)}{y}. This is OK when y_{m} is large enough so that the linear approximation is accurate, but if not the converted noise becomes nonGaussian.

You should regard your datapoints as sacrosanct and never apply any nonlinear function to them.

9. CHI-SQUARE FITTING, AND WEIGHTED FITTING: DISCUSSION IGNORING COVARIANCE

In least squares fitting the derived parameters minimize the sum of squares of residuals as in equation 13, which we repeat here

s^{2} = \frac{1}{M − N} ∑_{m=0}^{M−1}(y_{m} − \overline{y_{m}})^{2} .

Chi square fitting is similar except for two differences. One, we divide each residual by its intrinsic measurement error σ_{meas,m}; and two, we define χ^{2} as the sum

χ^{2} = ∑_{m=0}^{M−1} \frac{(y_{m} − \overline{y_{m}})^{2}}{σ_{meas,m}^{2}} .  (46a)

Along with χ^{2} goes the reduced chi squared \overline{χ}^{2} = \frac{χ^{2}}{M−N}

\overline{χ}^{2} = \frac{1}{M − N} ∑_{m=0}^{M−1} \frac{(y_{m} − \overline{y_{m}})^{2}}{σ_{meas,m}^{2}} .  (46b)

which is more directly analogous to the definition of s^{2}.

Chi-square fitting is very much like our least-squares fitting except that we divide each datapoint by its intrinsic measurement uncertainty σ_{meas,m}. Thus, the reduced Chi-square (χ^{2}) is equal to the ratio of the variance of the datapoint residuals (σ^{2}) to the adopted intrinsic measurement variances (σ_{meas}^{2}). So it should be obvious that in Chi-square fitting, you must know the measurement uncertainties σ_{meas} of the individual data points beforehand. If you want to give the various datapoints weights based on something other than σ_{meas}, then that is just like Chi-square fitting except that you can adopt an arbitrary scale factor for the uncertainties (section 9.4).

Chi-squared fitting treats uncertainties of the derived parameters in a surprising way. Getting the coefficient uncertainties with chi-squared fitting is a tricky business because 1. With the standard treatments, the errors in the derived parameters don’t depend on the

residuals of the datapoints from the fit (!).

2. The errors in the derived parameters can depend on their mutual covariances. This discussion requires a separate section, which we provide below in §10.

In this section we treat chi square fitting ignoring covariance. We begin by illustrating the difference between least squares and chi square fitting by discussing the simplest chi-square fitting case of a weighted mean; then we generalize to the multivariate chi square fitting case.

9.1. The weighted mean: the simplest chi-squared fit

First, recall the formulas for an ordinary unweighted average in which the value of each point is y_{m} and the residual of each point from the weighted mean is δy_{m}:

mean = \frac{∑ y_{m}}{M}  (47a)
s^{2} = \frac{∑ δy_{m}^{2}}{M − 1}  (47b)
s_{mean}^{2} = \frac{s^{2}}{M} = \frac{∑ δy_{m}^{2}}{M(M − 1)} ,  (47c)

s_{mean}^{2} is the variance of the mean and s^{2} is the variance of the datapoints around the mean. Recall that in this case the mean is the least squares fit to the data, so to use least squares jargon we can also describe s_{mean} as the error in the derived coefficient for this single-parameter least squares fit.

Now for a weighted average in which the weight of each point is w_{meas,m} = \frac{1}{σ_{meas,m}^{2}} = \frac{1}{σ_{meas}^{2}} . Applying maximum likelihood, in a nonweighted average the quantity that is minimized is ∑ δy_{m}^{2}; in a weighted average the quantity minimized is χ^{2} = ∑ \frac{δy_{m}^{2}}{σ_{meas,m}^{2}} = ∑ w_{meas,m}δy_{m}^{2}. So you’d think that the three equations corresponding to the above would become

mean_{w} = \frac{∑ w_{meas,m}y_{m}}{∑ w_{meas,m}} → \frac{∑ y_{m}}{M}  (48a)
s_{w}^{2} = \frac{M}{M − 1} \frac{∑ w_{meas,m}δy_{m}^{2}}{∑ w_{meas,m}} → \frac{∑ δy_{m}^{2}}{M − 1}  (48b)

s_{mean}^{2} = \frac{s_{w}^{2}}{M} = \frac{∑ w_{meas,m}δy_{m}^{2}}{(M − 1) ∑ w_{meas,m}} → \frac{∑ δy_{m}^{2}}{M(M − 1)} ,  (48c)

where the arrow to the rightmost expressions use our assumption that σ_{meas,m} are all identical. In fact, after a formal derivation, the first two equations are correct. However, after a formal derivation, the last of these equations becomes, and is always written (e.g. BR equation 4.19; Taylor equation 7.12)

s_{mean}^{2} = \frac{1}{∑ w_{meas,m}} → \frac{σ_{meas}^{2}}{M} ,  (49)

This is a problem, for the following reason.

Note the excruciatingly painful difference between equation 48c and equation 49: the former depends on the variance of the datapoint residuals s_{w}^{2}, as you’d think it should, while the latter depends on only the adopted intrinsic measurement variance of the data σ_{meas}^{2}, which is chosen by the guy doing the fit. If you do an unweighted average, and derive a certain variance, and next do a weighted average in which you choose some value for σ_{meas} that happens to be wrong, the two fits give different results for s_{mean}^{2}. This is crazy.

To get around this difficulty, we follow the procedure in BR equations 4.20 to 4.26. This introduces an arbitrary multiplicative factor for the weights and goes through the ML calculation to derive, instead of equation 49, the far superior

s_{mean}^{2} = \frac{χ^{2}}{∑ w_{meas,m}} = \frac{χ^{2}}{(M − 1) ∑ w_{meas,m}} → \frac{s_{w}^{2}}{M} ,  (50)

which is precisely the same as our intuitive guess, equation 48c. The difference between equations 50 and 49 is the numerator, which contains the reduced chi-squared χ^{2}; for this case where all σ_{meas,m} are identical, χ^{2} = \frac{s_{w}^{2}}{σ_{meas}^{2}} . Here, following equation 46, we write

χ^{2} = ∑ w_{meas,m}δy_{m}^{2} = ∑ \frac{δy_{m}^{2}}{σ_{meas,m}^{2}}  (51)

9.2. Another way to find s_{mean}^{2}

There is another way to find the variance in the derived mean. This works only if your specifications of σ_{meas,m} are correct.

1. Calculate χ^{2} from equation 46. In particular, you use the predicted values for \overline{y_{m}}, which are calculated from the derived parameter values (in this case, it’s only the derived mean, which we call a). Denote this value of χ^{2} by the symbol χ_{a}^{2}.

2. Perturb the derived mean by varying amounts δa and recalculate χ^{2}; denote this by the symbol χ_{δa}^{2}.

3. Form the difference Δχ^{2} = χ_{δa}^{2} − χ_{a}^{2}. The value of δa for which Δχ^{2} = 1 is the standard deviation of a (i.e., that value squared is the variance of a).

4. If you are doing a multivariate chi-square fit, then the Δχ^{2} isn’t unity, but larger. For example, when fitting two variables, Δχ^{2} must equal 2.3. We will discuss this situation below. The above procedure is invaluable when specifying confidence ellipses (or multidimensional ellipsoids) in multiparameter fits, as we shall discuss ad nauseum in section 10.

9.3. The multivariate Chi-square fit

In this case, Chi-square fitting is just like least-squares fitting except for the following:

1. In the matrix X_{σ} of equation 8a, each row m is a different measurement with a different intrinsic variance σ_{meas,m}. You are generating a new matrix X, which is identical to X_{σ} except that each row m (which contains a particular equation of condition) is divided by σ_{meas,m}. This new matrix is the same as NR’s design matrix (Figure 15.4.1), which they denote by A.

2. Divide each datapoint y_{m} in equation 8b by σ_{meas,m}. You are generating a new matrix Y, which is identical to Y_{σ} except that each row is divided by σ_{meas}. This new matrix is the same as NR’s vector b.

3. Carry through the matrix calculations in equations 52 below (using the matrices not subscripted with σ).

You’ve divided each row, i.e. the equation of condition for each row m, by a common factor, so the solution of that particular equation of condition is unchanged.

In equations 52 just below, we include also the expressions written in terms of the σ-subscripted least square matrix versions, under the supposition that all datapoints have the same intrinsic measurement uncertainty σ_{meas}. In the full matrix equation 9, which now reads X · a = Y, you’ve divided both sides by the same quantity, σ_{meas}. Obviously, this doesn’t change the derived results for a. (If the σ_{meas,m} differ from each other, then a is affected, of course.)

Let’s summarize the relationships between the least-squares (σ-subscripted) and χ^{2} matrices for this specific case where all σ_{meas,m} are identical:

X = \frac{X_{σ}}{σ_{meas}}  (52a)

Y = \frac{Y_{σ}}{σ_{meas}}  (52b)

X · a = Y  (52c)

[α] = X^{T} · X = \frac{[α_{σ}]}{σ_{meas}^{2}}  (52d)

[β] = X^{T} · Y = \frac{[β_{σ}]}{σ_{meas}}  (52e)

[α]^{−1} = σ_{meas}^{2}[α_{σ}]^{−1}  (52f)

a = [α]^{−1} · [β]  (52g)

χ^{2} = δY^{T} · δY  (52h)

χ^{2} = \frac{δY^{T} · δY}{M − N}  (52i)

For the σ-subscripted matrices, when you calculate the predicted values \overline{Y} in equation 16, using X instead of X_{σ}, the values are divided by σ_{meas}. Because Y is also divided by σ_{meas}, the residuals δY are as well. Thus when you calculate the square of the sum of the residuals as in equation 18 you get \frac{s^{2}}{σ_{meas}^{2}} , which is the reduced chi-square, denoted as χ^{2} :

χ^{2} = \frac{δY_{σ}^{T} · δY_{σ}}{σ_{meas}^{2}(M − N)}  (53a)

and, obviously, you could calculate the non-reduced χ^{2}

χ^{2} = \frac{δY_{σ}^{T} · δY_{σ}}{σ_{meas}^{2}}  (53b)

This is precisely analogous to equation 51 for the weighted mean example above. The reduced chi-square χ^{2} is equal to the ordinary chi-square χ^{2} except that it is divided by the number of degrees of freedom, ordinarily denoted by ν, which is equal to (M − N).

Finally, we have the analogy of equation 50 expressed in matrix form as in equation 19:

s_{a}^{2} = χ^{2}diag{[α]^{−1}}  (54)

This intuitively-derived result is in contrast to the result derived from a formal derivation, which is the analogy to equation 49; again, it omits the χ^{2} factor:

s_{a}^{2} = diag{[α]^{−1}}  (55)

This formally-deriuved result is what’s quoted in textbooks (e.g. NR equation 15.4.15, BR equation 7.25). It provides parameter errors that are independent of the datapoint residuals, and leads to the same difficulties discussed above for the weighted mean case..

We recommend that you use equation 54. Our result, equation 54, is very reasonable. Suppose, for example, that the least-squares fit model is perfect and the only deviations from the fitted curve result from measurement error. Then by necessity we have s^{2} ≈ σ_{meas}^{2} and χ^{2} ≈ 1. (We write “≈” instead of “=” because different experiments produce somewhat different values of s^{2} because of statistical fluctuations; an average gives σ^{2} = ⟨s^{2}⟩). In this situation, though, equations 54 and 55 are identical. However, if the least-squares fit model is not correct, meaning that it doesn’t apply to the data, then the residuals will be larger than the intrinsic measurement errors, which will lead to larger values of χ and χ^{2}.

9.4. Persnickety Diatribe on Choosing σ_{meas}

9.4.1. Choosing and correcting σ_{data,m}

In the previous section, equation 55 taught us that—formally, at least—the variances in the derived fit parameters (or their uncertainties, which are the square roots) depend only on the adopted uncertainties σ_{meas,m} and not on the actual variance of the datapoints.

Are you bothered by the fact that the variances of the derived parameters s_{a} are independent of the data residuals? You should be: it is obvious that the residuals should affect s_{a}.

Formally, s_{a} depends only on the adopted uncertainties σ_{meas,m}, which are chosen beforehand by you—you’re supposed be such a good experimentalist that you really do know the intrinsic uncertainty in your measured values. Moreover, you are assuming that there are no other sources of uncertainty—such as “cosmic scatter” or an inappropriate model to which you are fitting the data. Suppose your adopted values of σ_{meas,m} are off by a common scale factor, i.e. if σ_{meas,adopted} = fσ_{meas,true}. Then χ^{2} ≈ f^{−2} instead of χ^{2} ≈ 1. And to obtain the parameter errors from δχ^{2}, you must find the offset δx such that Δχ^{2} = f^{−2} ≈ χ^{2}.

You can correct for this erroneous common factor f by dividing your adopted values of σ_{meas,m} by f. Of course, you don’t know what this factor f is until you do the chi squared fit. Dividing them by f is equivalent to multiplying them by χ̂. And, of course, the same as multiplying σ_{meas,m}^{2} by χ^{2}.

9.4.2. Always use this procedure

To be kosher, after having run through the problem once with the adopted σ_{meas,m}, calculate the χ^{2}; multiply all σ_{meas,m} by χ̂; and redo the problem so that the new χ^{2} = 1. Then the derived variance s_{a} is also correct. You can obtain it either as the corresponding diagonal to the covariance matrix (equations 54 and 55, which are identical in this case) or by finding what departure from x_{0} is necessary to make Δχ^{2} = 1. This redoing the fit may seem like unnecessary work, but when we deal with multiparameter error estimation it’s the best way to go to keep yourself from getting confused.^{4}

9.4.3. Think about your results!

In the case Δχ ≈ 1 (and χ^{2} ≈ 1) the dispersion of the observed points σ_{m} is equal to the intrinsic dispersion of the datapoints σ_{meas} and the mathematical model embodied in the leastsquares fit is perfect. That, at least, is the theoretical conclusion. In practice, however, your obtaining such a low, good value for χ^{2} might mean instead that you are using too large a value for σ_{meas}: you are ascribing more error to your datapoints than they really have, perhaps by not putting enough faith in your instrument.

But there is another way you can get artificially small values for χ^{2}. This will occur if your measurements are correlated. Suppose, for example, that by mistake you include the same measurements several times in your fit. Then your measurements are no longer independent. Cowan discusses this possibility his §7.6.

High values of χ^{2} indicate that the model is not perfect and could be improved by the use of a different model, such as the addition of more parameters—or, alternatively, that you think your equipment works better than it really does and you are ascribing less error to your datapoints than they really have. And in this case, using equation 55 instead of 54 is disastrous.

Think about your results.

9.4.4. When your measurements are correlated…

One more point, a rather subtle one. There are circumstances in which your datapoints are not independent. Then the formulation of chi-square fitting (and least-squares fitting, for that matter) is more complicated. You need to calculate the covariance matrix for the measured values y_{m}; call this covariance matrix V. If this matrix is not unitary, then χ^{2} is no longer given by equation 52h.

4The Δχ^{2} = 1 condition only applies for the single-variable case. For more variables, the variances need to make Δχ^{2} > 1; for example, if you are fitting for two parameters, then the variances make Δχ^{2} = 2.3, as we discuss in §10.

Rather, it is given by

χ^{2} = δY^{T} · V^{−1} · δY  (56a)

Of course, this leads to a different expression for a, which replaces equation 52g,

a = (X^{T} · V^{−1} · X)^{−1}·X^{T} · V^{−1} · Y ,  (56b)

and also to a different equation for the covariance matrix, which replaces equation 52f,

[α]^{−1} = (X^{T} · V^{−1} · X)^{−1} .  (56c)

Correlated datapoints can occur when the measured y_{m} are affected by systematic errors or instrumental effects. Cowan §7.6 discusses this case. For example, suppose you take an image with a known point-spread function (psf) and want to fit an analytic function to this image. Example: a background intensity that changes linearly across the field plus a star. Here the independent variables in the function are the (x, y) pixel positions and the data are the intensities in each pixel. You’d take the intensity in each individual pixel and fit the assumed model. But here your data values are correlated because of the psf. Because you know the psf, you know the correlation between the various pixels. Such a formulation is required for CBR measurements because of the sidelobes of the radio telescope (which is just another way of saying “pdf”)..

Another case of correlated measurements occurs when your assumed model is incorrect. This is the very definition of correlation, because the residual δy_{m} is correlated with the data value y_{m}. But how do you calculate V? If you could do a large number J of experiments, each with M datapoints producing measured values y_{m,j}, each measured at different values of x_{m}, then each element of the covariance matrix would be V_{mn} = ∑_{j}(y_{m,j} − \overline{y_{m}})(y_{n,j} − \overline{y_{n}}). You don’t normally have this opportunity. Much better is to look at your residuals; if the model doesn’t fit, use another one!

Normally, and in particular we assume everywhere in this tutorial, the measurements are uncorrelated, so one takes V = I (the unitary matrix).

9.5. The case in which datapoints have different dispersions, or different weights: like a weighted average

Here the σ_{meas,m} are all different. The m^{th} row of the equation-of-condition matrix X_{σ} and the m^{th} element of the data vector Y_{σ} get divided by their corresponding σ_{meas,m}. The equation embodied in each row of the matrix equation 9 remains unchanged, but the different rows are weighted differently with respect to each other.

Consider two measurements with intrinsic measurement uncertainties (σ_{meas,1}, σ_{meas,2}); suppose σ_{meas,1} < σ_{meas,2}. After being divided by their respective σ_{meas}’s, all of the numbers in row 1 are larger than those in row 2. In all subsequent matrix operations, these larger numbers contribute more to all of the matrix-element products and sums. Thus, the measurement with smaller uncertainty has more influence on the final result, as it should.

Suppose that the above two measurements were taken under identical conditions except that measurement 1 received more integration time than measurement 2; we have \frac{σ_{meas,1}}{σ_{meas,2}} = (\frac{τ_{1}}{τ_{2}} )^{−1/2}, so the rows of X are weighted as τ^{1/2}. This means that during the computation of [α] = X^{T} · X, the self-products of row 1 are weighted as τ_{1}. This means that each datapoint is weighted as τ, which is exactly what you’d expect! Note that this is also exactly the same weighting scheme used in a weighted average, in which the weights are proportional to (\frac{1}{σ_{meas}} )^{2}. We conclude that weighting scheme of the first two steps in section 9.3 agrees with common sense.

Suppose you don’t know the intrinsic measurement dispersion σ_{meas}, but you do know the relative dispersion of the various measurements. For example, this would be the case if the datapoints were taken under identical conditions except for integration time; then σ_{meas} ∝ τ^{−1/2}. In this case, multiply each row by its weight w ∝ \frac{1}{σ_{meas}} and proceed as above.

9.6. The explicit IDL for the weighted and chi square cases

Here we give the explicit IDL commands required to do chi squared fitting, as we did for least squares fitting in §5.2. First, define the X and Y arrays as before in §5. However, here we need to weight the matrices. Let us denote the weighted versions by XW and YW. Then each row of XW and element of YW are divided by the corresponding σ_{meas}. Let us define W as the diagonal matrix^{5} whose elements are equal to \frac{1}{σ_{meas,m}}, i.e. w_{m} = \frac{1.}{σ_{meas,m}}

W = fltarr(M, M)  (57a)
W[indgen(M) ∗ (M + 1)] = [\frac{1}{σ_{meas,m}} ]  (57b)

5W has M × M elements but only M numbers. If M is large and you are impacting machine memory, then you can save memory by: reformulating equations 58 and 64 to use loops; not keeping XW and X as separate matrices; and using IDL’s temporary function wherever possible. The relevant for loop is for mr= 0, m do xw[*,mr]=w[mr]*x[*,mr], etc.

Note that w_{m} is the square root of the weights w_{meas,m} that we discussed in §9.1. Specifically, the current weights w_{m} are not those in the sense of the weighted average of §9.1; we use w_{m} as a matter of convenience only. Then

XW = W ## X  (58a)
YW = W ## Y  (58b)

All steps up to obtaining the uncertainties follow exactly the same steps as in §5.2, which we repeat here. In IDL we calculate the normal equation matrices and denote the [α] in equation 11a by XXW:

XXW = transpose(XW)##XW ,  (59a)

and we denote the [β] in equation 11b by XYW:

XYW = transpose(XW)##YW .  (59b)

Next we take the inverse of [α] (same as XXW). The standard way is

XXWI = invert(XXW) .  (60)

Optional comment: To take the inverse you might wish to revert to Singular Value Decomposition method; see §14 for a description. In my experience it is rare, but sometimes you will encounter a problem for which the standard inversion of XXW is unstable. In this case, XXW##XXWI ≠ I (I is the unit matrix). You can get a rough indication of this problem by calling the invert function with the status keyword; if status ≠ 0, you’ve definitely got problems^{6}.

The least-square fitted quantities are in the matrix a (equation 12), which we obtain in IDL with

a = XXWI ## XYW .  (61)

In IDL we denote the matrix of weighted predicted values by YBARW, which is

6However, the status keyword of IDL’s invert function is not fully reliable.

YBARW = XW ## a .  (62)

and we can also define the weighted residuals as

DELYW = YW − YBARW  (63)

The true (unweighted) predicted values and residuals are

YBAR = invert(W) ##XW ## a = X ## a .  (64a)

DELY = invert(W) ## (YW − YBARW) = Y − YBAR  (64b)

To calculate the chi square we want the weighted sum ∑ δw_{m}^{2}y_{m}^{2}:

chi sq = transpose(DELYW) ## DELYW .  (65a)

and the reduced chi square

red chi sq = transpose(DELYW) ## DELYW/(M − N) .  (65b)

The full covariance matrix is just XXWI and, corresponding to equation 30, the variances of the derived coefficients are the diagonal elements—unless the reduced chi squared differs significantly from unity, in which case you should seriously look at your assumed weights derived from σ_{meas,m} and/or the quality with which the model fits the data. That is,

vardc official = XXWI[(M + 1) ∗ indgen(M)]  (66a)

or, better,

vardc realworld = red chi sq ∗ XXWI[(M + 1) ∗ indgen(M)]  (66b)

Finally, the normalized covariance matrix is of course independent of chi squared and is the same as equation 33

ncov = XXWI/sqrt(vardc official#vardc official)  (67)

You should, of course, always look at the residuals from the fit. But here you need to look at the weighted residuals DELYW.

10. CHI-SQUARE FITTING, AND WEIGHTED FITTING: DISCUSSION INCLUDING COVARIANCE

10.1. Phenomenological description

Consider the first two coefficients in our example of §6.2. In this example, the fit gives y = A_{0} +A_{1}t+A_{2}t^{2}, where the numerical values are given in vector form by equation 31. The coefficient A_{0} is the y-intercept and A_{1} is the slope. They have derived values A_{0} = 96 ± 34 and A_{1} = 4 ± 9.

Remember what these uncertainties really mean: in an infinity of similar experiments, you’ll obtain an infinity of values of (A_{0}, A_{1}) that are normally distributed with dispersions (34,9). Loosely speaking, this means that A_{0} lies between (96 − 34 = 62) and (96 + 34 = 130) and A_{1} lies between −4.5 and 13.5.

Suppose you are interested in knowing about A_{0} without regard to A_{1}. By this we mean that as A_{0} is varied from its optimum value of 96, χ^{2} increases from its minimum value. As we vary A_{0}, if we allow A_{1} to take on whatever value it needs to for the purpose of minimizing χ^{2}, then this is what we mean by “knowing about A_{0} without regard to A_{1}”. For this case, the uncertainty of A_{0} is indeed 34. Ditto for A_{1}. In other words, equations 19, and 54 apply.

However, if you are interested in knowing about both, you must include their covariance. In our example, the large negative covariance follows logically just from looking at a graph: if you fit some points, all of which lie at positive t, then a more negative derived slope will raise the y-intercept.

Specifically, the large negative covariance means that positive departures of A_{0} are associated with negative departures of A_{1}. So even though the individual values δA_{0} = +34 and δA_{1} = +9 are acceptable, you cannot conclude that the pair of values (δA_{0}, δA_{1}) = (+34, +9) is acceptable, because this pair has both positive. In contrast, what is acceptable here would be something like (δA_{0}, δA_{1}) = (+34, −9).

We stress that the acceptable ranges of values depend on what you are interested in. This is sort of like the observer’s influence in quantum mechanics. If you are interested in A_{1} alone, then you can say A_{1} = 4 ± 9 and, in making this statement, you have to realize that, as A_{1} varies over this range, A_{0} can vary over (formally, at least) the range (∞ → −∞): you just don’t give a damn what happens to A_{0} because you’re not interested. But the moment you become interested and restrict its possible range, that influences the possible range for A_{1}, too.

There is no simple relationship between the covariance matrix elements and the acceptable ranges. For two variables, the best way to express this is to construct the ellipses that define the loci of constant Δχ^{2} and present them on a graph with axes (δa_{0}, δa_{1}) as in BR Figure 11.2 or NR Figure 14.5.4. For three variables, these ellipses become ellipsoids; for four, they become four-dimensional volumes, etc.

Fig. 2.— Illustrating the meaning of variance and covariance between (a1, a2) for our numerical example. See text for discussion.

We illustrate these concepts for the (a_{1}, a_{2}) parameters for our numerical example. We subtracted 7.75 from all times so that the covariance would be small enough to illustrate the difference between the tangents to the ellipses and the end points of the ellipses. Contours are calculated as described in §10.5 and are at Δχ^{2} = 1 and 2.3. The dashed horizontal and vertical lines are at δ_{a} = ±σ_{a}.

First consider the pair of vertical lines, which are drawn at δa_{1} = ±σ_{a_{1}}, where σ is the square root of the variance of the parameters as described in equations 19, 30, 54, 55, and 66. If the datapoints were projected downward, i.e. look at the marginal pdf of δa1 by taking small strips δa_{1} and integrating over δa_{2}, the marginal pdf of δa_{1} is Gaussian; ditto for the other coordinate.

Thus, 68% of the points lie between these dashed lines. This is what we mean by the phrase “being interested in knowing about a_{1} without regard to a_{2}”. If we allow a_{2} to vary so as to minimize χ^{2} as we consider departures δa_{1}, then the pdf of δa_{1} has dispersion σ_{a_{1}}. Alternatively, we can say that in a large number of experiments, the pdf of δa_{1} follows a chi squared pdf with one degree of freedom if we don’t care what happens to δa_{2}.

If, however, we are concerned about the pair, then we must look not at the projection down one axis or the other, but rather at the two-dimensional distribution. This is characterized by the tilted ellipses. Here, for a large number of experiments, the pair (a_{1}, a_{2}) follows a chi square distribution with 2 degrees of freedom (if we don’t care about a_{0}; if we do, it’s 3 degrees of freedom and the ellipse becomes an ellipsoid, but this is very hard to plot!). For ν = 2, 68.3% of the points lie within Δχ^{2} = 2.3, where we have drawn the outer contour in Figure 2. The points inside this ellipse are darker; 68.3% of the points lie within that ellipse.

The best description of the specifics of calculating these ellipsoids is in BR §11.5 (Confidence Intervals, Confidence Levels for Multiparameter Fits). To describe it, we’ll talk specifically about our numerical example, which has M = 4 measurements and N = 3 unknowns. The unknowns are a = [a_{0}, a_{1}, a_{2}]. We’ll first begin by discussing the case of single parameter; then we’ll generalize.

10.2. Calculating the uncertainties of a single parameter—gedanken experiment

First, suppose we want to know the value σ_{a_{0}} without regard to the values of a_{1} and a_{2}. Having already done the solution, we know the chi-square value of a_{0} so we consider variations δa_{0} around this best value.

Pick a value of δa_{0} and redo the least squares solution for [a_{1}, a_{2}]. This gives a new value for χ^{2} which is, of course, larger than the minimum value that was obtained with δa_{0} = 0. Call this difference Δχ_{δa_{0}}^{2}. Determine the dependence of Δχ_{δa_{0}}^{2} upon δa_{0} and find the value of δa_{0} such that Δχ_{δa_{0}}^{2} = 1. This is the desired result, namely the value σ_{a_{0}} without regard to the values of a_{1} and a_{2}.

This value is σ_{a_{0}}^{2} = [α]_{00}^{−1}, the same result quoted in equation 55.

Consider now what you’ve done in this process. For each least squares fit you used a trial value of δa_{0}. In specifying δa_{0} you had exactly one degree of freedom because you are fixing one and only one parameter. Having done this, you could do a large number of experiments (or Monte Carlo trials) to determine the resultant distribution of Δχ_{δa_{0}}^{2}. It should be clear that this distribution follows a chi squared distribution with one degree of freedom (ν = 1). So the uncertainty σ_{a_{0}} is that value for which Δχ_{δa_{0}}^{2} = 1. (The chi square fit for the other two parameters has M −2 degrees of freedom, but this is irrelevant because—by hypothesis—you don’t care what happens to those variables.)

10.3. Calculating the uncertainties of two parametersr—gedanken experiment

Suppose we want to know the value (σ_{a_{0}}, σ_{a_{2}}) without regard to the value of a_{1}. Now we consider variations (δa_{0}, δa_{2}) around the best values (a_{0}, a_{2}).

Pick values for (δa_{0}, δa_{2}) and redo the least squares solution for a_{1}. This gives a new value for χ^{2} which is, of course, larger than the minimum value that was obtained with (δa_{0}, δa_{2}) = 0. Call this difference Δχ_{(δa_{0},δa_{2})}^{2}. As above, this follows a chi square distribution, but now with ν = 2. Determine the dependence of Δχ_{(δa_{0},δa_{2})}^{2} upon (δa_{0}, δa_{2}) and find the set of values of (δa_{0}, δa_{2}) such that Δχ_{(δa_{0},δa_{2})}^{2} = 2.3. This is the desired result, namely the ellipse within which the actual values (δa_{0}, δa_{2}) lie with a probability of 68.3%, without regard to the value of a_{1}.

These values can be defined in terms of the curvature matrix [α], as we discuss below.

Consider now what you’ve done in this process. For each least squares fit you used trial values of (δa_{0}, δa_{2}). In specifying them you had exactly two degrees of freedom because you are fixing two parameters. This distribution follows a chi squared distribution with two degree of freedom (ν = 2). So the uncertainty σ_{a_{0}} is that value for which Δχ_{δa_{0}}^{2} = 2.3, which follows from the integrated probability for the chi squared distribution for ν = 2. (The chi square fit for the third parameter a_{1} has M − 1 degrees of freedom, but again this is irrelevant.)

One can expand this discussion in the obvious way. Consider finally…

10.4. Calculating the uncertainties of three parametersr—gedanken experiment

Suppose we want to know the values of all three parameters (or, generally, all N parameters). Then we pick trial values for all three. There is no least squares fit for the remaining parameters, because there are none. For each combination of the three (or N) parameters we obtain Δχ_{a}, which defines a 3- (or N-) dimensional ellipsoid. This follows a chi squared distribution with ν = 3 (or N). We find the (hyper)surface such that Δχ_{a} is that value within which the integrated probability is 68.3%. This defines the (hyper)surface of σ_{a}.

10.5. Doing these calculations the non-gedanken easy way

The obvious way to do the gedanken calculations described above is to set up a grid of values in the parameters of interest (δa_{n}); perform the chi squared fit on the remaining variables, keeping track of the resulting grid of χ^{2}; and plot the results in terms of a contour plot (for two parameters of interest) or higher dimensions.

There’s an easier way which is applicable unless you are doing a nonlinear fit and the parameter errors are large^{7}. The curvature matrix ([α] of equation 52d, the same as XXW) contains the matrix of the second derivatives of χ^{2} with respect to all pairwise combinations of δa_{n}, evaluated at the minimum χ^{2}; it’s known as the curvature matrix for this reason. Clearly, as long as the Taylor expansion is good we can write

Δχ_{a}^{2} = δa^{T} · [α] · δa  (68)

Knowing the curvature matrix, we don’t have to re-do the fits as we described above. Rather, we can use the already-known matrix elements. Use the following procedure [see NR §15.6 (Probability Distribution of Parameters in the Normal Case)]. Follow the steps

1. Decide which set of parameters you are interested in; call this number N_{i} and denote their vector by a_{i}. Here we use the above example and consider N_{i} = 2 and a_{i} = [a_{0}, a_{2}].

2. From the N × N covariance matrix [α]^{−1}, extract the rows and columns corresponding to the N_{i} parameters and form a new N_{i} × N_{i} covariance matrix [α]_{i}^{−1}; in our case the original covariance matrix is

[α]^{−1} = XXI = [ \begin{matrix}1156.8125 & −303.000 & 18.4375 \\ −303.000 & 81.000 & −5.000 \\ 18.4375 & −5.000 & 0.31250 \end{matrix} ]  (69a)

and it becomes

[α]_{i}^{−1} = XXI = [ \begin{matrix}1156.8125 & 18.4375 \\ 18.4375 & 0.31250 \end{matrix} ] .  (69b)

3. Invert this new covariance matrix to form a new curvature matrix [α]_{i}. The elements differ from the those in the original curvature matrix.

4. As usual, we have

Δχ_{a_{i}}^{2} = δa_{i}^{T} · [α]_{i} · δa_{i}  (70)

so find the locus of a_{i} such that the integrated probability of Δχ_{a_{i}} for ν = N_{i} contains 68.3% of the space; e.g. for ν = 2 this is Δχ_{a_{i}} = 2.3.

7In which case you use the gedanken technique!

You may well wonder why, in steps 2 and 3, you need to derive a new curvature matrix from the extracted elements of the original covariance matrix. Why not just use the extracted elements of the original curvature matrix? To understand this, read the NR’s discussion surrounding equation (15.6.2); this is not very clear, in my opinion. I find it easier to recall the way covariance matrices propagate errors according to the derivatives of the quantities derived, as in Cowan’s equation 1.54. I could explain this here but, quite frankly, don’t have the time; maybe in the next edition of these notes!

10.6. Important comments about uncertainties

Having said all the above, we offer the following important Comments:

• The easiest way to calculate these (hyper)surfaces is to set up a grid in N_{i}-dimensional space of trial values for δa_{i} and use a contour plot or volume plot package to plot the loci of constant Δχ_{a_{i}}^{2}.

• The procedure described in §10.5 works well for linear fits, or nonlinear fits in which the σ_{a} are small so that Δχ^{2} is well-approximated by the second derivative curvature matrix. This is not necessarily the case; an example is shown in BR Figure 11.2. Here, the higher-order curvature terms are important and it’s better to actually redo the fit for the grid of trial values of a_{i} as described above in §10.2, 10.3, and 10.4. In other words, use the gedanken technique.

• The variance (i.e., uncertainty squared) of the derived parameters a depends only on the elements in the covariance matrix [α]^{−1}. These, in turn, depend only on the curvature matrix [α]. These, in turn, depend only on X. This matrix X is the matrix of the quantities that are known exactly. For example, we began with the example in which the elements of X were the times at which the measurements were taken. Generally, then, the curvature and covariance matrix elements depend on the locations of the datapoints (the ensemble of t_{m} in equation 1) but not on the measured values (the ensemble of y_{m} in equation 1). And on your adopted values for σ_{meas,m}. Because of this…

• Think before making your measurements about the covariance matrix and how to minimize the off-diagonal elements. By taking measurements at well-chosen times, or well-chosen values of the independent variable x_{m} whatever it is, you can really optimize the accuracy-to-effort ratio! For example, in our numerical example if you can get a few measurements at negative times you’ll efforts will be repaid in terms of much better accuracy for the y-intercept.

11. REJECTING BAD DATAPOINTS II.: STETSON’S METHOD PLUS CHAUVENET’S CRITERION

Chauvenet’s criterion is an on-off deal: either you include the datapoint or you don’t. This makes sense from a philosophical point of view: either a datapoint is good or not, so you should either include it or exclude it.

However, when doing a nonlinear fit this presents a problem. As you iterate, the solution changes, and a given datapoint can change from being “bad” to “good”. Or vice-versa. You can imagine being in a situation in which the iteration oscillates between two solutions, one including a particular datapoint and the other excluding it; the solution never converges, it just keeps chasing its tail.

Enter Stetson’s beautiful technique^{8}. Stetson reasons that we shouldn’t have an on-off criterion. Rather, it should relieve a datapoint of its influence adiabatically: as its residual gets larger and larger, its weight gets smaller and smaller. With this, in nonlinear fitting all datapoints are always included and their weights automatically adjust as the fit parameters home into their correct values. And you can’t get into the chasing-tail syndrome that can happen with the strict on-off inclusion.

11.1. Stetson’s sliding weight

Stetson recommends using a sliding weight. To explain this, we review the ML concept of chi square fitting. We write equations for the case of a single parameter to keep things simple, but it works for multiparameter fits too. In chi square fitting for the single parameter a, we define chi square as

χ^{2} = ∑_{m=0}^{M−1} \frac{(y_{m} − af_{(x_{m})})^{2}}{σ_{m}^{2}}  (71a)

and we minimize χ^{2} by setting its derivative with respect to each parameter a equal to zero:

\frac{dχ^{2}}{da} = −2 ∑_{m=0}^{M−1} \frac{f_{(x_{m})}Δy_{m}}{σ_{m}^{2}}  (71b)

Here Δy_{m} = (y_{m} − af_{(x_{m})}). This gives

8Stetson is one of those anomalies, a true expert on fitting. He invented many of the stellar photometry routines used in daophote, all of which use least squares techniques. He provides a lively, engaging discussion of many fascinating and instructive aspects in his website: http://nedwww.ipac.caltech.edu/level15/Stetson/Stetson4.html.

∑_{m=0}^{M−1} \frac{f_{(x)}Δy_{m}}{σ_{m}^{2}} = 0  (71c)

Now we wish to modify this equation by introducing a weight w_{(|Δy_{m}|)} that makes datapoints with large |Δy_{m}| contribute less, so it reads like this:

∑_{m=0}^{M−1} \frac{w_{(|Δy_{m}|)}f_{(x_{m})}Δy_{m}}{σ_{m}^{2}} = 0  (72)

It’s clear that we need the following properties for w_{(Δy_{m})}:

1. w_{(Δy_{m})} = w_{(|Δy_{m}|)}, meaning simply that it should depend on the absolute value of the residual and not bias the solution one way or the other.

2. w_{(|Δy_{m}|)} → 1 as |Δy_{m}| → 0, meaning that datapoints with small residuals contribute their full weight.

3. w_{(|Δy_{m}|)} → 0 as |Δy_{m}| → ∞, meaning that datapoints with large residuals contribute nothing.

Stetson recommends

w_{(|Δy_{m}|)} = \frac{1}{1 + (\frac{|Δy|}{ασ})^{β}}  (73)

This function w_{(|Δy_{m}|)} has the desired properties. Also, for all β it equals 0.5 for |Δy_{m}| = ασ. As β → ∞ the cutoff gets steeper and steeper, so in this limit it becomes equivalent to a complete cutoff for |Δy_{m}| > ασ.

Stetson recommends α = 2 to 2.5, β = 2 to 4 on the basis of years of experience. Stetson is a true expert and we should take his advice seriously; he provides a vibrant discussion to justify these choices in real life, including an interesting set of numerical experiments.

However, for large M I see a problem with the choice α = 2 to 2.5. For large β, for which the cutoff is sharp, it seems to me that the cutoff should duplicate Chauvenet’s criterion. Referring to equation 42, this occurs by setting

α = erf^{−1} (1 − \frac{1}{2M} )  (74)

and I recommend making this change, at least for problems having reasonably large M; this makes α larger than Stetson’s choice. I’m more of a purist than Stetson, probably because I’m a radio astronomer and often fit thousands of spectral data datapoints that are, indeed, characterized mainly by Gaussian statistics. Stetson is an optical astronomer and probably sees a lot more departures from things like cosmic rays. Nevertheless, in a CCD image with millions of pixels, of which only a fraction are characterized by nonGaussian problems such as cosmic ray hits, it seems to me only reasonable to increase α above Stetson’s recommended values by using equation 74.

11.2. Implementation of the weight in our matrix equations

Clearly, implementing Stetson’s method requires a weighted fit, so you have to use the chi squared technique discussed in §9.6. There equation 57 defines a matrix of weights (which is diagonal) in which

W_{m,m} = \frac{1}{σ_{m}}  (75)

Comparing this with equation 71c, it’s clear what to do: we modify this equation to read

W_{m,m} = \frac{w_{m}^{1/2}}{σ_{m}}  (76)

where the weight w_{m} is defined in equation 73.

Now you must not forget here that the solution depends on the weights w_{m}, which in turn depend on the solution. Thus when you implement this technique you must iterate until the solution converges by not changing.

12. MEDIAN, INSTEAD OF LEAST SQUARE, FITTING

Least square fitting minimizes the square of the residuals. This means that datapoints having large residuals contribute importantly to the fit. If these datapoints are really bad, you’ve got problems; this is why it’s important to get rid of outliers! Sometimes you’re faced with a set of datapoints that look bimodal: most datapoints have a Gaussian-like pdf, and many lie outside the main distribution; sometimes it’s difficult to decide “where to draw the line”. Or you might have nonGaussian statistics. In these cases, using least squares is not a great idea because least squares gives greatest weight to the datapoints having the largest residuals, but you don’t know what the residuals are until after you’ve done the fit—and the fit is influenced, and maybe even dominated, by the outliers!

In these cases most astronomers use the median^{9}. The median is the solution for which there as many positive as negative residuals, irrespective of how big they are. As long as the discrepant datapoints are symmetrically distributed with respect to sign, the median works well. In fact, for any pdf that is symmetrically distributed with respect to sign, the median is a good solution because it doesn’t depend at all on the details of the pdf. If the statistics are nonGaussian then the error of the median is somewhat greater than that of the mean (by a factor of something like \frac{π}{2} ; I forget, but it’s straightforward to calculate) but, when the datapoints are contaminated by nonGaussian statistics, this is often less than the error introduced by the least squares technique.

12.1. The median and the double-sided exponential pdf

In fact, there is a specific pdf for which the median is the theoretically correct solution: the double-sided exponential. Here the pdf of the measured datapoints is

p_{(Δy_{m})} = \frac{e^{−|Δy_{m}|/σ_{m}}}{2σ_{m}}  (77)

where, again, Δy_{m} = (y_{m} −af_{(x_{m})}). For this, the logarithm of the likelihood function is (we exclude the term involving log Π_{m=0}^{M−1} \frac{1}{σ_{m}} for simplicity)

L_{(Δy_{m})} = log(L_{(Δy_{m})}) = − ∑_{m=0}^{M−1} [\frac{|y_{m} − af_{(x_{m})})|}{σ_{m}}]  (78a)

The absolute value signs are horrible to deal with, so we rewrite this as

L_{(Δy_{m})}) = ∑_{Δy_{m}>0} \frac{y_{m} − af_{(x_{m})}}{σ_{m}} − ∑_{Δy_{m}<0} \frac{y_{m} − af_{(x_{m})}}{σ_{m}}  (78b)

Now we take the derivative of L with respect to a and set it equal to zero to find the maximum. This gives

\frac{dL}{da} = ∑_{Δy>0} \frac{f_{(x_{m})}}{σ_{m}} − ∑_{Δy<0} \frac{f_{(x_{m})}}{σ_{m}} = 0  (79)

Amazing—neither Δy_{m} nor a appears!

9Stetson recommends, instead of the median, his sliding weight technique. You should read his website

Consider first the case f_{(x_{m})} = 1; this is like an average. If all σ_{m} are identical, this says: the proper solution is the median, for which half the datapoints have positive residuals! More generally, we have to define a generalized median where each datapoint is weighted by its \frac{f_{(xm)}}{σ_{m}}. Can you show that the median of the residuals is zero?

12.2. Doing a “chi square” fit for the double-sided exponential

It’s not obvious how to take a weighted median! Fortunately, we can use our standard chi square technique to do a median fit by adjusting the weights. To see how, we rewrite equation 72 here; this is the equation that has Stetson’s additional weight w_{(|Δy_{m}|)}

∑_{m=0}^{M−1} \frac{w_{(|Δy_{m}|)}f_{(x_{m})}Δy_{m}}{σ_{m}^{2}} = 0  (80)

We need to devise weight w_{(|Δy_{m}|)} so that this equation looks like equation 79. By inspection, we require

w_{m} = \frac{σ_{m}}{|Δy_{m}|}  (81)

so that in equations 72 and 80 we have

\frac{w_{m}}{σ_{m}^{2}} = \frac{1}{σ_{m}|Δy_{m}|}  (82)

and in weight matrix W of equations 57 and 76 we use the square root of the above quantity.

This works for multivariate fits using any set of functions f_{(x_{m})}. However, we caution again that the weights w_{m} depend on the parameters a and the parameters depend on the weights, so you must iterate. In the experiments I’ve performed the iteration is slow (but sure). And at the end you will certainly find that the reduced chi square χ^{2} will be nowhere near unity; nevertheless, you should use errors calculated χ^{2}, not forcing χ^{2} to equal unity (so use equation 55 instead of 54). [THIS STATEMENT IS ALMOST CERTAINLY CORRECT BUT NEEDS TO BE CHECKED WITH NUMERICAL EXPERIMENTS…].

12.3. IDL’s resources for median fitting

There are two cases for which IDL provides median fitting. One is IDL’s median function, which takes the median of a bunch of datapoints. It’s great if there’s no functional dependence on x_{m}. The other is ladfit (“least absolute deviation fit”). These are quick and you don’t have to do any iterations. And they give a slightly more accurate answer than the above technique; this is because IDL is careful about the exact definition of the median, e.g. when there is an even number of datapoints.

However, with IDL’s routines there are no errors given; there is no possibility for differing σ_{m}; there is no covariance matrix. And, of course, you can’t fit anything more complicated than a straight line.

13. FITTING WHEN ALL VARIABLES HAVE UNCERTAINTIES

We’ve mentioned that one of the essential assumptions of least squares is that the independent variables are known with high precision and the errors occur only in the measured data. Suppose you’re fitting two variables, t and y, as in equation 1. This essential assumption means that t is known with high precision and all the uncertainty is in y, and you are minimizing the squares of the residuals in the y-direction only. If both variables have uncertainties, then you have to be careful because the essential assumption is violated. If you go ahead with a standard least squares fit when there are errors in both coordinates, the slope will be systematically too small, as we discuss briefly below.

Thanks to Stetson, the ML formulation of this problem is straightforward. Nevertheless, as far as I know the proper formulation is rather recent; Stetson’s treatment is the first I’ve seen. Before reviewing Stetson’s formulation, let’s see two approaches:

1. Taylor §8.4 argues that you can account for x-variance σ_{x_{m}}^{2} by increasing the y-variance by the usual error propagation, i.e. define an equivalent y-variance σ_{y_{m}}^{2}(equiv) = [σ_{y_{m}}^{2} + (a_{1}σ_{x_{m}})^{2}], where a_{1} is the slope. This is equivalent to our results below.

2. Isobe et al (1990, ApJ 364, 104) discuss the case incorrectly. Look in particular at their Section V, where they make 5 numbered recommendations. Two of these are incorrect:

(a) Number 3 says, in essence, that if you have measurement errors in y but not in x, and want to predict x from y in some future dataset, that you should least squares fit the x values (which have no errors) to the y. This is flat wrong. Again, it leads to a slope that is systematically too small. The proper procedure is to fit y to x in the standard way, which is consistent with the ML formulation and gives the right answer; then use the resulting parameters, whose errors you know about, to predict x from y in the future.

(b) Number 4 says that if both x and y have errors, and your main focus is finding the true slope, you should use their “bisector” method. I won’t explain this because this concept is wrong.

Stetson has a beautiful discussion of the problem, which is obviously correct—both because it makes sense and because numerical experiments show that the derived parameters don’t depend on whether you consider x or y the independent variable. With his technique you are not restricted to fitting a straight line (a one-degree polynomial). NR §15.3 provides another method, which is more complicated and is restricted to a straight line. Here we review Stetson’s method.

13.1. A preliminary: Why the slope is systematically small

Why is the derived slope systematically too small if you use the standard least squares technique when both variables have errors? To see this, take a look back at equation 5, where we explicitly write the normal equations for fitting a straight line of the form As_{m} + Bt_{m} = y_{m}. To focus the discussion and make it easy, replace that problem with a single-parameter solution for only the slope B, and use the usual variables (x, y) in place of (t, y). Then we are fitting the set of M equations

Bx_{m} = y_{m}  (83a)

the set of two normal equations becomes just the single equation

B[x^{2}] = [xy]  (83b)

or, writing out the sums explicitly,

B = \frac{∑_{m=0}^{M−1} x∗_{m}y_{m}}{∑_{m=0}^{M−1} x∗_{m}^{2}}  (83c)

Here we use the star to designate the perfectly-known independent variable x∗_{m}. It is important to realize that the x_{m} that appear in this equation are the perfectly-known ones x∗_{m}; this is a fundamental tenet of least squares fitting, which comes from the concept and principle of maximum likelihood ML.

Because B is defined by the x∗_{m} and we are asking what happens when we use the imperfectly known x_{m} instead, let us reduce the problem to the essence and imagine that y_{m} is perfectly known, i.e. y_{m} = y∗_{m}; and that

x∗_{m} = x_{m} − δx_{m}  (84)

where δx_{m} is the observational error in point m. If we do standard least squares on this situation, then we (incorrectly) rewrite equation 83c to read

B = \frac{∑_{m=0}^{M−1} x_{m}y_{m}}{∑_{m=0}^{M−1} x_{m}^{2}}  (85a)

that is, using x_{m} instead of x∗_{m} (because we don’t know what x∗_{m} is!). Substituting equation 84, and remembering that y_{m} = y∗_{m} = Bx∗_{m}, we have

B = \frac{∑_{m=0}^{M−1} x∗_{m}(x∗_{m} + δx_{m})}{∑_{m=0}^{M−1}(x∗_{m}^{2} + 2x∗_{m}δx_{m} + δx_{m}^{2})}  (85b)

Now all terms having δx_{m} sum to zero because the errors are distributed symmetrically around zero. But the denominator contains δx_{m}^{2}. The denominator is irrevocably increased by this term, which decreases the derived value of B from its true value. Yes, this is only a second-order effect, but it matters—after all, χ^{2} is a second-order quantity! Try some numerical experiments!

13.2. Stetson’s method

13.2.1. Philosophy

To make things clear, we discuss the case of a straight line. Thus we assume

y∗ = a_{0} + a_{1}x∗  (86)

where the ∗ indicates the true quantity without errors. The measured points differ from the true ones:

δx_{m} = x_{m} − x∗_{m} ; δy_{m} = y_{m} − y∗_{m}  (87)

and we define, as usual, the apparent residual in y_{m} as

Δy_{m} = y_{m} − a_{1}x_{m}  (88)

Note that this differs from the true residual in y_{m}, which is δy_{m}; the apparent residual Δy_{m} assumes that x_{m} has no error.

Look at Figure 3. There (x_{m}, y_{m}) (the square plot symbol) is a measured point. The ellipse centered on it has (x, y) axial ratio (σ_{x}, σ_{y}) and traces out a the locus where the point’s contribution to χ^{2} is constant; call this contribution χ_{m}^{2}. The straight line is the fitted straight line, assumed to be perfectly correct, i.e. the line given by equation 86. The most probable place where the observed Fig. 3.— Illustrating the discussion of fitting when both variables have errors. Here, as usual in this document, differences are “measured minus true”, e.g. Δy_{m} = y_{m} − y∗_{m}, etc.; in the figure, Δy_{m} < 0, δy_{m} < 0, δx_{m} > 0.

point (x_{m}, y_{m}) would lie if it had no errors is (x∗_{m}, y∗_{m}); it lies on the smallest ellipse centered on (x_{m}, y_{m}), as pictured. This is the ellipse having the smallest possible χ_{m}^{2}. The true differences between the observed and most probable points are (δx_{m}, δy_{m}) given by equation 87 above, where these are measured from the center to the starred point on the ellipse.

It’s easy to calculate (δx_{m}, δy_{m}) in terms of Δy_{m}. These are the distances from the measured to the starred point. We realize that the starred point must satisfy two criteria:

1. The straight line fit is tangent to the ellipse, so the slopes much match. The ellipse has equation

\frac{δx_{m}^{2}}{σ_{x_{m}}^{2}} + \frac{δy_{m}^{2}}{σ_{y_{m}}^{2}} = const = χ_{m}^{2}  (89)

so the slope is

\frac{dδy_{m}}{dδx_{m}} = −\frac{σ_{y_{m}}^{2}}{σ_{x_{m}}^{2}} \frac{δx_{m}}{δy_{m}}  (90)

and this, evaluated at (x∗_{m}, y∗_{m}) must equal a_{1}.

2. The actual values must match. Rather than solve for these, we can solve in terms of Δy_{m} directly.

From the figure it’s clear that

Δy_{m} = δy_{m} − a_{1}δx_{m}  (91)

so combining equations 90 and 91 we get

δx_{m} = \frac{−\frac{Δy_{m}}{a_{1}}}{1 + \frac{σ_{ym}^{2}}{a_{1}^{2}σ_{xm}^{2}}}  (92a)

and

δy_{m} = \frac{Δy_{m}}{1 + \frac{a_{1}^{2}σ_{xm}^{2}}{σ_{ym}^{2}}}  (92b)

You could minimize χ^{2} by brute force. You would begin with the standard least squares solution to get an estimate of the correct parameter values. Then you’d make up a grid in the (a_{0}, a_{1}) plane to find the location of the minimum in χ^{2}, which is

χ^{2} = ∑_{m=0}^{M−1} χ_{m}^{2} = ∑_{m=0}^{M−1} \frac{δx_{m}^{2}}{σ_{x_{m}}^{2}} + \frac{δy_{m}^{2}}{σ_{y_{m}}^{2}}  (93)

which gives the parameters (a_{0}, a_{1}). The quantities (δx_{m}^{2}, δy_{m}^{2}) are given by equation 92. The errors are given by the usual Δχ^{2} contours as discussed in §10.

13.2.2. Direct solution

Using the principle of maximum likelihood (ME), we find the minimum χ^{2} by taking the derivative of χ^{2} with respect to each parameter a_{n} and setting the derivative equal to zero. To illustrate, we work this through for the slope. We are fitting for the slope in the equation

y_{m} = ax_{m}  (94)

for which χ^{2} is given by equation 93. Express (δx_{m}, δy_{m}) in terms of Δy_{m}, as in equations 92, and rewrite equation 93 to read

χ^{2} = ∑_{m=0}^{M−1} \frac{(y_{m} − ax_{m})^{2}}{a^{2}σ_{x_{m}}^{2} + σ_{y_{m}}^{2}} = ∑_{m=0}^{M−1} \frac{Δy_{m}^{2}}{a^{2}σ_{x_{m}}^{2} + σ_{y_{m}}^{2}}  (95)

Differentiate with respect to a and set equal to zero. After cancelling extraneous factors, and using equation 92a, we get

∑_{m=0}^{M−1} \frac{(y_{m} − ax_{m})(x_{m} − δx_{m})}{a^{2}σ_{x_{m}}^{2} + σ_{y_{m}}^{2}} = 0  (96)

The term ax_{m}(x_{m} −δx_{m}) defines the relevant diagonal matrix element in the curvature matrix, and corresponds to the lower right element (B[t^{2}]) in equation 5; the term y_{m}(x_{m}−δx_{m}) corresponds to [ty] in that equation. In other words, if we had standard least squares and no errors in x, equation 96 would read

∑_{m=0}^{M−1} \frac{(y_{m} − ax_{m})(x_{m})}{σ_{y_{m}}^{2}} = 0  (97)

Looking at these and, also, referring back to §3, it’s easy to see how to generalize this and write it in matrix notation. First, note that each datapoint has a net total error given by the denominator of equation 95, and in general—i.e., when we fit something other than a straight line—the slope a varies from point to point so that, even if the (σ_{x}, σ_{y}) are constant, the net total error is not. This means you need to use the weighted fit technique of §9.3 and 9.6 with

σ_{meas,m}^{2} = (a_{m}^{2}σ_{x_{m}}^{2} + σ_{y_{m}}^{2})  (98a)

where the local slope a_{m} is evaluated at (x_{m} − δx_{m})

a_{m} = \frac{∂y}{∂x}|_{(x_{m}−δx_{m}),(y_{m}−δy_{m})}  (98b)

We repeat all this to be perfectly explicit:

1. You are fitting

X · a = Y  (99)

In X, each column n contains a different function of x called f_{n}(x); each row contains a different value of x, which is x_{m}. ^{10}

2. In our normal equations we have the factor (x_{m} − δx_{m}) instead of x_{m}, which means that in forming the matrix normal equations we need to define a new matrix

X MOD = X − δX  (100)

where the elements of δX are evaluated using equation 92. That is, δX_{m,n} = f_{n}(x_{m} − δx_{m}). If it’s a straight-line fit, the slope is the same everywhere. If it’s a more general fit, the slope is the local slope derived using equation 98; this slope can be calculated either analytically or numerically.

3. Define the weighted versions of X, X MOD, and Y

XW = W ## X = W · X  (101a)
X MODW = W ## X MOD = W · X MOD  (101b)
YW = W ## Y = W · Y  (101c)

where W_{m,m} = \frac{1}{σ_{meas,m}} (from equation 98a).

4. Then the matrix normal equation is

[α mod]·a = [β mod]  (102)

where

[α mod] = X MODW^{T} · X  (103a)

[β mod] = X MODW^{T} · Y  (103b)

10In mathematics, the entry for (row, col) = (m, n) = f_{n}(x_{m}). In IDL, it’s reversed: the entry for (col, row) = (n, m) = f_{n}(x_{m}).

But, you ask, how can we form the matrix δX if we don’t know the slope a? The answer: iterate. That is,

1. Solve y_{m} = a_{0} + a_{1}x_{m} using standard least squares or chi square; or use any other technique to get an initial guess for the parameters a. Calculate the local slopes using (x_{m}, y_{m}).

2. Determine the Δy_{m} and calculate δx_{m} using equation 92a.

3. Calculate the matrices X MODW, [α mod], and [β mod] as above.

4. Solve for new parameters a.

5. Calculate the local slopes using (x_{m} − δx_{m}, y_{m} − δy_{m}).

6. Loop back and repeat steps 2, 3, 4,and 5 until convergence.

7. After convergence, calculate the uncertainties, as usual, by the equivalent of equation 54, i.e.

s_{a}^{2} = χ^{2}diag{[α mod]}  (104)

13.3. Commentary on Stetson’s solution

13.3.1. This solution is general and makes sense

This solution makes sense. Consider standard least squares in which σ_{x_{m}} → 0; here we should revert to standard least squares. Also, when the slope is zero it’s like taking an average of y_{m}, so it shouldn’t matter what σ_{x} is. Equation 92a satisfies these expectations because the denominator → ∞: the ellipse becomes tall and thin so δx_{m} → 0.

Similarly, if σ_{y_{m}} → 0 or a → ∞ then the ellipse becomes very short and fat: δx_{m} → \frac{Δy_{m}}{a_{1}} , meaning that all of the correction goes to δx_{m}. This is equivalent to y being the independent variable, and you can solve it this way, but of course it is much more efficient to use standard least squares with y actually being the independent variable.

Finally, there’s no reason why we had to use a straight-line fit. It could have been any function; in an arbitrary case, we replace the slope by the local slope, namely the derivative of the function evaluated at the current value x_{m} + δx_{m}. So this technique is general. In his website, Stetson provides an example of fitting a circle to points in (x, y)!

13.3.2. Yet another point I do not understand

First, I found Stetson’s website discussion charming and stimulating. However, I found it very difficult to glean just exactly how he suggests solving the problem! His formulation involves using differences from one iteration to the next and differs from the one I’ve given above.

Second, Stetson spends several pages discussing the importance of the difference between two quantities. One is Δy_{m}. The other is

ϵ = −δF_{m} = \frac{∂F}{∂x} |_{x∗,y∗} δx_{m} + \frac{∂F}{∂y} |_{x∗,y∗} δy_{m}  (105)

where

F (x, y) = a_{0} + a_{1}x − y  (106)

You could put contours of constant F on Figure 3. The zero contour is where F (x∗, y∗) = 0, which is what you get with (δx_{m}, δy_{m}) = (0, 0). This zero-level contour is the straight line solution of the problem. Moving away from zero gives contours of constant ϵ; these are offset from the zero contour and have the same shape.

Stetson emphasizes that the proper thing to do is use ϵ_{m} in place of Δy_{m} in equation 92a. However, after scratching my head for quite a while, I conclude that in fact these two quantities are always equal, i.e. ϵ_{m} = Δy_{m}. I do not understand Stetson’s emphasis on this point.

13.4. Generalization to Multivariate Case

The subsections above deal with the simple case of two variables. It is straightforward to generalize to the multivariate case. As in section 13.2.1, suppose you are fitting the multivariate analogy of equation 86

z∗ = a_{0} + a_{1}x∗ + a_{2}y∗ + …  (107a)

which is equivalent to

dz∗ = a_{1}dx∗ + a_{2}dy∗ + …  (107b)

Then the discussion follows the same track except that here, instead of the surface of χ^{2} being an ellipse on a two-dimensional space, it’s an ellipsoid in a three-dimensinal space. And for the case of N variables it’s an N-dimensional ellipsoid. We can’t draw this, but the ideas embodied in Figure 3 remain: we match the slopes and values to find the smallest ellipse, i.e. the smallest χ^{2}, for each measured point. Equations 90 and 91 become

dδz_{m} = −\frac{σ_{z_{m}}^{2}}{σ_{x_{m}}^{2}} \frac{δx_{m}}{δz_{m}} dδx_{m} − \frac{σ_{z_{m}}^{2}}{σ_{y_{m}}^{2}} \frac{δy_{m}}{δz_{m}} dδy_{m} + …  (108)

and

Δz_{m} = δz_{m} − a_{1}δx_{m} − a_{2}δy_{m} + …  (109)

Comparing equations 107b and 109, we see that

a_{1} = −\frac{σ_{z_{m}}^{2}}{σ_{x_{m}}^{2}} \frac{δx_{m}}{δz_{m}} → δz_{m} = −\frac{σ_{z_{m}}^{2}}{a_{1}σ_{x_{m}}^{2}} δx_{m}  (110a)
a_{2} = −\frac{σ_{z_{m}}^{2}}{σ_{y_{m}}^{2}} \frac{δy_{m}}{δz_{m}} → δz_{m} = −\frac{σ_{z_{m}}^{2}}{a_{2}σ_{y_{m}}^{2}} δy_{m}  (110b)
⋮  (110c)

There are N such equations. We then generate N equations that express Δz_{m} in terms of δx_{m}, δy_{m}, … by successively substituting for δz_{m} the right-hand-side versions of equations 110 into equation 109:

Δz_{m} = −a_{1} [\frac{σ_{z_{m}}^{2}}{(a_{1}σ_{x_{m}})^{2}} + 1] δx_{m} − a_{2}δy_{m} − …  (111a)
Δz_{m} = −a_{1}δx_{m} − a_{2} [\frac{σ_{z_{m}}^{2}}{(a_{2}σ_{y_{m}})^{2}} + 1] δy_{m} − …  (111b)
⋮  (111c)

or, in matrix form,

[ \begin{matrix}−a_{1} [\frac{σ_{zm}^{2}}{(a_{1}σ_{xm})^{2}} + 1] & −a_{2} & … \\ −a_{1} & −a_{2} [\frac{σ_{zm}^{2}}{(a_{2}σ_{ym})^{2}} + 1] & … \\ ⋮ & ⋮ & ... \end{matrix} ] [ δx_{m} ; δy_{m} ; ⋮ ] = [ Δz_{m} ; Δz_{m} ; ⋮ ]  (112a)

or, symbolically,

A_{m} · δ_{m} = Δz_{m}  (112b)

so that, as in equation 92a, we obtain (δx_{m}, δy_{m}, …) in terms of Δz_{m} from

δ_{m} = A_{m}^{−1} · Δz_{m}  (113)

Finally, as in equations 93 and 95, we compute χ^{2} from

χ^{2} = ∑_{m=0}^{M−1} \frac{(z_{m} − a_{1}x_{m} − a_{2}y_{m} − …)^{2}}{σ_{z_{m}}^{2} + a_{1}^{2}σ_{x_{m}}^{2} + a_{2}^{2}σ_{y_{m}}^{2} + …} = ∑_{m=0}^{M−1} \frac{Δz_{m}^{2}}{σ_{z_{m}}^{2} + a_{1}^{2}σ_{x_{m}}^{2} + a_{2}^{2}σ_{y_{m}}^{2} + …}  (114)

and equation 98a becomes

σ_{meas,m}^{2} = σ_{z_{m}}^{2} + a_{1}^{2}σ_{x_{m}}^{2} + a_{2}^{2}σ_{y_{m}}^{2} + …  (115)

which is used, as above, to calculate W_{m,m} = \frac{1}{σ_{meas,m}}.

14. USING SINGULAR VALUE DECOMPOSITION (SVD)

Occasionally, a normal-equation matrix [α] = XX = XXW is sufficiently ill-posed that inverting it using standard matrix inversion doesn’t work. In its invert function, IDL even provides a keyword called status to check on this (although I find that it is not perfectly reliable; the best indicator of reliability is to check that the matrix product XX##XXI = I, the unitary matrix). In these cases, Singular Value Decomposition (SVD) comes to the rescue.

First, we reiterate the solution to our least-squares problem in terms of the normal equations (equation 10 or its equivalent derived from equations 52)

[α]·a = [β]  (116)

for which the solution is, of course,

a = [α]^{−1}·[β]  (117)

We need to find the inverse matrix [α]^{−1}. Above in this document we did this using simple matrix inversion. SVD provides a bombproof and interestingly informative way to do this.

You don’t need to slog through the phenomologicl deacription of SVD below in §14.1 to use the SVD technique in IDL. Implementing SVD in our least-squares solutions is trivially easy, and we provide the IDL prescription below in §14.3. But be sure to look at §14.2!!

14.1. Phenomenonological description of SVD

For a discussion of the details of SVD, see NR §2.6; for SVD applied to least squares, see NR §15.4; if you are rusty on matrix algebra, look at NR §11.0. For our purposes, SVD is a technique for inverting a square symmetric N × N matrix A (SVD is more general; see NR §2.6)^{11}.

As a preliminary to the disussion of SVD, we need to refer to the matrix equation

A · x = b  (118)

The two vectors x and b both have N dimensions. The matrix A is N × N, so if it is nonsingular then, no matter what the values of b, an x can be found that satisfies this equation. But if A is singular, this is no longer true. For example, if two rows of A are identical, then there is one dimension of b that cannot be reproduced by x, no matter how x is chosen. The term nullity, which we denote by NL, is the number of degenerate rows in A. The term rank, which we denote by R, is the number of independent rows in A. Clearly, R + L = N.

SVD says that the inverse of a matrix can be expressed as a product of three matrices:

A^{−1} = V· [\frac{1}{w} ] ·U^{T}  (119)

where

1. w is the set of N weights (also: the N singular values) of A, and the N × N matrix [\frac{1}{w}] is diagonal, with the diagonal elements being \frac{1}{w_{n}}. Some of these w_{n} might be zero, in which case A is singular. The number of nonzero w_{n} is the rank R, i.e. the number of independent rows in A. The number of zero w_{n} is the nullity NL, i.e. the number of degenerate rows in A.

2. U is a square matrix of N column vectors. R of these are associated with nonzero w_{n}; these are an orthonormal set of basis vectors that, when multiplied by an infinite variety of vectors x, can reproduced a set of R dimensions in vectors b.

11In this subsection §14.1, we replace [alpha] by the A to achieve consistency with NR’s notation in their §2.6.

3. V is a square matrix of N column vectors. NL of these are associated with zero w_{n}; these are an orthonormal set of basis vectors that, when multiplied by an infinite variety of vectors x, can reproduce the set of NL dimensions in b that cannot be attained using U.

The inverse matrix A^{−1} is defined by the equation

A · A^{−1} = I  (120)

and, after a little thought, it is clear that the process of evaluating the inverse of our matrix A can be thought of as solving for the N column vectors x_{n}, separately, using the N equations

A · x_{0} = [1, 0, 0...]  (121a)

A · x_{1} = [0, 1, 0...]  (121b)

…

A · x_{N−1} = [0, 0, ..., 1]  (121c)

and putting the N x vectors as columns in A^{−1}. If A is degenerate, then there will be only R equations that have a solution; these constitute the corresponding basis vectors in U. They are associated with the R nonzero weights; think of these weight as the projections of x into the R dimensions of the right hand side above.

Because of the degeneracy, NL of the above equations have no solution—in other words, no projection onto the NL dimensions with which they are associated. This means that their associated weights are zero.

When calculating A^{−1}, equation 119 uses the inverses of the weights [\frac{1}{w_{n}}]. It’s easy to see why: in order to reproduce a dimension on which x always has a small projection, i.e. one that is represented by a small weight, you need to compensate by multiplying the basis vectors by a big number. And if the projection is zero (i.e., w_{n} = 0), you need to multiply by infinity!

The SVD technique provides the weights w_{n} as well as the sets of vectors U and V. If any w_{n} is small, then you know that the particular row in A that it respresents is degenerate with a different row.

14.2. Important Conclusion for Least Squares!!!

For our task of least-squares fitting, the curvature matrix [α] (the same as A above) is symmetric. Suppose two rows of this matrix are degenerate; then because of the matrix symmetry this degeneracy of rows is also a degeneracy of columns. The degeneracy means that two derived parameters that are associated with those two columns have huge covariance.

This, in turn, means that the formulation of the least squares model is faulty: two parameters represent the same physical quantity.

This is a big advantage for the SVD technique: the model parameter that is associated with any small weight is highly covariant, or degenerate, with another parameter. The SVD technique tells you which parameters have problems. You should then take another look at your physical model and reformulate it to get rid of the covariance. And, of course, SVD not only tells about true degeneracy; it also tells about approaches to degeneracy.

Just exactly how small is “small”? Strictly speaking, this is governed by numerical machine accuracy. Measure this by the ratio of the largest w_{n} to the smallest. This ratio is the condition number. Single precision 32-bit floats have accuracy ∼ 10^{−6}, so if the condition number is larger than 10^{6} you certainly have problems. However, smaller condition numbers are a warning of an approach to degeneracy and should put you on the lookout for problems with your model.

14.3. Obtaining [α]^{−1} in IDL

Obtaining [α]^{−1} in IDL is easy. As above, we represent [α] by XX. Above we used the the usual matrix inversion technique, which in IDL is

XXI = invert(XX)  (122)

For SVD, you simply replace this by a short series of statements. First, decompose XX using SVD:

svdc, XX, w, U, V  (123a)

which derives w, U, and V from XX using SVD. These three derived quantities are the same as discussed above in §14.1. Next, form the N × N matrix WM from the vector of weights w:

WM = fltarr(N, N)  (123b)

WM[(N + 1 ∗ intarr(N)] = w  (123c)

and finally, the inverse is

XXI = V · WM · transpose(U)  (123d)

Then proceed as usual.

Having done the SVD, be sure to check the weights! See §14.2!

15. BRUTE FORCE CHI SQUARED AND THE CURVATURE MATRIX

15.1. Parameter Uncertainties in Brute Force Chi square fitting

There are times when “brute force” least squares is appropriate. For example, if you have a nonlinear problem in which taking derivatives is complicated, and if the number of unknown coefficients is small, then it might be easier to search through the coefficient parameter space, calculate the s^{2} for each combination of parameters, and find the minimum. This provides the best-fit parameter values.

For a least-squares fit, the above process provides the sample variance s^{2} (equation 13). You can evaluate the diagonal elements of the curvature matrix by examining the numerical behavior of s^{2} as the parameter values vary and then use equation 19. For example, in taking a straight average, the curvature matrix [α_{σ}] = M, the number of measurements, so [α_{σ}]^{−1} = M^{−1} and s_{a} = \frac{s^{2}}{M}.

Similarly, for a chi-square fit, the above process provides χ^{2} and \overline{χ}^{2} (equation 46; also equations 52i and 52h). Then use equation 54 (or, if you dare, equation 55). For example in taking a weighted average, the curvature matrix [α] = ∑ \frac{1}{σ_{m}^{2}} so [α]^{−1} = (∑ \frac{1}{σ_{m}^{2}} )^{−1}.

16. NOTATION COMPARISON WITH NUMERICAL RECIPES

I learned least squares from Appendix A of Chauvenet (1863). He didn’t use χ^{2} and didn’t use matrix techniques, but §1 and 2 follows his development quite closely. I wrote the first version of this document before knowing of NR’s treatment, which explains my orientation towards least squares instead of chi-square. I’m fortunate in this approach because it made me realize the pitfalls one can get into with chi-square, as I discuss in §9.

On the other hand, NR describe the least squares approach with some disdain in the discussion of equation (15.1.6) and warn that it is “dangerous” because you aren’t comparing the residuals to the intrinsic inaccuracies of the data. In astronomy, though, more often than not you don’t have an independent assessment of σ_{meas}. But you might know the relative weights, and this is a plus for Chi-square fitting. In any case, heed our warnings about chi-square fitting in §9.

In this writeup I have revised my old notation to agree, partially, with NR’s. This effort wasn’t completely successful because I didn’t read NR very carefully before starting. To make it easier to cross-reference this document with NR, I provide the following table of correspondences (left of the double arrow is ours, to the right is theirs):

X ⟷ A  (124a)

Y ⟷ b  (124b)

X^{T} · X = XX = [α] ⟷ A^{T} · A = [α]  (124c)

XX^{−1} = XXI = [α]^{−1} ⟷ [α]^{−1} = [C] = C  (124d)

I use M for the number of measurements and N for the number of unknown coefficients; NR uses the opposite, so we have

N ⟷ M  (124e)

M ⟷ N  (124f)

REFERENCES

Bevington, P.R. & Robinson, D. 1992, Data Reduction and Error Analysis for the Physical Sciences

(WCB/McGraw-Hill). Chauvenet, W. 1863, A Manual of Spherical and Practical Astronomy, Dover Press. Cowan, G. 1998, Statistical Data Analysis, Clarendeon Press. Press, W.H., Flannery, B.P., Teukolsky, S.A., & Vetterling, W.T. 2001, Numerical Recipes (second

edition), Cambridge University Press. Taylor, J.R. 1997, An Introduction to Error Analysis, University Science Books.

This preprint was prepared with the AAS L^{A}TEX macros v5.0.