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

\begin{matrix}c \\ ◯ \end{matrix}Carl Heiles April 12, 2007

14. FITTING WHEN MORE THAN ONE MEASURED PARAMETERS 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 ??. 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.

Thanks to Jefferys (1980), the ML formulation of this problem is straightforward. Nevertheless, there is a lot of confusion on such fitting and not-inconsiderable propagation of myth. Before reviewing Jefferys’ 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.

14.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 ??, 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} .  (14.1a)

The set of two normal equations becomes just the single equation

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

or, writing out the sums explicitly,

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

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} ,  (14.2)

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

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

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

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

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!

14.2. Jefferys’ Method: Introduction

Elsewhere we have regarded y as the dependent variable with measured errors and x (and its multiple counterparts) as independent variables that have no measured errors. But sometimes all of the measured parameter’s have uncertainties. In this case the distinction between “dependent” and “independent” variables disappears and we need to treat them symmetrically. Moreover, we can have an arbitrarily large number of variables. This means we need to generalize our notation—and our mathematical technique. Our treatment follows Jefferys (1980) and we will adopt his notation in some measure. To explain the method we will follow his example and present the general case and also show the application to a specific example.

Figure 14.1 compares a conventional fit to y = a_{0} + a_{1}t with a proper fit when both variables have uncertainties. The right-hand panel is the conventional fit in which the measurement uncertainties in t are set equal to zero; the left-hand panel includes the uncertainties in t. Some important differences are apparent:

1. Errors on the left panel are completely specified by errorbars. On the right panel we could use vertical errorbars for y and horizontal ones for t, but these would be sufficiently only if the errors were uncorrelated (as they are for points 0, 2, and 4). However, when they are correlated we must use error ellipses that specify χ^{2} = 1.

2. For the left panel, the best-fit points have the same t values as the datapoints but different y values. For the right panel the best-fit values differ from the data values for both variables.

3. If you look carefully, you’ll see that the fitted slope on the left is a bit smaller than that on the right. This systematic bias arises from ignoring the errors in t.

14.3. The Data Matrix and Vector

First, a word on notation. We distinguish between scalars, vectors, and matrices as follows: for Roman letters, vectors are lower-case bold and matrices are UPPER-CASE BOLD. For Greek Fig. 14.1.— Comparing a conventional fit for to y = a_{0} + a_{1}t (left panel) to a proper one when both measured variables have errors. On the right, the ellipses denote the correlated errors in (t, y); these are the generalization of the errorbars on the left. The right-hand slope is a bit steeper than the left-hand one.

letters, scalars are the letter itself (e.g. α), vectors are single-underlined (e.g. \underline{ϕ}), and matrices are double-underlined (e.g. \underline{\underline{σ}}).

For the general case we have M experiments. In each experiment we measure J parameters. We want to combine these (M × J) measurements to derive N coefficients. We denote the matrix of measured parameters by X, which is (M × J) [using conventional matrix notation, not IDL notation, in which in the vertical dimension (number of rows) is M and the horizontal dimension (number of columns) is J]. We will want to concatenate this matrix to a one-dimensional vector x of length (MJ). Specifically, the M × J datapoint matrix is

X = [ \begin{matrix}x_{0,0} & x_{0,1} & x_{0,2} & … & x_{0,J−1} \\ x_{1,0} & x_{1,1} & x_{1,2} & … & x_{1,J−1} \\ x_{2,0} & x_{2,1} & x_{2,2} & … & x_{2,J−1} \\ . & . & . & … & . \\ . & . & . & … & . \\ . & . & . & … & . \\ x_{M−1,0} & x_{M−1,1} & x_{M−1,2} & … & x_{M−1,J−1} \end{matrix} ]  (14.4a)

We don’t use this big matrix in the solution. Instead, we turn it into a vector in which the first J elements are the data for m = 0, the second J elements are for m = 1, etc. So the vector has dimensions (MJ) × 1, like this:

x = [ \begin{matrix}x_{0,0} \\ x_{0,1} \\ x_{0,2} \\ . \\ . \\ . \\ x_{0,J−1} \\ x_{1,0} \\ x_{1,1} \\ x_{1,2} \\ . \\ . \\ . \\ x_{1,J−1} \\ x_{2,0} \\ x_{2,1} \\ x_{2,2} \\ . \\ . \\ . \\ x_{2,J−1} \\ … \\ … \\ … \\ x_{M−1,0} \\ x_{M−1,1} \\ x_{M−1,2} \\ . \\ . \\ . \\ x_{M−1,J−1} \end{matrix} ]  (14.4b)

One important reason for writing the whole set of data as a vector instead of a matrix is to make it possible to write the covariance matrix for all measured data in the conventional form, as we now discuss.

14.4. The Data Covariance Matrix and Defining Chi-Square

The whole object of fitting is to minimize the chi-square. When all measured parameters have uncertainties, their uncertainties can be correlated. We have to generalize the definition of chi-square accordingly.

First, suppose that the observational errors in the datapoints are uncorrelated. Then the intrinsic variance of each datapoint is described by a single number. In our example, for uncorrelated observational errors we’d have the variances in the y values be (σ_{y0}^{2}, σ_{y1}^{2}, …), and similarly for the t values; this would give

χ^{2} = ∑_{m} \frac{δy_{m}^{2}}{σ_{ym}^{2}} + \frac{δt_{m}^{2}}{σ_{tm}^{2}}  (14.5)

However, it is common that errors are correlated. For example, if we were fitting y to a polynomial in t, then the errors in the various powers of t would certainly be correlated. More generally, then, the covariances among the different measured values are nonzero. These covariances are the off-diagonal terms in the covariance matrix. Thus, if we denote the covariance matrix for the measured x values by the (MJ × (MJ) matrix \underline{\underline{σ}}, then the most general case has \underline{\underline{σ}} with no nonzero elements.

Less general, but much more common, is the situation shown in Figure 14.1. Here, the covariances among the J measured parameters nonzero are for a particular experiment m, but the covariances from one experiment m to another are zero; in other words, each experiment is completely independent of the others. In this less general but very common case, the covariance matrix looks like

\underline{\underline{σ}} = [ \begin{matrix}\underline{\underline{σ_{0}}} & 0 & 0 & … \\ 0 & \underline{\underline{σ_{1}}} & 0 & … \\ 0 & 0 & \underline{\underline{σ_{2}}} & … \\ . \\ . \\ . \end{matrix} ]  (14.6)

Here, each element (including the 0 elements) is itself a J × J matrix. For our specific example of §14.7, J = 2 so \underline{\underline{σ_{0}}} is a covariance matrix of the form

\underline{\underline{σ_{0}}} = [ \begin{matrix}σ_{yy} & σ_{yt} \\ σ_{yt} & σ_{tt} \end{matrix} ]  (14.7)

Generally, the chi-square is given by (e.g. Cowan equation 2.39)

χ^{2} = δx^{T} · \underline{\underline{σ^{−1}}} · δx  (14.8)

14.5. Formulation of the Problem and its Solution with Lagrange Multipliers

We will be referring to various versions of the data parameters x and derived parameters a: measured, best-fit, and (for the iterative solution) guessed. The subscript d denotes the set of measured datapoints, of which there are (JM). The subscript ∗ denotes the set of best-fit quantities; these parameters include not only the datapoints x, but also the derived parameters a. We will be doing an iterative fit using guessed values of both the data and derived parameters, represented by the subscript g.

We begin by writing exact equations for each measurement. The fitted vales, subscripted with stars, satisfy the exact equations of condition

f(x_{∗}, a_{∗}) = 0  (14.9a)

This is an M-long vector of functions f(x, a) = 0 (one row for each measurement). This set of M equations doesn’t do us much good because we don’t know the best-fit (starred) values. Consequently, for the datapoints we define the difference between the best-fit and measured data values

δx = x_{d} − x_{∗}  (14.9b)

This is the negative of Jefferys’ definition of the corresponding quantity v̂ in his section II. With this, the equation 14.9a becomes

f(x_{d} − δx, a_{∗}) = 0 .  (14.9c)

Our goal is to solve these M equations for the (MJ) differences δx and the N parameters a_{∗} and, simultaneously, minimize χ^{2}.

This is a classic minimization problem: we minimize χ^{2} with respect to the (MJ +N) values of δx and a, subject to the M constraints of equation 14.9c. Such problems are solved using Lagrange multipliers. Here, the M Lagrange multipliers form the vector \underline{λ}. We define the Lagrangian L as

L = [\frac{1}{2}δx^{T} · \underline{\underline{σ^{−1}}} · δx] + [f^{T}(x_{d} − δx, a) · \underline{λ}] ;  (14.10)

the \frac{1}{2} arises because, for a Gaussian pdf for the errors, the residuals are distributed as e^{− \frac{χ^{2}}{2}} (e.g. Cowan equation 2.28). We differentiate L with respect to each of the unknowns δx and a. The solution provides a, the vector of N derived coefficients (so a is defined as elsewhere in this tutorial) together with the MJ fitted datapoints. To proceed, we need the derivatives of f in equation 14.9a with respect to the vectors x and a, as we now discuss.

14.6. The Derivative Matrices

For the analysis we will need the derivatives of f with respect to the vectors x and a. The derivatives will always be evaluated at the guessed values x_{g} and a_{g}. The derivative with respect to x is a M × MJ matrix and looks like this [We take the case (M = 4, J = 3) for transparency; subscripts of x are in the order (m, j)]:

\frac{∂f(x, a)}{∂x} |_{x_{g},a_{g}} = [ \begin{matrix}\frac{∂f}{∂x_{0,0}} & \frac{∂f}{∂x_{0,1}} & \frac{∂f}{∂x_{0,2}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{∂f}{∂x_{1,0}} & \frac{∂f}{∂x_{1,1}} & \frac{∂f}{∂x_{1,2}} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{∂f}{∂x_{2,0}} & \frac{∂f}{∂x_{2,1}} & \frac{∂f}{∂x_{2,2}} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{∂f}{∂x_{M−1,0}} & \frac{∂f}{∂x_{M−1,1}} & \frac{∂f}{∂x_{M−1,2}} \end{matrix} ] (14.11)

where all the derivatives are evaluated at (x_{g}, a_{g}). Much easier is the derivative with respect to a, which is a M × N matrix and looks like

\frac{∂f(x, a)}{∂a}|_{x_{g},a_{g}} = [ \begin{matrix}\frac{∂f}{∂a_{0}} |_{x_{g,0},a_{g}} & \frac{∂f}{∂a_{1}} |_{x_{g,0},a_{g}} & \frac{∂f}{∂a_{2}} |_{x_{g,0},a_{g}} & ... & \frac{∂f}{∂a_{N−1}} |_{x_{g,0},a_{g}} \\ \frac{∂f}{∂a_{0}} |_{x_{g,1},a_{g}} & \frac{∂f}{∂a_{1}} |_{x_{g,1},a_{g}} & \frac{∂f}{∂a_{2}} |_{x_{g,1},a_{g}} & ... & \frac{∂f}{∂a_{N−1}} |_{x_{g,1},a_{g}} \\ \frac{∂f}{∂a_{0}} |_{x_{g,2},a_{g}} & \frac{∂f}{∂a_{1}} |_{x_{g,2},a_{g}} & \frac{∂f}{∂a_{2}} |_{x_{g,2},a_{g}} & ... & \frac{∂f}{∂a_{N−1}} |_{x_{g,2},a_{g}} \\ & & . \\ & & . \\ & & . \end{matrix} ]  (14.12)

where again, all the derivatives are evaluated at (x_{g}, a_{g}).

14.7. The Specific Example

We illustrate the above with this specific example, for which we fit a first-order polynomial to (t, y) of the form [y = a_{0} + a_{1}t]; the coefficients are (a_{0}, a_{1}) and [f(x, a) = y − a_{0} − a_{1}t]. For this example, the matrix of measured parameters is

X = [ \begin{matrix}y_{d0} & t_{d0} \\ y_{d1} & t_{d1} \\ y_{d2} & t_{d2} \\ . \\ . \\ . \end{matrix} ]  (14.13)

and the concatenated vector version is

x = [ \begin{matrix}y_{d0} \\ t_{d0} \\ y_{d1} \\ t_{d1} \\ y_{d2} \\ t_{d2} \\ . \\ . \\ . \end{matrix} ]  (14.14)

The exact vector-of-functions equation f(x_{∗}, a_{∗}) = 0 uses the starred (best-fit) values, and is

f(x_{∗}, a_{∗}) = [ \begin{matrix}y_{∗0} − a_{∗0} − a_{∗1}t_{∗0} \\ y_{∗1} − a_{∗0} − a_{∗1}t_{∗1} \\ . \\ . \\ . \end{matrix} ] = [ \begin{matrix}0 \\ 0 \\ . \\ . \\ . \end{matrix} ]  (14.15)

The derivative matrix of f with respect to the vector x is

\frac{∂f(x, a)}{∂x}|_{x_{g},a_{g}} = [ \begin{matrix}1 & −a_{g1} & 0 & 0 & 0 & 0 & … \\ 0 & 0 & 1 & −a_{g1} & 0 & 0 & … \\ 0 & 0 & 0 & 0 & 1 & −a_{g1} & … \\ & & & . \\ & & & . \\ & & & . \end{matrix} ] ,  (14.16)

always evaluated at the guessed values of the parameters (subscript g). The derivative matrix of f with respect to the vector a is

\frac{∂f(x, a)}{∂a}|_{x_{g},a_{g}} = [ \begin{matrix}−1 & −t_{g0} \\ −1 & −t_{g1} \\ −1 & −t_{g2} \\ . \\ . \\ . \end{matrix} ] ,  (14.17)

again always evaluated at the guessed values of the parameters (subscript g).

14.8. The Solution to the Lagrangian: Two Matrix Equations

Jefferys does all this this^{1} for us and, after some algebraic manipulation, provides the two following matrix equations (from the set of four in his equation 10):

\underline{\underline{σ^{−1}}} · δx + \frac{∂f^{T}(x, a)}{∂x}|_{x_{d},a_{∗}} · \underline{λ} = 0  (14.18a)

This single matrix equation embodies MJ individual equations. Here we write \frac{∂f^{T}}{∂x} instead of \frac{∂f^{T}}{∂δx} for clarity, and use the fact that they are the negative of each other.

\frac{∂f^{T}(x, a)}{∂a}|_{x_{d},a_{∗}} · \underline{λ} = 0  (14.18b)

This matrix equation embodies N individual equations. These two matrix equations embody MJ +N individual equations, which is equal to the number of unknowns, so we can solve for them! In both of these equations, the dimensions of the factors are (Note the transposes!):

\underline{\underline{σ}} : (MJ) × (MJ)  (14.19a)
δx : MJ × 1  (14.19b)
f^{T} : 1 × M  (14.19c)
\frac{∂f^{T}(x, a)}{∂x} : (MJ) × M  (14.19d)

1And more: he includes the possibility for additional constraints, the second line in his equation (7). This introduces additional complications and we ignore this possibility in the interest of pedagogical simplification. For example, his set of equations (10) has four equations, not two as we write here.

\underline{λ} : M × 1  (14.19e)
\frac{∂f^{T}(x, a)}{∂a} : N × M  (14.19f)

Note that in 14.19d above the dimension is (MJ) × M: M elements in f, differentiated by (MJ) different variables x; similarly for 14.19f above.

14.9. Solving Equations 14.18a and 14.18b Iteratively

Generally, these equations are nonlinear and we solve them iteratively using guessed solutions. We denote the guessed values with subscript g, so we have

a_{g} = guessed values for a_{∗}  (14.20a)

x_{g} = guessed values for x_{∗}  (14.20b)

ITERATION STEP 1: We define the difference between the measured data quantities and their currently-guessed counterparts

Δx_{g} ≡ x_{d} − x_{g}  (14.20c)

Above, our Δx_{g} is the negative of Jefferys’ v̂. The nonlinear fit solves for corrections to these guesses, which we denote by Δx_{new} (the negative of Jefferys’ v̂_{new}) and Δa_{new} (which is identical to Jefferys’ ˆ \underline{δ}).

ITERATION STEP 2: We define the following:

1. The M × M weight matrix W is from Jefferys’ equation (15):

W ≡ [ \frac{∂f(x, a)}{∂x}|_{x_{g},a_{g}} · \underline{\underline{σ}} · \frac{∂f^{T}(x, a)}{∂x}|_{x_{g},a_{g}}]^{−1}  (14.21a)

2. The [0, 0] element of his equation (17), which is equivalent to our \underline{\underline{α}} (the curvature matrix) elsewhere in this document:

\underline{\underline{α}} ≡ \frac{∂f^{T}(x, a)}{∂a}|_{x_{g},a_{g}} · W · \frac{∂f(x, a)}{∂a}|_{x_{g},a_{g}}  (14.21b)

\underline{\underline{α}} is N × N.

3. The modified equations of condition from his equation (18) (we have a + instead of his − because Δx_{g} = −v̂)

\underline{ϕ_{g}} ≡ f(x_{g}, a_{g}) + ( \frac{∂f(x, a)}{∂x}|_{x_{g},a_{g}}) · Δx_{g}  (14.21c)

\underline{ϕ_{g}} is M × 1.

ITERATION STEP 3: The solutions follow directly. The matrix equation for Δa_{new} is Jefferys equation (17)

\underline{\underline{α}} · Δa_{new} = − \frac{∂f^{T}(x, a)}{∂a}|_{x_{g},a_{g}} · W · \underline{ϕ_{g}}  (14.22)

which is solved for Δa_{new} in the conventional way by premultiplying both sides by \underline{\underline{α}}^{−1}. The matrix equation for the new, corrected Δx_{g} is Jefferys equation (19)

Δx_{new} = −\underline{\underline{σ}} · \frac{∂f^{T}(x, a)}{∂x}|_{x_{g},a_{g}} · W · (\underline{ϕ_{g}} + \frac{∂f}{∂a}|_{x_{g},a_{g}} · Δa_{new})  (14.23)

ITERATION STEP 4: The previous two equations provide corrections to the values x_{d} and a_{g}. One applies them (Jefferys equation 20):

a_{g,new} = a_{g} + Δa_{new}  (14.24a)
x_{g,new} = x_{d} + Δx_{new}  (14.24b)

Note that Δx_{new} is applied to x_{d}, not to x_{g}.

ITERATION STEP 5: Return to Iteration Step 1, using these new values a_{g,new} and x_{g,new} in place of a_{g} and x_{g}. Iterate until convergence. Convergence occurs when all derived corrections Δa_{new} and Δx_{new} become negligible. [How do you define “negligible”…]

14.10. Taking all those derivatives!

Why not use numerical derivatives? It’s easier with complicated functions!

14.11. The Initial Guess

Getting started is the hard part unless the problem is relatively easy. One way is to use conventional least squares to derive the initial guesses a_{g}. To do this, you designate one particular variable as the dependent variable and all the others as the independent ones (for which you assume no errors). Do the conventional least squares solution for the parameters a and use these as the initial guesses. A good choice for the dependent variable is the one with the largest errors, if there is one; otherwise the choice is more-or-less arbitrary. For the initial guesses of the data parameters x_{g}, just use the measured data values x_{d}. If this scheme doesn’t work, you’re on your own!

14.12. The Covariance Matrix (and errors) of the Derived Parameters

The matrix \underline{\underline{α}}, which is defined above in equation 14.21b, is the conventionally-defined curvature matrix and its inverse \underline{\underline{α^{−1}}} is the conventionally-defined covariance matrix. Because we have formulated this problem as a chi-squared one, the elements of this matrix give the covariance directly. Thus, the errors in the derived parameters can be taken as the square-root of the diagonal elements of \underline{\underline{α^{−1}}}.

One usually wants to calculate the chi-squared of the solution. This involves not only the best-fit parameters a_{∗} but also the best fit datapoints x_{∗}. To do this, use equation 14.8 using Δx_{g} in place of δx. This is the same as Jefferys’ equation (43).

The expectation value of χ^{2} is the number of degrees of freedom. In one sense this is just like the usual least-squares solution: it’s equal to the number of datapoints minus the number of derived parameters. Here the number of datapoints is not the number of experiments M; rather, it’s JM. So the number of degrees of freedom is JM − N.