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

\begin{matrix}c \\ ◯ \end{matrix}Carl Heiles March 1, 2002

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, NM) and Taylor (1997; T97), and we update the notation to partially conform with NM.

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

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

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. 7 5 A NUMERICAL EXAMPLE AND ITS SOLUTION IN IDL. 9

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

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

7 REJECTING BAD DATAPOINTS: CHAUVENET’S CRITERION 16

8 NONLINEAR LEAST SQUARES 18

9 CHI-SQUARE FITTING, AND WEIGHTED FITTING: DISCUSSION IG-NORING COVARIANCE 21 9.1 The weighted mean: the simplest chi-squared fit . . . . . . . . . . . . . . . . . . . . 22 9.2 The multivariate Chi-square fit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 9.3 Persnickety Diatribe on Choosing σ_{meas} . . . . . . . . . . . . . . . . . . . . . . . . 25

9.3.1 Choosing and correcting σ_{data,m} . . . . . . . . . . . . . . . . . . . . . . . . . 25 9.3.2 Come back to reality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

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

weighted average . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 9.5 The explicit IDL for the weighted and chi square cases . . . . . . . . . . . . . . . . . 27

10 CHI-SQUARE FITTING, AND WEIGHTED FITTING: DISCUSSION IN-CLUDING COVARIANCE 30 10.1 Phenomenological description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 10.2 Calculating the uncertainties of a single parameter . . . . . . . . . . . . . . . . . . . 32 10.3 Calculating the uncertainties of two parameter . . . . . . . . . . . . . . . . . . . . . 33 10.4 Calculating the uncertainties of three parameter . . . . . . . . . . . . . . . . . . . . 33 10.5 Doing these calculations the easy way . . . . . . . . . . . . . . . . . . . . . . . . . . 33

10.6 Important comments about uncertainties . . . . . . . . . . . . . . . . . . . . . . . . 35

11 BRUTE FORCE CHI SQUARED AND THE CURVATURE MATRIX 35

11.1 Parameter Uncertainties in Brute Force Chi square fitting . . . . . . . . . . . . . . . 35

11.2 Example: the single-parameter case . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

12 USING SINGULAR VALUE DECOMPOSITION (SVD) IN IDL 36

12.0.1 Step 0: Introduce a common block that will contain the X matrix . . . . . . 37

12.0.2 Step 1: Populate the X matrix and Y vector almost as before . . . . . . . . 37

12.0.3 Step 2: Call svdfit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

13 NOTATION COMPARISON WITH NUMERICAL RECIPES 39

14 ADDENDUM: FITTING A LINE WHEN BOTH X AND Y HAVE UNCER-TAINTIES 40

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 θ_{m} be the m^{th} measurement of the observed quantity (in this example, θ_{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 θ_{m} has uncertainties in its measurement.

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

A + Bt_{m} = θ_{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} = ∑ θ_{m}  (2a)
A ∑ t_{m} + B ∑ t_{m}^{2} = ∑ t_{m}θ_{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} = θ_{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}θ_{m}  (4a)
A ∑ s_{m}t_{m} + B ∑ t_{m}^{2} = ∑ t_{m}θ_{m} .  (4b)

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

A[s^{2}] + B[st] = [sθ]  (5a)
A[st] + B[t^{2}] = [tθ] .  (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} = θ_{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 θ_{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}[sθ] \\ [tθ] \\ [uθ] \\ [vθ] \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}θ_{0} \\ θ_{1} \\ θ_{2} \\ θ_{3} \\ θ_{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 NM’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 NM’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 NM’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)

(NM 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}(θ_{m} − \overline{θ_{m}})^{2} ,  (13)

where \overline{θ_{m}} are the values for θ_{m} predicted by the derived quantities a. Note the difference: θ_{m} are the observed values, while \overline{θ_{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

δθ_{m} = θ_{m} − \overline{θ_{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{θ_{m}}, which is derived from the data, instead of the truly correct value θ^{∗}. (In formal statistical language, θ^{∗} is the mean of the parent population from which your set of measurements is drawn). If you used the truly correct value θ^{∗}, then the sum would be larger than when using \overline{θ_{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 [θ^{∗} and M] instead of [\overline{θ_{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 [θ^{∗} and M] in equation 14? The reason is obvious: we don’t know θ^{∗}! (If we did, we wouldn’t be doing this analysis!)

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

\overline{Y_{σ}} = [ \overline{θ_{0}} ; \overline{θ_{1}} ; \overline{θ_{2}} ; \overline{θ_{3}} ; \overline{θ_{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 §9.

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 θ 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}:

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.

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 θ are (equation 8c)

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

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{θ_{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)

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.

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.

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

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

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. The true covariance matrix C (as defined in Numerical Recipes for example) is

C = σ_{meas}^{2}XXI  (34)

because C is defined for χ^{2} instead of σ^{2} (see §9).

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.

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 in 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, so much so that we won’t even carry it through for our example.

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: 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!

• 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 (NM §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 NM’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} = θ_{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 θ_{0,m}:

sin(A_{0}t_{m}) + B_{0}t_{m} = θ_{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} = θ_{m} − θ_{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 NM’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 (NM §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 NM’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 sacroscanct 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}(θ_{m} − \overline{θ_{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{(θ_{m} − \overline{θ_{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{(θ_{m} − \overline{θ_{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.3).

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, however, the last of these equations 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)

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

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

9.2. 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 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 θ_{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.

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

Now suppose 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)

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}{M − N} = \frac{δY_{σ}^{T} · δY_{σ}}{σ_{meas}^{2}(M − N)}  (53)

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

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

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}}  (55)

This is in contrast to the result quoted in textbooks (e.g. NM equation 15.4.15, BR equation 7.25), which omits the χ^{2} factor:

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

and, as in the standard textbook solution for the weighted mean case, provides parameter errors that are independent of the datapoint residuals.

Our result, equation 55, 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 55 and 56 are identical.

9.3. Persnickety Diatribe on Choosing σ_{meas}

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

In the previous section, equation 56 taught us that—formally, at least—the variance in the derived fit parameter (or its uncertainty, which is the square root) depends only on the adopted uncertainties σ_{meas,m} and not on the actual variance of the datapoints.

Are you bothered by the fact that the uncertainty of the mean σ_{x0} is independent of the data values? You should be: it is obvious that the data values affect σ_{mean}.

Formally, σ_{x0} 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}.

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 σ_{x_{0}}^{2} is also correct. You can obtain it either as the corresponding diagonal to the covariance matrix (equations 55 and 56, 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.^{3}

9.3.2. Come back to reality

In the case χ^{2} ≈ 1 the dispersion of the observed points σ is equal to the intrinsic dispersion of the datapoints σ_{meas} and the mathematical model embodied in the least-squares 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. Specifically, if the datapoints really have intrinsic measurement dispersion σ_{meas} and you specify that they have 2σ_{meas} when computing (X, Y), then you’ll find χ^{2} = 0.25.

There are, in fact, two ways you can get artificially small values for χ^{2}. One is by overestimating your adopted intrinsic values of σ_{meas}. The other is 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; in other words there is covariance in your measured points. Cowan includes this possibility in his equation 7.4 and also example 7.6.1.

3The Δχ^{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 discussed in §10.

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 56 instead of 55 is disastrous.

Think about it.

9.4. 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.2 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.5. 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^{4} 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)

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)

Optional comment: At this point you might wish to revert to Singular Value Decomposition method using IDL’s SVDFIT function; see §12 for a cookbook. It is rare, but sometimes you will encounter a problem for which the inversion of XXW is unstable. In this case, XXW##XXWI ≠ I (I is the unit matrix). You can tell this with the invert function by calling it with the status keyword; if status ≠ 0, you’ve got problems.

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)

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

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

XXWI = invert(XXW) .  (60)

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

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 above example, which we discussed a bit above in §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 55 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, 55, 56, and 66. If the datapoints were projected downward, i.e. if we take small strips δa_{1} and integrate over δa_{2}, the 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

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

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 parameter

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 parameter

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 easy way

The obvious way to do the 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. 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 NM §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.

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 NM’s discussion surrounding equation (15.6.2). Good luck.

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.

• 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 θ_{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. BRUTE FORCE CHI SQUARED AND THE CURVATURE MATRIX

11.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 σ^{2} for each combination of parameters, and find the minimum.

This procedure certainly provides the least squares solution. However, getting the uncertainties requires a bit more work. You can get the full set of uncertainties if you know the curvature matrix, which is the matrix of second derivatives of the pairwise parameters evaluated at the minimum χ^{2}. Having found the minimum, you can derive these curvatures (equal to half the second derivatives) numerically. Then, to find the uncertainties in the parameters, follow the procedure in §10.5.

11.2. Example: the single-parameter case

It is particularly easy to discuss this if there is only one parameter. Let us consider this case and assume that it is a simple weighted average (which is the least squares solution for a constant). Suppose the best fit χ^{2} value for the unknown coefficient (i.e., the weighted average) is x_{0} and that its variance is σ_{x_{0}}^{2}. Now consider offsets from x_{0}, namely x_{0} + δx. Because χ^{2} is minimized at x_{0}, it is clear that

Δχ^{2} = χ_{x_{0}+δx}^{2} − χ_{x_{0}}^{2} = \frac{1}{2} \frac{d^{2}χ^{2}}{dδx^{2}} δx^{2}  (71)

Generally, the curvature is the relevant element of the curvature matrix [α]. In the particular case of a weighted average, there is just one element so

\frac{1}{2} \frac{d^{2}χ^{2}}{dδx^{2}} = [α_{σ}] = [α_{σ}]_{00} = ∑ (\frac{1}{σ_{meas,m}^{2}} )  (72)

Because we have a single parameter x_{0}, its variance σ_{x0}^{2} is the offset δx_{0}^{2} from x_{0} that gives Δχ^{2} = 1. If we had more than one variable, we’d define ellipsoids in the manner discussed in §10.5.

Are you bothered by the fact that the curvature is independent of the data values? Then go back and read §9.3.

12. USING SINGULAR VALUE DECOMPOSITION (SVD) IN IDL

You don’t need to read this section to use the SVD technique in IDL. Rather, you can just use IDL’s documentation for svdfit. The only problem is that you need to fit to a function that you define in an IDL function. If you are fitting to tabular values, or to a function that is not easily expressed analytically, then read on.

Here we present a technique that works on the X matrix we have previously defined^{5}. Often you populate this matrix with values that are not easily expressible analytically, so the following trick is handy. Fundamentally, what this method does is to replace the analytic function by values in the X matrix and replace the values of the independent variable (t_{m} in our examples) by the measurement number m.

12.0.1. Step 0: Introduce a common block that will contain the X matrix

Let’s call this common block svdinput. The relevant IDL command is common svdinput,x. The only catch is that, in IDL, you can’t use a variable before defining a common block that contains it. So it is probably best to exit IDL and re-enter it, defining this common block as the first statement.

12.0.2. Step 1: Populate the X matrix and Y vector almost as before

Here you populate the X and Y matrices just as before in §5.2 and 9.5, with some small differences. One, the measurement uncertainties are an input to svdfit so we don’t create the weighted versions of X and Y. Two, we must make Y a row vector instead of a column vector. Here we repeat all the IDL commands as in §9.5, but without the verbal commentary, and include this slight modification for Y. We use the numbers in the numerical example of §5.2 to be absolutely unambiguous about what to do. In this example, the number of unknowns N = 3 and the number of measurements M = 4.

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

Y = [142, 168, 211, 251] .  (73b)

and we define ME, the vector of measurement errors, to be

ME = [5., 3., 7., 4.] .  (73c)

5It is a pleasure to thank Jason Wright for telling me about this technique.

where we have invented a four arbitrary values for illustrative purposes. You can, of course, use unity for all of the elements of ME; also, you can call svdfit without specifying any measurement errors, in which case it becomes a standard least squares fit. In the following we assume that you are doing a chi-square fit, i.e. following §9.5.

12.0.3. Step 2: Call svdfit

Now in in §9.5 we replace 59 to 67 (and in §5.2, equation 23 to 30) by the following:

dataindex = indgen(M)  (74a)

where M, we remind you, is the number of measurements; and

a = svdfit(dataindex, Y, N, measure errors = ME, function name =′ svdfcn′,

yfit = ybar, chisq = chi sq, variance = vardc official, sigma = sigdc official,

covar = XXWI, singular = singular)  (74b)

ncov = XXWI/sqrt(vardc official#vardc official)  (74c)

Finally, you have to define the following function.

function svdfcn, dataindex, MIDL
common svdinput
return, x_svd[ *, dataindex]
end

If you are doing least squares fitting and setting all elements of ME equal to unity, then of course chisq is the same as the sample variance s^{2}.

This technique is simpler because it requries fewer IDL statements. However, you have to define the additional function and, if you are using the X SVD matrix, the common block, which must be done at the beginning of the IDL session. This technique is good in those rare instances where inversion of the curvature matrix fails, but it works for all cases, of course. NM §15.4 extols its advantages and mentions its one disadvantage, which is requiring more memory; they recommend that you use it always. We got along without it for 40 years, perhaps because we never did any truly challenging problems.

13. 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 NM’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, NM 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 NM’s. This effort wasn’t completely successful because I didn’t read NM very carefully before starting. To make it easier to cross-reference this document with NM, I provide the following table of correspondences (left of the double arrow is ours, to the right is theirs):

X ⟷ A  (75a)

Y ⟷ b  (75b)

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

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

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

N ⟷ M  (75e)

M ⟷ N  (75f)

14. ADDENDUM: FITTING A LINE WHEN BOTH X AND Y 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 θ, as in equation 1. This essential assumption means that t is known with high precision and all the uncertainty is in θ, and you are minimizing the squares of the residuals in the θ-direction only.

If both variables have uncertainties, then you have to be careful because the essential assumption is violated. Before doing anything, read (lightly) Isobe, Feigelson, Akritas, and Babu (1990, ApJ 364, 104). Then adopt one of the two following procedures:

(1) You may be able to reformulate the problem by defining two new variables, one of which has much less fractional uncertainty than the other. The most straightforward reformulation is to define x = (θ + t) and y = (θ − t). The sum, x, has smaller fractional uncertainty than y so its dispersion might not matter as much. In such a redefinition, you should include a scaling factor on one of the original variables to make the two original dispersions about equal.

Before proceeding with this option, re-read Isobe et al!

(2) More generally, use the prescription in Isobe et al. Lucky for me, I haven’t had to delve into this method—but their formulation looks good!

REFERENCES

Bevington, P.R. & Robinson, D. 1992, Data Reduction and Error Analysis for the Physica 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.