5.4 Fitting an arbitrary function ...  5. Least-Squares Methods ...  5.2 Fitting a straight line ...

5.3 Fitting functions linear in parameters

Consider the case where the function $f({\bf x})$ consists of a linear combination of other functions, $g({\bf x})$:

$$f({\bf x}) = \sum_k a_k g_k({\bf x}) .$$ (5.27)

 

The functions g_k may be polynomials, such as $\{1,x,x^2$...$\}$ or in some applications the Legendre polynomials, but the procedure developed here applies to any functions. The functions are not necessarily linear in x; the following development only requires that the relationship between the measurements and the general functions $g_k({\bf x})$ be linear as expressed in (5.27). The chisquare function for this case, if deviations from the measurements $\{y_i\}$ corresponding to the values of the parameter $\{{\bf x}_i\}$ are characterized by a Gaussian error distribution, is

$$\chi^2 = \sum_i {{\bigl(y_i-\sum_ka_kg_k({\bf x}_i)\bigr)^2}\over{\sigma_i^2}}$$ (5.28)

 

where $\sigma_i$ is the standard deviation of the measurement error. The minimum value of $\chi^2$ satisfies

$${{\partial\chi^2}\over{\partial a_\ell}} = \sum_i {{2\bigl(y_......g_k({\bf x}_i)\bigr)g_\ell({\bf x}_i)}\over{\sigma_i^2}} = 0 .$$ (5.29)

 

From (5.3), the information matrix for this case is

$${{1}\over{2}} {{\partial^2\chi^2}\over{\partial a_\ell\partia...... \sum_i{{g_\ell({\bf x}_i)g_m({\bf x}_i)}\over{\sigma_i^2}}.$$ (5.30)

 

If we define

$$M_\ell=\sum_i{{y_ig_\ell({\bf x}_i)}\over{\sigma_i^2}}$$ (5.31)

 

then

$${\bf M} = {\bf H a}$$ (5.32)

 

or

$${\bf a} = {\bf H}^{-1} {\bf M} .$$ (5.33)

 

The fit procedure is then:

1.

Find the information matrix H from (5.30), and the column matrix M from (5.31).

2.

Invert the matrix H to obtain the error matrix, H^-1.

3.

Multiply H^-1 and M, according to (5.33), to obtain the coefficients a that represent the least-squares fit.

Alternately, there are procedures that obtain the solution a directly from (5.32), perhaps while also calculating the inverse of H. See, for example, the Gauss-Jordan inversion method described by Press et al. (1992), p. 36, and other methods for the solution of (5.32) described there.

This procedure will fail if the matrix H is not invertible. This will occur if the functions g_k are not linearly independent, because then the determinant of the matrix H will be zero and it cannot be inverted. Because fits are often obtained using computer routines, rounding errors may cause the calculated determinant to differ from zero even in cases that are not invertible, so it is usually useful to check the value of the determinant and flag cases where that value is near the minimum resolvable value.
 

Exercise 5.1: Show that the chisquare for the fit is given by

$$\chi^2 = \sum_i {{y_i^2}\over{\sigma_i^2}} - 2 {\bf a}^t{\bf M}+ {\bf a}^t{\bf H}{\bf a} \ .$$ (5.34)

 

This is a useful relationship because, if the summation in the first term on the right side is accumulated while the matrices H and M are calculated for the fit, the chisquare for the result can be determined without another pass through the data. Beware of numerical precision problems when using this in computer routines, because the chisquare may be a small difference between very large numbers.

   
 5.4 Fitting an arbitrary function ...  5. Least-Squares Methods ...  5.2 Fitting a straight line ...