\documentclass[psfig,preprint]{aastex}
\begin{document}

\title{LEAST- AND CHI-SQUARES FOR THE BUDDING AFICIONADO: \\ 
ART AND PRACTICE } 

\author{\copyright Carl Heiles \today }

	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 \S\ref{numexample},  and the {\it
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 (\S \ref{sectionone} to
\ref{chauvenetsection}) 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 $\chi^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 $\chi^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 $\chi^2$ one. In later sections we generalize to
the $\chi^2$ and this other weighted-observations case.

\tableofcontents

\section{LEAST-SQUARE FITTING FOR TWO PARAMETERS, AS WITH A STRAIGHT
LINE.} \label{sectionone}

\subsection{The closed-form expressions for a straight-line fit}

	First consider the least squares fit to a straight line.  Let
$\theta_m$ be the $m^{th}$ measurement of the observed quantity (in this
example, $\theta_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 least-squares technique,
quantities such as $t_m$ are regarded to be known with high accuracy
while the quantity $\theta_m$ has uncertainties in its measurement. 

	We expect the zenith distance $\theta_m$ to change linearly with
time as follows:

\begin{equation} \label{one}
A + B t_m = \theta_m \; .
\end{equation}

\noindent 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 {\it
normal equations}

\begin{mathletters} \label{normalone}
\begin{equation}
AN + B \ \sum t_m = \sum \theta_m
\end{equation}
\begin{equation}
A \ \sum t_m + B \ \sum t_m^2 = \sum t_m \theta_m
\end{equation}
\end{mathletters}

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

\subsection{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
\ref{one}, $s_m = 1$. Then we rewrite equation \ref{one} to include this
extra dependence

\begin{equation} \label{two}
A s_m + B t_m = \theta_m \; .
\end{equation}

\noindent 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 \ref{two}, one for each
measurement.  They are known as the {\it 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 \ref{two} we form the {\it
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 \ref{normalone}; the result is

\begin{mathletters} \label{ones}
\begin{equation}
A\ \sum s_m^2 + B \ \sum s_m t_m = \sum s_m\theta_m
\end{equation}
\begin{equation}
A \ \sum s_m t_m + B \ \sum t_m^2  = \sum t_m \theta_m \ .
\end{equation}
\end{mathletters}

\noindent We can rewrite these equations using the notation $[st] =
\sum s_m t_m$, etc.:

\begin{mathletters} \label{normaltwo}
\begin{equation}
A [ s^2 ] + B [ s t ] = [ s \theta ]
\end{equation}
\begin{equation}
A [ s t ] + B [ t^2 ] = [ t \theta ] \ .
\end{equation}
\end{mathletters}

\noindent This is, of course, precisely analogous to equation
\ref{normalone}. And now it's clear how to generalize to more
parameters!

\section{LEAST-SQUARE FITTING FOR MANY PARAMETERS, AS WITH A CUBIC}
\label{sectiontwo}

	With this notation it's easy to generalize to more ($N$)
unknowns: the method is obvious because in each equation of condition
(like equation \ref{two}) we simply add equivalent additional terms such
as $C u_m$, $D v_m$, etc; and in the normal equations (equation
\ref{normaltwo}) 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 \ref{ones}
is the $M$ equations

\begin{equation} \label{eqcond}
A s_m + B t_m + C u_m + D v_m = \theta_m \; ,
\end{equation}

\noindent with $m = 0 \rightarrow (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
$\theta_m$.  We then form the four normal equations; the generalization
of equation \ref{normaltwo} written in matrix format is:

\begin{eqnarray} \label{smeqn} 
\left[ 
\begin{array}{cccc} 
{[ ss ]} & {[ st ]} & {[ su ]} & {[ sv ]} \\ 
{[ ts ]} & {[ tt ]} & {[ tu ]} & {[ tv ]} \\ 
{[ us ]} & {[ ut ]} & {[ uu ]} & {[ uv ]} \\ 
{[ vs ]} & {[ vt ]} & {[ vu ]} & {[ vv ]} \\
\end{array} 
\; \right] 
\left[ 
\begin{array}{c} 
A \\ 
B \\ 
C \\ 
D \\
\end{array} 
\; \right] 
\; = 
\left[ 
\begin{array}{c} 
{[ s \theta ]} \\ 
{[ t \theta ]} \\ 
{[ u \theta ]} \\ 
{[ v \theta ]} \\ 
\end{array} 
\; \right]
\end{eqnarray} 

\noindent The $N \times N$ matrix on the left is symmetric. With $N$
equations and $N$ unknowns, you can actually {\it solve} for $A, B, C,
D$! 

\section{FAR, FAR BEST AND EASIEST: MATRIX ALGEBRA.} \label{matrixmethod}

	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 {\bf matrix} by {\bf
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

\begin{mathletters} \label{matrixdefinition}
\begin{eqnarray} 
{\bf X_\sigma} = \left[
\begin{array}{cccc} 
{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{array} 
\; \right] 
\end{eqnarray}
%\end{mathletters}

%\begin{mathletters}
\begin{eqnarray}
{\bf a_\sigma} = \left[
\begin{array}{c}
A \\
B \\
C \\
D \\
\end{array} \; \right] 
\end{eqnarray}
%\end{mathletters}

%\begin{mathletters}
\begin{eqnarray}
{\bf Y_\sigma} = \left[
\begin{array}{c} 
{\theta_0} \\ 
{\theta_1} \\ 
{\theta_2} \\ 
{\theta_3} \\ 
{\theta_4} \\ 
\end{array} 
\; \right]
\end{eqnarray} 
\end{mathletters} 

\noindent so, in matrix form, the equations of condition (equation
\ref{eqcond}) reduce to the single matrix equation

\begin{equation}
\label{equationofcondition}
{\bf X_\sigma \cdot a} = {\bf Y_\sigma} \; .
\end{equation}

\noindent The notation for these equations corresponds to NM's. We write
them with subscripts $\sigma$ to emphasize that they are calculated
without dividing by $\sigma_{meas}$, i.e.\ that we are doing least
squares instead of chi square fitting. For chi square fitting, see \S
\ref{chisqsection} and \ref{chicoeffs}.

	Our matrix $\mathbf{X_\sigma}$ corresponds to NM's ``design
matrix'' $\bf A$ of Figure 15.4.1, except that our elements are not
divided by $\sigma_{meas,m}$, and the matrix equation of condition
(equation \ref{equationofcondition}) 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
$\sigma_{meas,m}$, the intrinsic measurement uncertainties (\S
\ref{chisqsection}).

	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
\ref{smeqn}) reduce to the single equation 

\begin{equation} \label{matrixnormal}
\mathbf{[\alpha_\sigma] \cdot a} = \mathbf{[\beta_\sigma]} \; ,
\end{equation}

\noindent (NM equation 15.4.10), where

\begin{mathletters} 
\label{ssdef}
\begin{equation} 
[\alpha_\sigma] = {\bf X_\sigma^T \cdot X_\sigma} 
\end{equation}
\begin{equation}
[\beta_\sigma] = {\bf X_\sigma^T \cdot Y_\sigma}
\end{equation}
\end{mathletters}

\noindent The matrix $[\alpha_\sigma]$ is known as the {\it curvature
matrix} because each element is twice the curvature of $\sigma^2$ (or
$\chi^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 $[\alpha_\sigma]$ (that is, by
$[\alpha_\sigma]^{-1})$, which gives

\begin{equation}
\label{adef}
{\bf a} = [\alpha_\sigma]^{-1} \mathbf{ \cdot [\beta_\sigma]} \; .
\end{equation}

\noindent All of these matrix operations are trivially easy in IDL (\S
\ref{numexample}). 

\section{UNCERTAINTIES IN THE DERIVED COEFFICIENTS.} \label{sigmas}

	How about the uncertainties in the derived quantities contained
in the matrix ${\bf a}$?

       The first thing to do is derive the {\it 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 = [ \sum_0^{M-1} (x_m -
\overline{ x_m})^2/(M-1)]$.  The generalization is, simply, to replace   
the $M-1$ in the denominator by $\nu = M-N$.  In the straight-average
case, $N=1$ so this fits. Here $\nu$ is known as the number of {\it
degrees of freedom} and $N$, the number of unknown coefficients, is
known as the number of {\it constraints}. So we have

\begin{equation} \label{samplevarianceone}
s^2 = {1 \over M - N} \sum_{m=0}^{M-1} (\theta_m -
\overline{\theta_m})^2 \; ,
\end{equation}  

\noindent where $\overline{\theta_m}$ are the values for $\theta_m$ {\it
predicted by the derived quantities ${\bf a}$}. Note the difference:
$\theta_m$ are the {\it observed} values, while $\overline{\theta_m}$ are the
values {\it predicted by the least squares fit}.  The predicted values
are those that are computed from the derived coefficients $A, B, C$\dots
The $M$ quantities

\begin{equation} \label{sequation}
\delta \theta_m = \theta_m - \overline{\theta_m}
\end{equation}

\noindent are called the {\it residuals} or {\it 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 ${M \over 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 {\theta_m}$, which is derived from
the data, instead of the {\it truly} correct value $\theta^*$. (In
formal statistical language, $\theta^*$ is the mean of the parent
population from which your set of measurements is drawn). If you used
the truly correct value $\theta^*$, then the sum would be larger than
when using $\overline{ \theta_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 $\sigma^2$, which you'd get by
using [$\theta^*$ and $M$] instead of [$\overline{\theta_m}$ and
$(M-N)$] in equation \ref{sequation}; 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 [$\theta^*$ and $M$]
in equation \ref{sequation}? The reason is obvious: we don't know
$\theta^*$! (If we did, we wouldn't be doing this analysis!)

	It's easy to calculate the $\overline{\theta_m}$ with matrices.
First define the matrix ${\bf \overline{ Y_\sigma}}$ that contains these
values:

\begin{eqnarray}
{\bf \overline{ Y_\sigma}} = \left[
\begin{array}{c} 
\overline{\theta_0} \\ 
\overline{\theta_1} \\ 
\overline{\theta_2} \\ 
\overline{\theta_3} \\ 
\overline{\theta_4} \\ 
\end{array} 
\; \right]
\end{eqnarray} 

{\noindent} Calculating ${\bf \overline{ Y_\sigma}}$ is simple:

\begin{equation} \label{ybardefinition}
{\bf \overline{ Y_\sigma}} = {\bf  X_\sigma \cdot a} \; . 
\end{equation}

\noindent Note that ${\bf X_\sigma}$ is already defined (equation
\ref{matrixdefinition}) and ${\bf a}$ was solved for in equation
\ref{adef}. It's convenient to define the residual matrix

\begin{equation}
{\bf \delta Y_\sigma } = {\bf Y_\sigma - \overline{ Y_\sigma}}
\end{equation}

\noindent so we can write

\begin{equation} \label{samplevariancetwo}
s^2 = {1 \over M - N} {\bf \delta Y_\sigma^T \cdot \delta Y_\sigma} \; . 
\end{equation}

	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 \S \ref{sectionone}. We won't go through
the derivation here; you can find it in Taylor \S 8.4 and equation 8.16,
8.17. The result is

\begin{equation} \label{coeffvarianceone}
{\bf s_a}^2 = s^2 diag \{[\alpha_\sigma]^{-1} \}
\end{equation}

\noindent Or, to put it simply in words: to get the variance of 
coefficient $n$ in the matrix ${\bf a}$, multiply $s^2$ by the $n^{th}$
diagonal element of $[\alpha_\sigma]^{-1}$. 

	Although the above equation for ${\bf 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 \S \ref{chisqsection}. 

\section{A NUMERICAL EXAMPLE AND ITS SOLUTION IN IDL.} \label{numexample}

	If the following sounds like Greek to you, take a look at \S
\ref{matrixmethod} and \ref{sigmas}.

\subsection{Generation of the numerical example}

	Suppose that we make four measurements of the angle $\theta$ and
we want to fit to a parabolic function in time $t$.  In the notation of
equation~\ref{eqcond}, $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
\rightarrow 3$.  Suppose that the four values of time are 5, 7, 9, 11. 

\begin{figure}[h!]
\begin{center}
\leavevmode
%\includegraphics[height=7.5in, width=6.0in]{twolvl3b.ps}
\includegraphics[scale=.55, angle=90]{lsfitfig.ps}
%\includegraphics{twolvla.ps}
\end{center}
\caption{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.\label{lsfitfig}}
\end{figure}

	First we create the matrix {\bf X} in IDL

\begin{equation}
{\bf X = fltarr(N,M) = fltarr(3,4)}
\end{equation}
	
\noindent 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 \ref{matrixdefinition}a we have $s_m = 1$,
$t_m = time_m$, $u_m = time_m^2$:

\begin{mathletters}
\begin{eqnarray}
{\bf X} = \left[
\begin{array}{ccc}
1  &  5   &  25  \\
1  &  7   &  49  \\
1  &  9   &  81  \\
1  &  11  &  121 \\
\end{array} \; \right] 
\end{eqnarray}
\end{mathletters}

\noindent Suppose that the four measured values of $\theta$ are
(equation~\ref{matrixdefinition}c)

\begin{mathletters}
\begin{eqnarray}
{\bf Y} = \left[
\begin{array}{c}
142 \\
168 \\
211 \\
251 \\
\end{array} \; \right] \; .
\end{eqnarray}
\end{mathletters}

\noindent Figure \ref{lsfitfig} 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 $\bf Y$ above, you have to define
${\bf 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 ${\bf Y} = [142,168,211,251]$; you can make it into the
necessary column vector by taking its transpose, i.e.\ ${\bf Y =
transpose(Y)}$. 

	
\subsection{Solution of the Numerical Example in IDL} \label{idlsolution}

	In IDL we calculate the normal equation matrices and denote the
$[\alpha]$ in equation \ref{ssdef}a by ${\bf XX}$:

\begin{mathletters} \label{eqn24}
\begin{equation}
{\bf XX} = {\bf transpose(X) \# \# X} \; ,
\end{equation}

\noindent and we denote the $[\beta]$ in equation \ref{ssdef}b by {\bf
XY}:

\begin{equation}
{\bf XY} = {\bf transpose(X) \# \# Y} \; .
\end{equation}
\end{mathletters}

\noindent In IDL we take the inverse of $[\alpha]$ (same as $\bf XX$) by

\begin{equation}
{\bf XXI} = {\bf invert(XX)} \; .
\end{equation}

	The least-square fitted quantities are in the matrix {\bf a}
(equation \ref{adef}), which we obtain in IDL with

\begin{equation}
{\bf a} = {\bf XXI \ \#\# \ XY} \; .
\end{equation}

	In IDL we denote the matrix of predicted values $\overline{
\theta_m}$ by ${\bf YBAR}$, which is 

\begin{equation}
\label{bteqn}
{\bf YBAR} = {\bf X \ \#\# \ a} \; .
\end{equation}

\noindent and we can also define the residuals in $\bf Y$ as

\begin{equation}
{\bf DELY} = {\bf Y - YBAR}
\end{equation}

\noindent In IDL we denote $s^2$ in equations \ref{samplevarianceone}
and \ref{samplevariancetwo} by $s\_sq$ and write 

\begin{mathletters}
\label{sigsq}
\begin{equation}
s\_sq = {\bf transpose(DELY) \#\# DELY}/(M-N) \; . 
\end{equation}

\noindent or

\begin{equation}
s\_sq = total({ \bf DELY} \wedge 2)/(M-N) \; . 
\end{equation}
\end{mathletters}

\noindent It is {\it always} a good idea to plot all three quantities
(the measured values {\bf Y}, the fitted values {\bf YBAR}, and the
residuals ${\bf 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 \times N$
matrix has $N$ diagonal elements; a convenient way to get them is

\begin{equation}
diag \ elements \ of \ {\bf XXI}= {\bf XXI[(N+1) * indgen(N)]}
\end{equation}

\noindent In IDL, we define the variances of
the $N$ derived coefficients by {\bf vardc} (think of ``{\bf var}iances of
{\bf d}erived {\bf c}oefficients''). You can get this as in equation
\ref{coeffvarianceone} from\footnote{If you used equation \ref{sigsq}a
instead of \ref{sigsq}b, then IDL considers $s\_ sq$ an array and you
need to use a $\#$ instead of a $*$ in this equation.}

\begin{equation} \label{lscoefferrorone}
{\bf vardc} = s\_sq * {\bf XXI[(N+1)*indgen(N)]} \; .
\end{equation}

\subsection {Discussion of the numerical example}

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

\begin{mathletters}
\label{avalues}
\begin{eqnarray}
{\bf a} = \left[
\begin{array}{c}
96.6250 \\
4.5000 \\
0.8750 \\
\end{array} \; \right] 
\end{eqnarray}

\noindent and the errors in the derived coefficients [the square root of
the $\sigma^2$'s of the derived coefficients, i.e.  $[\sigma_n^2]^{1/2}$
or, in IDL, $sqrt({\bf vardc})$ in equations \ref{lscoefferrorone}] are:

\begin{eqnarray}
{\bf \sigma_A} = \left[
\begin{array}{c}
34.012 \\
9.000 \\
0.5590 \\
\end{array} \; \right] \; .
\end{eqnarray}
\end{mathletters}

\noindent These results look {\it 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! 

\section{THE COVARIANCE MATRIX AND ITS NORMALIZED COUNTERPART}

\label{ncov}

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

\subsection{ Definition of the normalized covariance matrix}

	The variances in the derived coefficients are obtained from the
diagonal elements of ${\bf XXI}$. The off-diagonal elements represent
the {\it covariances} between the derived coefficients. Covariance means,
specifically, the degree to which the {\it uncertainty} in {\it one}
derived coefficient affects the uncertainty in {\it 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 {\it per se}, we first discuss its normalized
counterpart. 

	The normalized covariance
matrix\footnote{It is a pleasure to thank Doug Finkbiener for
introducing me to this concept.} {\bf ncov} is derived from ${\bf XXI}$
by dividing each element by the square root of the product of the
corresponding diagonal elements. Let ${\bf ncov}$ be the normalized
covariance matrix; then

\begin{equation}
ncov_{ik} = {XXI_{ik} \over \sqrt {XXI_{ii} \ XXI_{kk}}}
\end{equation}

\noindent This is the same normalization that one does with the Pearson
linear correlation coefficient.  In fact, the elements of the normalized
covariance matrix {\it are} the correlation coefficients.  In IDL, you
do the following:
 
\begin{mathletters} \label{idlcovariance}
\begin{equation}
{\bf dc = XXI[(N+1)*indgen(N)]}
\end{equation}
\begin{equation}
{\bf ncov = XXI/{\it sqrt}(dc \# dc)}
\end{equation}
\end{mathletters}

\noindent In the above, ${\bf dc \# dc}$ is an $N \times N$ matrix
consisting of products of the diagonals of ${\bf XXI}$, so dividing
${\bf XXI}$ by $sqrt({\bf dc \# dc})$ generates the normalized version. 

	Because {\bf ncov} is a {\it normalized} covariance matrix, you
might think that it's non-normalized parent is {\bf XXI}---and you'd be
{\it almost} right.  The true covariance matrix ${\bf C}$ (as defined in
{\it Numerical Recipes} for example) is 

\begin{equation}
{\bf C} = \sigma_{meas}^2 {\bf XXI}
\end{equation}

\noindent because ${\bf C}$ is defined for $\chi^2$ instead of
$\sigma^2$ (see \S \ref{chisqsection}). 

	In ${\bf 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 \rightarrow 1$. Each matrix element is the correlation coefficient
between the {\it 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 {\it true} values $(A_0^*, A_1^*)$ by
amounts $\delta A_0$ and $\delta A_1$.  If their covariance is zero,
then in the infinity of data sets you'd find that $\delta A_0$ is
uncorrelated with $\delta 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 \S \ref{chicoeffs}, 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 {\it 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.

\subsection{ The covariance in our numerical example} \label{covarianceone}

	Apply equation \ref{idlcovariance} to obtain the covariance
matrix for our numerical example:

\begin{eqnarray}
{\bf ncov} = \left[
\begin{array}{ccc}
1  &  -.989848   &  .969717  \\
-.989848  &  1   &  -.993808  \\
.969717  & -.993808 &  1  \\
\end{array} \; \right] \; .
\end{eqnarray}

	The off-diagonal elements are {\it 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 \ref{lsfitfig} 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

\begin{mathletters}
\begin{eqnarray}
{\bf A} = \left[
\begin{array}{c}
188.625 \\
18.500 \\
0.87500 \\
\end{array} \; \right] 
\end{eqnarray}

\begin{eqnarray}
{\bf \sigma_A} = \left[
\begin{array}{c}
3.58 \\
1.00 \\
0.559 \\
\end{array} \; \right] \; .
\end{eqnarray}
\end{mathletters}

\noindent 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

\begin{eqnarray}
{\bf ncov} = \left[
\begin{array}{ccc}
1  &     0       &   -0.78086881  \\
   0  &  1   &    0  \\
 -0.78086881 &  0  &  1  \\
\end{array} \; \right] \; .
\end{eqnarray}

\noindent which is nice, but not perfect: Our step is {\it partial}
because the second-order coefficient $A_2$ affects $A_0$ because, over
the range of $[(time - 8) = -3 \rightarrow +3]$, the quantity $[A_2 \;
\Sigma(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 \rightarrow 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 {\it 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 \S \ref {chicoeffs}.
Also, you might also want to try out another example in Taylor's \S 8.5.

\section{REJECTING BAD DATAPOINTS: CHAUVENET'S CRITERION}
\label{chauvenetsection}

	Least squares fitting is derived from the maximum likelihood
argument assuming the datapoint residuals $\delta y_m$ have a Gaussian
pdf. This means that the errors are distributed as 

\begin{equation}
p_{(\delta y; \sigma)} = {1 \over \sqrt{2\pi} \sigma} 
    e^{-\left( \delta y^2 \over 2 \sigma^2 \right)}
\end{equation}

\noindent where $\sigma^2$ is the true variance of the datapoints, i.e.\
$s^2$ in equation \ref{samplevarianceone} (to be precise, $s^2$ needs to
be averaged over many experiments). 

	More importantly, the probability of finding datapoints inside
the limits $\pm \Delta y$ is 

\begin{equation}
P_{(|\delta y| < \Delta y)} = \int_{-\Delta y}^{+\Delta y} 
  p_{(\delta y; \sigma)} d(\delta y) = 
  {\rm erf}\left({\Delta y \over \sqrt{2} \sigma} \right)
\end{equation}

\noindent where we use the commonly-defined error function ${\rm
erf}_{(X)} = {1 \over \sqrt{\pi}} \int_{-X}^{+X} e^{-x^2} dx$. A
particularly important value is for $\Delta y = \sigma$, for which 

\begin{equation}
P_{(|\delta y| < \sigma)} = 0.683
\end{equation}

	If we have an experiment with $M$ datapoints, then the number of
datapoints we expect to lie outside the interval $\pm \Delta y$ is

\begin{equation}
M_{(outside \ \Delta y)} = M \left[ 1 - 
  {\rm erf} \left({\Delta y \over \sqrt{2} \sigma} \right) \right]
\end{equation}

\noindent Chauvenet's criterion simply says: \begin{enumerate}

	\item Find $\Delta y$ such that $M_{(outside \ \Delta y)} =
0.5$. This is given by 

\begin{equation}
{\Delta y \over \sigma} = {\rm erf}^{-1}\left( {1 - {1 \over 2M}} \right)
\end{equation}

	\item Discard all datapoints outside this range.

\end{enumerate}

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

\begin{table} [!h]
\begin{center}
%\hline
%\caption{ Chauvenet's criterion versus $M$} \label{chauvenettable}
\begin{tabular}{cc} 
%\hline
\\ \hline \hline
\multicolumn{2}{c}{ Chauvenet's criterion versus $M$} 
\\ \hline
$M$ & $\Delta y \over \sigma$ \\ \hline
100 & 2.81 \\
1000 & 3.48 \\
$10^4$ & 4.06 \\
$10^5$ & 4.56 \\ \hline
\hline
\end{tabular}
\end{center}
\end{table}

\noindent which is a moderately interesting set of numbers.

	We offer the following important {\it Comments}:
\begin{itemize}

	\item 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 {\it outliers}.

	It is very important to get rid of the outliers because the
least squares process minimizes the {\it 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 {\it intrinsic} pdf, then you should think twice before using
least squares!

	\item You may wish to relax Chauvenet's criterion by {\it
increasing} the $\Delta 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 {\it never} make Chauvenet's criterion more stringent
by {\it decreasing} the $\Delta x$ beyond which you discard points. This
rule hardly needs elaboration: it means you are discarding datapoints
that follow the assumed pdf!

	\item 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 {\it only once} and certainly not
repeatedly.

	Being a real-life astronomer, our approach is different. There
{\it do} exist outliers. They increase the calculated value of $\sigma$.
When you discard them, you are left with a more nearly perfect
approximation to Gaussian statistics and the new $\sigma$ calculated
therefrom will be smaller than when including the outliers. Because the
original $\sigma$ 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.

	\item Many observers use the $3\sigma$ criterion: discard any
points with residuals exceeding $3\sigma$. This is definitely {\it not}
a good idea: the limit $3\sigma$ is Chauvenet's criterion for $M=185$
datapoints. Very often $M$ exceeds this, often by a lot. 

	\item 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 {\bf 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 \S 9.4 and 9.6)
together with the error function; both procedures exist in IDL as {\bf
newton} and {\bf errorf}. You at least ought to skim lightly some of
NM's chapter 9 about root finding, because some day you'll need it.

\end{itemize}

\section{NONLINEAR LEAST SQUARES} \label{nonlinearls}

	The least squares formulation requires that the data values
depend {\it linearly} on the unknown coefficients. For example, in
equation \ref{one}, 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 

\begin{equation} \label{basicone}
\sin(A t_m) + B t_m = \theta_m \; .
\end{equation}

\noindent 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 {\it
predicted} values $\theta_{0,m}$:

\begin{equation} \label{trialone}
\sin(A_0 t_m) + B_0 t_m = \theta_{0,m} \; .
\end{equation}

\noindent Subtract equation \ref{trialone} from \ref{basicone}, and
express the differences as derivatives. Letting $\delta A = A - A_0$ and
$\delta B = B - B_0$, this gives

\begin{equation} \label{difference}
\delta A [t_m \cos(A_0 t_m)] + \delta B t_m = \theta_m - \theta_{0,m} \; .
\end{equation}

\noindent This is linear in $(\delta A, \delta B)$ so you can solve for
them using standard least squares. Increment the original guessed values
to calculate $A_{0,new}=A_0 + \delta A$ and $B_{0,new}=B_0 + \delta 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.

	{\bf (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. 

	{\bf (1): Multiple minima:} Nonlinear problems often have
multiple minima in $\sigma^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 $\sigma^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 $\sigma^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 $\sigma^2$ weren't the true minimum.

	{\bf (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 {\it full} confidence) that the solution is unique.

	{\bf (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 ($\delta A$ and $\delta B$ above) by a factor
${\cal F}$ less than unity, for example ${\cal F} = 0.5$ or 0.75. This
increases the number of iterations required for convergence but often
allows convergence instead of producing instability.

	{\bf (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.~$\delta A \over 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 ${\cal F} \neq 1$, set
${\cal F} = 1$, do a few more iterations, and you're done. 

	{\bf (5): Numerical derivatives:} Sometimes the equations of
condition are so complicated that taking the derivatives, as in
obtaining equation \ref{difference}, 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.

	{\bf (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 ${\cal 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 \S 15.5); IDL function {\bf 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. 

	{\bf (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!

        {\bf (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 $\phi$
in the function $y = \cos( \theta + \phi)$.  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 \theta - B \sin \theta$,
where ${B \over A} = \tan \phi$.  Solve for $(A,B)$ using linear least
squares, calculate $\phi$, and propagate the uncertainties. 

        An example of the not-so-harmless case is in NM's \S 15.4
example: fit for $(A,B)$ with equations of condition $y_m = A
e^{-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
{\it observed} values $y_m$; thus the associated errors
$\sigma_{meas,m}$ are {\it also} affected.  This means you have to do
weighted fitting, which is discussed in \S \ref{chisqsection} below. 
Suppose that $A=1$, your datapoints all have $\sigma_{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 $\sigma_{meas_m}$, but what is the corresponding
value of $\sigma_{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 {\it negative}??? Even for $y_m$ not near zero,
you need to calculate new $\sigma_{meas,m}$ by error propagation; in
this case, you need to reassign $\sigma(\log y) = {d \log y \over dy}
\sigma(y) = {\sigma(y) \over 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.

\section{CHI-SQUARE FITTING, AND WEIGHTED FITTING: DISCUSSION IGNORING
COVARIANCE } \label{chisqsection}

	In least squares fitting the derived parameters minimize the sum of
squares of residuals as in equation \ref{samplevarianceone}, which we
repeat here

$$ s^2 = {1 \over M - N} \sum_{m=0}^{M-1} (\theta_m -
\overline{\theta_m})^2 \; . $$

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

\begin{mathletters}
\begin{equation}
 \chi^2 = \sum_{m=0}^{M-1} {(\theta_m - \overline{\theta_m})^2
     \over \sigma_{meas,m}^2 } \; . 
\end{equation}

\noindent Along with $\chi^2$ goes the {\it reduced} chi squared
$\overline{\chi}^2 = {\chi^2 \over M-N}$

\begin{equation}
\overline{\chi}^2 = {1 \over M - N} \sum_{m=0}^{M-1} {(\theta_m - \overline{\theta_m})^2
     \over \sigma_{meas,m}^2 } \; . 
\end{equation}
\end{mathletters}

\noindent 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 $\sigma_{meas,m}$. Thus, the reduced Chi-square
($\widehat{\chi}^2$) is equal to the ratio of the {\it variance of the
datapoint residuals} ($\sigma^2$) to the {\it adopted intrinsic
measurement variances} ($\sigma_{meas}^2$).  So it should be obvious
that in Chi-square fitting, you must know the measurement uncertainties
$\sigma_{meas}$ of the individual data points beforehand.  If you want
to give the various datapoints weights based on something other than
$\sigma_{meas}$, then that is just like Chi-square fitting except that
you can adopt an arbitrary scale factor for the uncertainties
(section~\ref{diatribe}). 

	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 \begin{enumerate}

\item With the standard treatments, the errors in the derived parameters
don't depend on the residuals of the datapoints from the fit (!).

\item The errors in the derived parameters can depend on their mutual
covariances. This discussion requires a separate section, which we
provide below in \S \ref{chicoeffs}. 

\end{enumerate}

\noindent 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. 

\subsection{ The weighted mean: the simplest chi-squared fit}

\label{chicoeffsnoco}

        First, recall the formulas for an ordinary {\it unweighted}  
average in which the value of each point is $y_m$ and the residual of
each point from the weighted mean is $\delta y_m$:

\begin{mathletters}
\begin{equation}
mean = {\sum y_m \over M}
\end{equation}
\begin{equation}
%\sigma^2 = {\sum \delta y_m^2 \over M-1}
s^2 = {\sum \delta y_m^2 \over M-1}
\end{equation}
\begin{equation} \label{meanvariance}
%\sigma_{mean}^2  = {\sigma^2 \over M} = {\sum \delta y_m^2 \over M(M-1)} \ ,
s_{mean}^2  = {s^2 \over M} = {\sum \delta y_m^2 \over M(M-1)} \ ,
\end{equation}
\end{mathletters}

\noindent $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} = {1 \over \sigma_{meas,m}^2} = {1 \over \sigma_{meas}^2}$.
Applying maximum likelihood, in a nonweighted average the
quantity that is minimized is $\sum \delta y_m^2$; in a weighted average
the quantity minimized is $\chi^2 = \sum {\delta y_m^2 \over
\sigma_{meas,m}^2} = \sum w_{meas,m} \delta y_m^2$.  So you'd
think that the three equations corresponding to the above would become

\begin{mathletters} \label{wgtideal}
\begin{equation}   
mean_w = {\sum w_{meas,m}y_m \over \sum w_{meas,m}} \rightarrow { \sum y_m \over M}
\end{equation}
\begin{equation} \label{wvariance}
%\sigma_w^2 = 
s_w^2 = 
  {M \over M-1} {\sum w_{meas,m} \delta y_m^2 \over \sum w_{meas,m}} 
    \rightarrow {\sum \delta y_m^2 \over M-1}
\end{equation}

\begin{equation} \label{wmeanvariancecc}
%\sigma_{mean}^2 = {\sigma_w^2 \over M} = 
s_{mean}^2 = {s_w^2 \over M} = 
  {\sum w_{meas,m} \delta y_m^2 \over (M-1) \sum w_{meas,m}} \rightarrow
   {\sum {\delta y_m^2} \over M(M-1) } \ ,
\end{equation}   
\end{mathletters}

\noindent where the arrow to the rightmost expressions use our
assumption that  $\sigma_{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)
        
\begin{equation} \label{wmeanvariancebb}
%\sigma_{mean}^2 = {1 \over \sum w_{meas,m}} \rightarrow { \sigma_{meas}^2 \over M} \ ,
s_{mean}^2 = {1 \over \sum w_{meas,m}} \rightarrow { \sigma_{meas}^2 \over M} \ ,
\end{equation}

        Note the excruciatingly painful difference between equation
\ref{wgtideal}c and equation \ref{wmeanvariancebb}: the former depends
on the {\it variance of the datapoint residuals} $s_w^2$, as you'd
think it should, while the latter depends on only the {\it adopted
intrinsic measurement variance of the data} $\sigma_{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 $\sigma_{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 \ref{wmeanvariancebb}, the far superior

\begin{equation} \label{wmeanvarianceaa}
%\sigma_{mean}^2 = {\widehat{\chi}^2 \over \sum w_{meas,m}} =
s_{mean}^2 = {\widehat{\chi}^2 \over \sum w_{meas,m}} =
    {\chi^2 \over (M-1) \sum w_{meas,m}} \rightarrow
%   { \sigma_{w}^2 \over M} \ ,
   { s_{w}^2 \over M} \ ,
\end{equation}

\noindent which is precisely the same as our intuitive guess, equation
\ref{wmeanvariancecc}.  The difference between equations
\ref{wmeanvarianceaa} and \ref{wmeanvariancebb} is the numerator, which
contains the reduced chi-squared $\widehat{\chi}^2$; for this case where
all $\sigma_{meas,m}$ are identical, $\widehat{\chi}^2 = {s_w^2
\over \sigma_{meas}^2}$.  Here

\begin{equation} \label{chisqone}
\chi^2 = \sum w_{meas,m} \delta y_m^2 = \sum {\delta y_m^2 \over \sigma_{meas,m}^2}
\end{equation}

\subsection{The multivariate Chi-square fit}
\label{sectionchisq}

	In this case, Chi-square fitting is just like least-squares
fitting except for the following: \begin{enumerate}

	\item In the matrix $\bf X_\sigma$ of equation
\ref{matrixdefinition}a, each row $m$ is a different measurement with a
different intrinsic variance $\sigma_{meas,m}$.  You are generating a
new matrix $\mathbf {X}$, which is identical to ${\bf X_\sigma}$ except
that each row $m$ is divided by $\sigma_{meas,m}$.  This new matrix is
the same as NR's {\it design matrix} (Figure 15.4.1), which they denote
by $\bf A$. 

	\item Divide each datapoint $\theta_m$ in equation
\ref{matrixdefinition}b by $\sigma_{meas,m}$. You are generating a new
matrix $\mathbf {Y}$, which is identical to ${\bf Y_\sigma}$ except that
each row is divided by $\sigma_{meas}$.  This new matrix is the same as
NR's vector {\bf b}. 
\end{enumerate}

\noindent 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 $\sigma_{meas}$.  In the full matrix
equation~\ref{equationofcondition}, which now reads ${\bf X \cdot a} =
{\bf Y}$, you've divided both sides by the same quantity,
$\sigma_{meas}$.  Obviously, this doesn't change the derived results for
{\bf a}. (If the $\sigma_{meas,m}$ differ from each other, then $\bf a$
is affected, of course.) 

	Let's summarize the relationships between the least-squares
($\sigma$-subscripted) and $\chi^2$ matrices for this specific case
where all $\sigma_{meas,m}$ are identical:

\begin{mathletters}
\begin{equation}
{\bf X} =  { {\bf X_\sigma} \over \sigma_{meas} }
\end{equation}

\begin{equation}
{\bf Y} =  { {\bf Y_\sigma} \over \sigma_{meas} }
\end{equation}

\begin{equation}
{\bf X \cdot a} = {\bf Y}
\end{equation}

\begin{equation} \label{curvmatrixalpha}
[\alpha]  =  {\bf X^T \cdot X} = { [\alpha_\sigma] \over \sigma_{meas}^2 }
\end{equation}

\begin{equation}
[\beta] = {\bf X^T \cdot Y} = { [\beta_\sigma] \over \sigma_{meas} }
\end{equation}

\begin{equation}
[\alpha]^{-1}  =  \sigma_{meas}^2  [\alpha_\sigma]^{-1} 
\end{equation}

\begin{equation}
{\bf a} = [\alpha]^{-1} \cdot [\beta]
\end{equation}

\end{mathletters}

	When you calculate the predicted values ${\bf \overline{ Y}}$ in
equation \ref{ybardefinition}, using $\bf X$ instead of $\bf X_\sigma$,
the values are divided by $\sigma_{meas}$. Because $\bf Y$ is also
divided by $\sigma_{meas}$, the residuals $\bf \delta Y$ are as well.
Thus when you calculate the square of the sum of the residuals as in
equation \ref{samplevariancetwo} you get $s^2 \over \sigma_{meas}^2$,
which is the reduced chi-square, denoted as $\widehat{\chi}^2$ :

\begin{equation}
\label{reducedchisq}
\widehat{\chi}^2 = { {\bf \delta Y^T \cdot \delta Y} \over M-N}
                       = { {\bf \delta Y_\sigma^T \cdot \delta Y_\sigma} \over \sigma_{meas}^2(M-N)}
\end{equation}

\noindent and, obviously, you could calculate the non-reduced $\chi^2$

\begin{equation}
\label{chisq}
\chi^2 =  {\bf \delta Y^T \cdot \delta Y}
  = { {\bf \delta Y_\sigma^T \cdot \delta Y_\sigma} \over \sigma_{meas}^2}
\end{equation}

\noindent This is precisely analogous to equation \ref{chisqone} for the
weighted mean example above. The reduced chi-square $\widehat{\chi}^2$
is equal to the ordinary chi-square $\chi^2$ except that it is divided
by the number of degrees of freedom, ordinarily denoted by $\nu$, which
is equal to $(M-N)$. 

	Finally, we have the analogy of equation \ref{wmeanvarianceaa}
expressed in matrix form as in equation \ref{coeffvarianceone}:

\begin{equation} \label{coeffvariancegood}
{\bf s_a}^2 = \widehat{\chi}^2 diag \{[\alpha]^{-1} \} 
\end{equation}

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

\begin{equation} \label{coeffvarianceawful}
{\bf s_a}^2 = diag \{[\alpha]^{-1} \}
\end{equation}

\noindent 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 \ref{coeffvariancegood}, 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 \approx
\sigma_{meas}^2$ and $\widehat{\chi}^2 \approx 1$. (We write
``$\approx$'' instead of ``='' because different experiments produce
somewhat different values of $s^2$ because of statistical fluctuations;
an average gives $\sigma^2 = \langle s^2 \rangle$). In this situation,
though, equations \ref{coeffvariancegood} and \ref{coeffvarianceawful}
are identical.

\subsection{ Persnickety Diatribe on Choosing $\sigma_{meas}$ }
\label{diatribe} 

\subsubsection{Choosing and correcting $\sigma_{data,m}$}

	In the previous section, equation \ref{coeffvarianceawful}
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 $\sigma_{meas,m}$ and not on the {\it actual
variance} of the {\it datapoints}. 

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

        Formally, $\sigma_{x0}$ depends only on the {\it adopted}
uncertainties $\sigma_{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 $\sigma_{meas,m}$ are off by a
common scale factor, i.e.\ if $\sigma_{meas,adopted} =
f\sigma_{meas,true}$.  Then $\widehat{\chi}^2 \approx f^{-2}$ instead of
$\widehat{\chi}^2 \approx 1$.  And to obtain the parameter errors from
$\delta \chi^2$, you must find the offset $\delta x$ such that $\Delta
\chi^2 = f^{-2} \approx \widehat{\chi}^2$. 

	You can correct for this erroneous common factor $f$ by dividing
your adopted values of $\sigma_{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 $\widehat{\chi}$.  And,
of course, the same as multiplying $\sigma_{meas,m}^2$ by
$\widehat{\chi}^2$. 

	To be kosher, after having run through the problem once with the
adopted $\sigma_{meas,m}$, calculate the $\widehat{\chi}^2$; multiply
all $\sigma_{meas,m}$ by $\widehat{\chi}$; and redo the problem so that
the new $\widehat{\chi}^2 = 1$.  Then the derived variance
$\sigma_{x_0}^2$ is also correct.  You can obtain it either as the
corresponding diagonal to the covariance matrix (equations
\ref{coeffvariancegood} and \ref{coeffvarianceawful}, which are
identical in this case) or by finding what departure from $x_0$ is
necessary to make $\Delta \chi^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.\footnote{The $\Delta \chi^2 = 1$ condition only applies for
the single-variable case.  For more variables, the variances need to
make $\Delta \chi^2 > 1$; for example, if you are fitting for two
parameters, then the variances make $\Delta \chi^2 = 2.3$, as we
discussed in \S \ref{chicoeffs}.}

\subsubsection{Come back to reality}

	In the case $\widehat{\chi}^2 \approx 1$ the dispersion of the
observed points $\sigma$ is equal to the intrinsic dispersion of the
datapoints $\sigma_{meas}$ and the mathematical model embodied in the
least-squares fit is perfect.  That, at least, is the {\it theoretical}
conclusion.  In practice, however, your obtaining such a low, good value
for $\widehat{\chi}^2$ might mean instead that you are using {\it too
large} a value for $\sigma_{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 $\sigma_{meas}$ and you specify that they have $2
\sigma_{meas}$ when computing $(\mathbf{X},\mathbf{Y})$, then you'll
find $\widehat{\chi}^2 = 0.25$. 

        There are, in fact, two ways you can get artificially small
values for $\widehat{\chi}^2$. One is by overestimating your adopted
intrinsic values of $\sigma_{meas}$. The other is if your {\it
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.

	High values of $\widehat{\chi}^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 {\it less} error to your datapoints than they really have. And
in this case, using equation \ref{coeffvarianceawful} instead of
\ref{coeffvariancegood} is disastrous.

	Think about it.

\subsection{The case in which datapoints have different dispersions,
or different weights: like a weighted average}

\label{weightedcase}

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

	Consider two measurements with intrinsic measurement
uncertainties $(\sigma_{meas,1}, \sigma_{meas,2})$; suppose
$\sigma_{meas,1} < \sigma_{meas,2}$.  After being divided by their
respective $\sigma_{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 ${\sigma_{meas,1} \over
\sigma_{meas,2}} = \left({\tau_1 \over \tau_2}\right)^{-1/2}$, so the
rows of $\mathbf {X}$ are weighted as $\tau^{1/2}$.  This means that
during the computation of $[\alpha] = \mathbf {X^T \cdot X}$, the
self-products of row 1 are weighted as $\tau_1$.  This means that each
datapoint is weighted as $\tau$, 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 $\left(1
\over \sigma_{meas}\right)^2$.  We conclude that weighting scheme of the
first two steps in section~\ref{sectionchisq} agrees with common sense. 

	Suppose you don't know the intrinsic measurement dispersion
$\sigma_{meas}$, but you {\it do} know the {\it 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 $\sigma_{meas} \propto \tau^{-1/2}$.  In this case, multiply
each row by its weight $w \propto {1 \over \sigma_{meas}}$ and proceed
as above. 

\subsection{The explicit IDL for the weighted and chi square cases}
\label{idlexplicitone}

	Here we give the explicit IDL commands required to do chi
squared fitting, as we did for least squares fitting in \S
\ref{idlsolution}.  First, define the {\bf X} and {\bf Y} arrays as
before in \S \ref{numexample}.  However, here we need to weight the
matrices.  Let us denote the weighted versions by {\bf XW} and {\bf YW}.
 Then each row of {\bf XW} and element of {\bf YW} are divided by the
corresponding $\sigma_{meas}$.  Let us define {\bf W} as the diagonal
matrix\footnote{ {\bf W} has $M \times M$ elements but only $M$ numbers.
If $M$ is large and you are impacting machine memory, then you can save
memory by: reformulating equations \ref{xwwx} and \ref{bteqnww} to use
loops; not keeping {\bf XW} and {\bf X} as separate matrices; and using
IDL's {\bf temporary} function wherever possible.} whose elements are
equal to $1 \over \sigma_{meas,m}$, i.e. $w_m = {1. \over
\sigma_{meas,m}}$

\begin{mathletters}
\begin{equation}
{\bf W} = fltarr( M, M)
\end{equation}
\begin{equation}
{\bf W[ indgen(M) * (M+1)]} = \left[ { 1 \over \sigma_{meas,m} } \right]
\end{equation}
\end{mathletters}

\noindent Note that $w_m$ is the {\it square root} of the weights
$w_{meas,m}$ that we discussed in \S \ref{chicoeffsnoco}.  Specifically,
the current weights $w_m$ are {\it not} those in the sense of the
weighted average of \S \ref{chicoeffsnoco}; we use $w_m$ as a matter of
convenience only. Then

\begin{mathletters} \label{xwwx}
\begin{equation}
{\bf XW} = {\bf W \ \#\# \ X}
\end{equation}
\begin{equation}
{\bf YW} = {\bf W \ \#\# \ Y}
\end{equation}
\end{mathletters}

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

	All steps up to obtaining the uncertainties follow exactly the
same steps as in \S \ref{idlsolution}, which we repeat here.  In IDL we
calculate the normal equation matrices and denote the $[\alpha]$ in
equation \ref{ssdef}a by ${\bf XXW}$:

\begin{mathletters} \label{eqn60}
\begin{equation}
{\bf XXW} = {\bf transpose(XW) \# \# XW} \; ,
\end{equation}

\noindent and we denote the $[\beta]$ in equation \ref{ssdef}b by {\bf
XYW}:

\begin{equation}
{\bf XYW} = {\bf transpose(XW) \# \# YW} \; .
\end{equation}
\end{mathletters}

\noindent In IDL we take the inverse of $[\alpha]$ (same as $\bf XXW$) by

\begin{equation}
{\bf XXWI} = {\bf invert(XXW)} \; .
\end{equation}

        The least-square fitted quantities are in the matrix {\bf a}
(equation \ref{adef}), which we obtain in IDL with

\begin{equation}
{\bf a} = {\bf XXWI \ \#\# \ XYW} \; .
\end{equation}

        In IDL we denote the matrix of {\it weighted} predicted values
by ${\bf YBARW}$, which is

\begin{equation}
\label{bteqnw}
{\bf YBARW} = {\bf XW \ \#\# \ a} \; .
\end{equation}

\noindent and we can also define the {\it weighted} residuals as

\begin{equation}
{\bf DELYW} = {\bf YW - YBARW}
\end{equation}

\noindent The {\it true} (unweighted) predicted values and residuals are

\begin{mathletters}
\label{bteqnww}
\begin{equation}
{\bf YBAR} = {\bf invert(W) \ \#\# XW \ \#\# \ a} = {\bf X \ \#\# \ a} \; .
\end{equation}

\begin{equation}
{\bf DELY} = {\bf invert(W)\ \# \# \ (YW - YBARW)} = {\bf Y - YBAR }
\end{equation}
\end{mathletters}

\noindent To calculate the chi square we want the weighted sum $\sum
\delta w_m^2 y_m^2$:

\begin{mathletters}
\begin{equation}
\label{sigsqw}   
chi\_sq = {\bf transpose(DELYW) \ \#\# \ DELYW} \; .
\end{equation}

\noindent and the reduced chi square

\begin{equation}
\label{sigsqww}   
red\_chi\_sq = {\bf transpose(DELYW) \ \#\# \ DELYW}/(M-N) \; .
\end{equation}
\end{mathletters}

\noindent The full covariance matrix is just {\bf XXWI} and,
corresponding to equation \ref{lscoefferrorone}, 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 $\sigma_{meas,m}$
and/or the quality with which the model fits the data.  That is, 

\begin{mathletters} \label{zyx}
\begin{equation} 
{\bf vardc\_ official} = {\bf XXWI[(M+1)*indgen(M)]} 
\end{equation}

\noindent or, better, 
\begin{equation} 
{\bf vardc\_ realworld} = red\_chi\_sq * {\bf XXWI[(M+1)*indgen(M)]} 
\end{equation}
\end{mathletters}

\noindent Finally, the normalized covariance matrix is of course
independent of chi squared and is the same as equation
\ref{idlcovariance}

\begin{equation} \label{idlcovariancew}
{\bf ncov = XXWI/{\it sqrt}( vardc\_official \# vardc\_official)}
\end{equation}

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

\section{CHI-SQUARE FITTING, AND WEIGHTED FITTING: DISCUSSION INCLUDING
COVARIANCE } \label{chicoeffs}

\subsection{ Phenomenological description}

	Consider the first two coefficients in our above example, which
we discussed a bit above in \S \ref{covarianceone}.  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~\ref{avalues}. The coefficient
$A_0$ is the $y$-intercept and $A_1$ is the slope. They have derived
values $A_0 = 96 \pm 34$ and $A_1 = 4 \pm 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$ {\it without
regard to $A_1$}. By this we mean that as $A_0$ is varied from its
optimum value of 96, $\chi^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 $\chi^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 \ref{coeffvarianceone}, and \ref{coeffvariancegood} apply. 

	However, if you are interested in knowing about {\it 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 {\it individual} values $\delta A_0=+34$ and $\delta
A_1=+9$ are acceptable, you {\it cannot} conclude that the {\it pair} of
values $(\delta A_0,\delta A_1) = (+34, +9)$ is acceptable, because this
pair has both positive.  In contrast, what {\it is} acceptable here
would be something like $(\delta A_0,\delta 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 \pm 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 $(\infty \rightarrow -\infty)$: you just don't give a
damn {\it 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 $\Delta \chi^2$ and present them on a graph with axes $(\delta
a_0, \delta 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. 

\begin{figure}[h!] 
\begin{center} 
\leavevmode
%\includegraphics[height=7.5in, width=6.0in]{twolvl3b.ps}
%\includegraphics[scale=.55, angle=90]{lsfitfig.ps}
\includegraphics{chisq.ps} 
\end{center} 
\caption{Illustrating
the meaning of variance and covariance between $(a1,a2)$ for our numerical
example. See text for discussion. \label{figchisq}}
\end{figure}

	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 \S \ref{easycalc} and are at
$\Delta \chi^2 = 1$ and 2.3. The dashed horizontal and vertical lines
are at $\delta_a = \pm \sigma_a$.

	First consider the pair of vertical lines, which are drawn at
$\delta a_1 = \pm \sigma_{a_1}$, where $\sigma$ is the square root of
the variance of the parameters as described in equations
\ref{coeffvarianceone}, \ref{lscoefferrorone},
\ref{coeffvariancegood}, \ref{coeffvarianceawful}, and \ref{zyx}.  If
the datapoints were projected downward, i.e.\ if we take small strips
$\delta a_1$ and integrate over $\delta a_2$, the pdf of $\delta 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 $\chi^2$ as we consider departures
$\delta a_1$, then the pdf of $\delta a_1$ has dispersion
$\sigma_{a_1}$. Alternatively, we can say that in a large number of
experiments, the pdf of $\delta a_1$ follows a chi squared pdf with one
degree of freedom if we don't care what happens to $\delta 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 $\nu=2$,
$68.3\%$ of the points lie within $\Delta \chi^2=2.3$, where we have
drawn the outer contour in Figure \ref{figchisq}. 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 \S11.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 ${\bf a} = [a_0, a_1, a_2]$. We'll first begin by
discussing the case of single parameter; then we'll generalize.

\subsection{ Calculating the uncertainties of a single parameter}
\label{gencalc0} 

	First, suppose we want to know the value $\sigma_{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 $\delta a_0$ around this best value.  

	Pick a value of $\delta a_0$ and redo the least squares solution
for $[a_1, a_2]$.  This gives a new value for $\chi^2$
which is, of course, larger than the minimum value that was obtained
with $\delta a_0=0$.  Call this difference $\Delta \chi_{\delta a_0}^2$. 
Determine the dependence of $\Delta \chi_{\delta a_0}^2$ upon $\delta
a_0$ and find the value of $\delta a_0$ such that $\Delta \chi_{\delta
a_0}^2 = 1$.  This is the desired result, namely the value
$\sigma_{a_0}$ without regard to the values of $a_1$ and $a_2$. 

	This value is $\sigma_{a_0}^2 = [\alpha]_{00}^{-1}$, the same
result quoted in equation \ref {coeffvarianceawful}.

	Consider now what you've done in this process.  For each least
squares fit you used a trial value of $\delta a_0$.  In specifying
$\delta 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 $\Delta \chi_{\delta a_0}^2$.  It should be
clear that this distribution follows a chi squared distribution with one
degree of freedom ($\nu = 1$).  So the uncertainty $\sigma_{a_0}$ is
that value for which $\Delta \chi_{\delta 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.)

\subsection{ Calculating the uncertainties of two parameter}
\label{gencalc1} 

	Suppose we want to know the value $(\sigma_{a_0},\sigma_{a_2})$
without regard to the value of $a_1$.  Now we consider variations
$(\delta a_0, \delta a_2)$ around the best values $(a_0, a_2)$. 

	Pick values for $(\delta a_0, \delta a_2)$ and redo the least
squares solution for $a_1$.  This gives a new value for $\chi^2$ which
is, of course, larger than the minimum value that was obtained with
$(\delta a_0,\delta a_2)=0$.  Call this difference $\Delta \chi_{(\delta
a_0,\delta a_2)}^2$.  As above, this follows a chi square distribution,
but now with $\nu = 2$. Determine the dependence of $\Delta \chi_{(\delta
a_0,\delta a_2)}^2$ upon $(\delta a_0, \delta a_2)$ and find the set of
values of $(\delta a_0, \delta a_2)$ such that $\Delta \chi_{(\delta
a_0, \delta a_2)}^2 = 2.3$.  This is the desired result, namely the
ellipse within which the actual values $(\delta a_0, \delta 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
$[\alpha]$, as we discuss below.  

	Consider now what you've done in this process.  For each least
squares fit you used trial values of $(\delta a_0,\delta 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 ($\nu = 2$).  So the uncertainty
$\sigma_{a_0}$ is that value for which $\Delta \chi_{\delta a_0}^2 =
2.3$, which follows from the integrated probability for the chi squared
distribution for $\nu=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\dots

\subsection{ Calculating the uncertainties of three parameter}
\label{gencalc2} 

	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 $\Delta \chi_{\bf a}$, which defines a 3- (or $N$-)
dimensional ellipsoid.  This follows a chi squared distribution with
$\nu = 3$ (or $N$).  We find the (hyper)surface such that $\Delta
\chi_{\bf a}$ is that value within which the integrated probability is
$68.3\%$.  This defines the (hyper)surface of $\sigma_{\bf a}$. 

\subsection{ Doing these calculations the easy way} \label{easycalc}

	The obvious way to do the calculations described above is to set
up a grid of values in the parameters of interest $(\delta a_n)$;
perform the chi squared fit on the remaining variables, keeping track of
the resulting grid of $\chi^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
($[\alpha]$ of equation \ref{curvmatrixalpha}, the same as {\bf XXW})
contains the matrix of the second derivatives of $\chi^2$ with respect
to all pairwise combinations of $\delta a_n$, evaluated at the minimum
$\chi^2$; it's known as the curvature matrix for this reason.  Clearly,
as long as the Taylor expansion is good we can write

\begin{equation}
\Delta \chi_{\bf a}^2 = {\bf \delta a^T \cdot [\alpha] \cdot \delta a}
\end{equation}

\noindent 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 \S15.6 (Probability
Distribution of Parameters in the Normal Case)]. Follow the steps
\begin{enumerate}

	\item Decide which set of parameters you are interested in; call
this number $N_i$ and denote their vector by ${\bf a_i}$.  Here we use
the above example and consider $N_i=2$ and ${\bf a_i} = [a_0,a_2]$. 

	\item From the $N \times N$ covariance matrix $[\alpha]^{-1}$,
extract the rows and columns corresponding to the $N_i$ parameters and
form a new $N_i \times N_i$ covariance matrix $[\alpha]_i^{-1}$; in our
case the original covariance matrix is 

\begin{mathletters}
\begin{eqnarray}
[\alpha]^{-1} = {\bf XXI} = \left[
\begin{array}{ccc}
    1156.8125  &   -303.000   &   18.4375  \\
     -303.000 &     81.000  &   -5.000  \\
      18.4375  &   -5.000   &  0.31250  \\
\end{array} \; \right] 
\end{eqnarray}

\noindent and it becomes

\begin{eqnarray}
[\alpha]_i^{-1} = {\bf XXI} = \left[
\begin{array}{cc}
    1156.8125  &   18.4375  \\
      18.4375  &  0.31250  \\
\end{array} \; \right] \ .
\end{eqnarray}
\end{mathletters}

	\item Invert this new covariance matrix to form a new curvature
matrix $[\alpha]_i$. The elements differ from the those in the original
curvature matrix. 

	\item As usual, we have

\begin{equation}
\Delta \chi_{\bf a_i}^2 = {\bf \delta a_i^T \cdot [\alpha]_i \cdot \delta a_i}
\end{equation}

\noindent so find the locus of ${\bf a_i}$ such that the integrated
probability of $\Delta \chi_{\bf a_i}$ for $\nu = N_i$ contains $68.3\%$
of the space; e.g.\ for $\nu=2$ this is $\Delta \chi_{\bf a_i} = 2.3$. 

\end{enumerate}

	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. 


\subsection{ Important comments about uncertainties}

	Having said all the above, we offer the following important {\it
Comments}: \begin{itemize}

	\item The easiest way to calculate these (hyper)surfaces is to
set up a grid in $N_i$-dimensional space of trial values for $\delta
{\bf a_i}$ and use a contour plot or volume plot package to plot the
loci of constant $\Delta \chi_{\bf a_i}^2$.

	\item The procedure described in \S \ref{easycalc} works well
for linear fits, or nonlinear fits in which the $\sigma_{\bf a}$ are
small so that $\Delta \chi^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 $\bf a_i$ as described above in \S
\ref{gencalc0}, \ref{gencalc1}, and \ref{gencalc2}.

	\item The variance (i.e., uncertainty squared) of the derived
parameters {\bf a} depends {\it only on the elements in the covariance
matrix $[\alpha]^{-1}$}.  These, in turn, depend only on the curvature
matrix $[\alpha]$.  These, in turn, depend only on ${\bf X}$.  This
matrix {\bf X} is the matrix of the quantities that are {\it known
exactly}.  For example, we began with the example in which the elements
of {\bf X} were the times at which the measurements were taken. 

	Generally, then, the curvature and covariance matrix elements
depend on the {\it locations} of the datapoints (the ensemble of $t_m$
in equation \ref{one}) but not on the {\it measured values} (the
ensemble of $\theta_m$ in equation \ref{one}). And on your adopted
values for $\sigma_{meas,m}$. Because of this\dots

	\item Think {\it before} making your measurements about the
covariance matrix and how to minimize the off-diagonal elements. By
taking measurements at {\it well-chosen} times, or {\it 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. 

\end{itemize}

\section{BRUTE FORCE CHI SQUARED AND THE CURVATURE MATRIX}

\subsection{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 $\sigma^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 $\chi^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 \S \ref{easycalc}. 

\subsection{ 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 $\chi^2$ value for the unknown coefficient (i.e.,
the weighted average) is $x_0$ and that its variance is
$\sigma_{x_0}^2$.  Now consider offsets from $x_0$, namely $x_0 + \delta
x$.  Because $\chi^2$ is minimized at $x_0$, it is clear that

\begin{equation} 
\Delta \chi^2 = \chi^2_{x_0 + \delta x} - \chi^2_{x_0} = 
  {1 \over 2} {d^2\chi^2 \over d \delta x^2} \delta x^2
\end{equation}

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

\begin{equation}
  {1 \over 2} {d^2\chi^2 \over d \delta x^2} = 
   [\alpha_\sigma] = [\alpha_\sigma]_{00} = \sum \left( {1 \over
\sigma_{meas,m}^2} \right)
\end{equation}

	Because we have a single parameter $x_0$, its variance
$\sigma_{x0}^2$ is the offset $\delta x_0^2$  from $x_0$ that gives
$\Delta \chi^2 = 1$. If we had more than one variable, we'd define
ellipsoids in the manner discussed in \S \ref{easycalc}. 

	Are you bothered by the fact that the curvature is independent
of the data values? Then go back and read \S \ref{diatribe}.

\section{USING SINGULAR VALUE DECOMPOSITION (SVD) IN IDL} \label{SVD}

        You don't need to read this section to use the SVD technique in
IDL. Rather, you can just use IDL's documentation for {\bf 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 {\bf X} matrix we
have previously defined\footnote{It is a pleasure to thank Jason Wright
for telling me about this technique.}.  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 {\bf X} matrix and  
replace the values of the independent variable ($t_m$ in our examples)
by the measurement number $m$.

\subsubsection{ Step 0: Introduce a common block that will contain the 
{\bf X} matrix}

        Let's call this common block {\bf svdinput}. The relevant IDL
command is {\bf 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.

\subsubsection{ Step 1: Populate the {\bf X} matrix and {\bf Y} vector
almost as before}

        Here you populate the {\bf X} and {\bf Y} matrices just as
before in \S \ref{idlsolution} and \ref{idlexplicitone}, with some small
differences.  One, the measurement uncertainties are an input to {\bf
svdfit} so we don't create the weighted versions of {\bf X} and {\bf Y}.
 Two, we must make {\bf Y} a row vector instead of a column vector.
Here we repeat all the IDL commands as in \S \ref{idlexplicitone}, but
without the verbal commentary, and include this slight modification for
{\bf Y}.  We use the numbers in the numerical example of \S
\ref{idlsolution} to be absolutely unambiguous about what to do. In this
example, the number of unknowns $N=3$ and the number of measurements
$M=4$.
\begin{mathletters}
\begin{eqnarray}
{\bf X\_SVD} = \left[
\begin{array}{ccc}
1  &  5   &  25  \\
1  &  7   &  49  \\
1  &  9   &  81  \\
1  &  11  &  121 \\
\end{array} \; \right]
\end{eqnarray}

\begin{eqnarray}
{\bf Y} = \left[142, 168, 211, 251 \right] \; .
\end{eqnarray}

\noindent and we define {\bf ME}, the vector of measurement errors, to
be 

\begin{eqnarray}
{\bf ME} = \left[5., 3., 7., 4. \right] \; .
\end{eqnarray}
\end{mathletters}
        
\noindent where we have invented a four arbitrary values for
illustrative purposes.  You can, of course, use unity for all of the  
elements of {\bf ME}; also, you can call {\bf 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 \S \ref{idlexplicitone}.

\subsubsection{ Step 2: Call {\bf svdfit}}

        Now in in \S \ref{idlexplicitone} we replace \ref{eqn60} to
\ref{idlcovariancew} (and in \S \ref{idlsolution}, equation \ref{eqn24}
to \ref{lscoefferrorone}) by the following:
        
\begin{mathletters} \begin{equation} dataindex = indgen( M) 
\end{equation}

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

\begin{equation}
{\bf a}= {\bf svdfit(dataindex, Y, N, measure\_errors= ME,
  function\_name='svdfcn',
} $$
    $$  {\bf yfit=ybar, chisq={\it chi\_sq}, variance=vardc\_official, 
      sigma=sigdc\_official,} $$
$$ {\bf covar=XXWI, singular=singular) }
\end{equation}

\begin{equation}
{\bf ncov = XXWI/{\it sqrt}( vardc\_official \# vardc\_official)}
\end{equation}
\end{mathletters}

\noindent Finally, you have to define the following function. \begin{verbatim}
function svdfcn, dataindex, MIDL
common svdinput
return, x_svd[ *, dataindex]
end
\end{verbatim}
\noindent If you are doing least squares fitting and setting all
elements of {\bf ME} equal to unity, then of course {\bf 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 {\bf 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 \S 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.

\section{ NOTATION COMPARISON WITH NUMERICAL RECIPES}

	I learned least squares from Appendix A of Chauvenet (1863). He
didn't use $\chi^2$ and didn't use matrix techniques, but \S
\ref{sectionone} and \ref{sectiontwo} 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
\S \ref{chisqsection}. 

	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 $\sigma_{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 \S \ref{chisqsection}. 

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

{\begin{mathletters}
\begin {equation}
{\bf X} \longleftrightarrow {\bf A}
\end{equation}

\begin{equation}
{\bf Y} \longleftrightarrow {\bf b}
\end{equation}

\begin{equation}
{\bf X^T \cdot X} = {\bf XX} = [\alpha] \longleftrightarrow 
   {\bf A^T \cdot A} = [\alpha]
\end{equation}

\begin{equation}
{\bf XX}^{-1} = {\bf XXI} = [\alpha]^{-1} \longleftrightarrow 
   [\alpha]^{-1} = [C] = {\bf C}  
\end{equation}

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

\begin{equation}
 N \longleftrightarrow M
\end{equation}

\begin{equation}
 M \longleftrightarrow N
\end{equation}

\end{mathletters}

\section{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 $\theta$, as in equation~\ref{one}.  This
essential assumption means that $t$ is known with high precision and all
the uncertainty is in $\theta$, and you are minimizing the squares of
the residuals in the $\theta$-direction only. 

	If {\it 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:

	{\bf (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=(\theta + t)$ and $y=(\theta - t)$.  The sum, $x$, has smaller
fractional uncertainty than $y$ so its dispersion {\it 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!

	{\bf (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!

\begin{references}

\reference {} Bevington, P.R. \& Robinson, D. 1992, Data Reduction and
Error Analysis for the Physica Sciences (WCB/McGraw-Hill).

\reference {} Chauvenet, W. 1863, A Manual of Spherical and Practical
Astronomy, Dover Press.

\reference {} Cowan, G. 1998, Statistical Data Analysis, Clarendeon
Press.

\reference {} Press, W.H., Flannery, B.P., Teukolsky, S.A., \&
Vetterling, W.T. 2001, Numerical Recipes (second edition), Cambridge
University Press.

\reference {} Taylor, J.R. 1997, An Introduction to Error Analysis,
University Science Books.

\end{references}
\end{document}


\subsection{The uncertainties of the derived coefficients:  discussion
ignoring covariance}

\label{chicoeffsnoco}

	In chi-squared fitting, you know beforehand the intrinsic
uncertainty in each measurement $\sigma_{meas,m}$. This is a function
{\it only} of the {\it instrument} and the {\it data-taking procedure}:
if you know the properties of the instrument, and if you understand how
it works, and if there are no systematic errors, then you should be able
to specify each $\sigma_{meas,m}$. Specifically, $\sigma_{meas,m}$ is
most certainly {\it not} the same thing as the dispersion of the points
from the fit. 



        If you read any standard textbook on least squares fitting, you
will discover that the error in a derived parameter is equal to its
diagonal element in the covariance matrix. That is, the formula is just
like our equation \ref{coeffvarianceone} but with two differences: the
$\sigma$ subscripts disappear and there is no $s^2$. See NM equation
14.3.15; BR 7.25, Taylor 8.16, 8.17. This is not as bad as it sounds
because, without the $\sigma$ subscripts, the covariance matrix
$[\alpha]^{-1}$ contains the factor $\sigma_{meas}^2$. If you've chosen
$\sigma_{meas}^2$ accurately, and if the model for the fit is good so
that the fitted curve accurately represents the points, then $s^2
\approx \sigma_{meas}^2$ so the results are equivalent.
        
        However, realize that in Chi-squared fitting the covariance
matrix $[\alpha]^{-1}$ depends only on $\bf X$, i.e.\ on the times at  
which you took the data---{\it not on the datapoint values}. Because
$[\alpha]^{-1}$ is {\it independent of the datapoint values}, 
your derived parameter errors are also independent of the datapoint   
values. So they don't depend on the datapoint variance  $\sigma^2$ or
the $\chi^2$.

        This means that the standardly-derived errors in the parameters
depend only on what you chose for $\sigma_{meas}$. If you have an   
inflated opinion of the quality of your instrument or observing prowess
and choose an unrealistically low value for $\sigma_{meas}$, then your
derived parameter errors are too low by the same factor. Or if your
estimate of $\sigma_{meas}$ is accurate but you are fitting a poor
model, i.e.\ fitting a quadratic when you should be fitting a sine wave,
then the residuals will be large and your parameter errors will also be
much larger than you calculate.
   
        How many observational papers in the ApJ provide parameter
errors based on Chi-squared fitting that are unrealistically low? Or,
for that matter, high?  [probably fewer!]

	To understand the derived-parameter uncertainty issue, we first
illustrate the situation with the simplest of weighted chi square fits,
the weighted mean. We also treat the case in which all assumed variances
of the datapoints are identical ($\sigma_{meas,m} = \sigma_{meas}$) so
that we can focus on the essentials of the problem; this restriction in
no way affects our derived results. This will provide insight that we
can extend to the multivariate chi-squared case.



