Introduction to Data Analysis

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2006

Initialization

ClearAll["Global`*"];
Off[General::spell, General::spell1]

$DefaultFont = {"Times", 12};
$TextStyle = {FontFamily -> "Times", FontSize -> 12, 
          FontSlant -> "Italic"};

Needs["Graphics`"];
Needs["Utilities`Notation`"];
Needs["Statistics`"];

Symbolize[(Subscript[σ, x])];Symbolize[(Subscript[σ, y])];
Symbolize[(Subscript[S, x])];Symbolize[(Subscript[S, y])];Symbolize[(Subscript[S, x, x])];Symbolize[(Subscript[S, x, y])];

Symbolize[\(y\&_\)]

It is also useful to append the directory containing this notebook and related files to $Path.

AppendTo[$Path, "C:\Documents and Settings\Jim\My Documents\My \
eBooks\EssentialMathematica\Data analysis"];

You will need to adjust this directory name for your computer.

Principle of maximum likelihood

Let x represent an independent variable, such as an applied voltage that we can vary freely and measure accurately, and let y represent the dependent variable, such as the current in a circuit.  We often think of the dependent variable as the response of a system to a variable that is controlled externally.  We usually have less control over the response and measure the dependent variable less accurately than the independent variable; therefore, we assume that the random errors, (Subscript[σ, x]), are negligible for the independent variable and consider only random errors, (Subscript[σ, y]), in the dependent variable.  Our data set then consists of the set of N measurements {{(Subscript[x, 1]), (Subscript[y, 1]), (Subscript[σ, 1])}, {(Subscript[x, 2]), 
          (Subscript[y, 2]), (Subscript[σ, 2])}, ⋯, {(Subscript[x, N]), 
          (Subscript[y, N]), (Subscript[σ, N])}}.

We imagine that our data set is merely one outcome out of an infinite number of possible outcomes and that there exists a parent probability function which determines the likelihood of our particular data set.  Furthermore, we assume that each measurement is independent such that the joint probability for the entire data set

P[{(Subscript[y, i])}]  =  (∏_(i = 1)^N) (Subscript[P, i])[
          (Subscript[y, i])]

factors into a product of probabilities for each individual measurement.  For counting experiments we can be pretty sure that (Subscript[P, i]) should be a Poisson distribution, but usually we do not know the proper (Subscript[P, i]) functions very accurately.  In the absence of additional information, we usually assume that (Subscript[P, i]) is a normal or Gaussian distribution of the form

\(TraditionalForm\`P\_i[y\_i]\  ∝ \ \(σ\_i\^\(-1\)\) 
          Exp[\(-\((\(y\_i - y[x\_i; a]\)\/σ\_i)\)\^2\)]\ \ \ ⟹\ \ \ Log[
              P]\  = \ \(-\(∑\+\(i = 1\)\%N\((\(y\_i - y[x\_i; a]\)\/σ\_i)\)\^2\)\)\  + \ constant\)

where y[x; a] is a fitting function that depends upon the independent variable x and a set of m parameters a = {(Subscript[a, μ]), μ = 1, m} that represent the properties of the system that we seek to determine from our data set.  Such a fitting function is based upon a physical model of the system.  Unless one is paranoid, it is reasonable to assume that our data set is close to the most likely result for an experiment of this type.  Therefore, the principle of maximum likelihood claims that the best estimates of the model parameters are those that maximize the joint probability for our data set.  Because Log is a monotonic function, maximizing Log[P] is equivalent to maximizing P, and is much simpler.  Thus, one must solve the system of equations

\(TraditionalForm\`∂Log[P]\/∂a\_μ == 
      0\ \ ⟹\ \ \ \(∑\+\(i = 1\)\%N\(∂y[
                    x\_i; a]\/∂a\_μ\) \(y\_i - y[x\_i; a]\)\/σ\_i\^2\)\  == \ 0\)

for the parameters that maximize the likelihood.  These solutions are denoted \(TraditionalForm\`a\&^ = {a\&^\_μ}\) and are described as maximum-likelihood estimators of the model parameters.  To minimize confusion, it is useful to employ latin indices (i, j) to enumerate data points and greek indices (μ, ν) to enumerate model parameters.

The method of maximum likelihood is quite general and can be used for a wide variety of sample probability distributions (Subscript[P, i]).  When these are assumed to be Gaussian, the most common situation, maximizing Log[P] becomes equivalent to minimizing (χ^2)[a] where

\(TraditionalForm\`\(χ\^2\)[
        a]\  = \ ∑\+\(i = 1\)\%N\((\(y\_i - y[x\_i; a]\)\/σ\_i)\)\^2\)

is a weighted sum of the squares of the difference between observed and fitted quantities.  Thus, for normal distributions we use the method of least-squares.  Although not strictly correct, one often uses least squares to fit data based upon the Poisson distribution also.

It is important to recognize that the values we obtain for \(TraditionalForm\`{a\&^\_μ}\) are merely estimates of the physical parameters based upon a particular finite set of measurements.  We expect that these estimates will become increasingly accurate as the number and/or precision of the measurements increases, but we can never be absolutely certain that measured values are exactly equal to true values.  Therefore, we describe the fitted parameters as maximum-likelihood estimators and often adorn the symbols \(TraditionalForm\`{a\&^\_μ}\) with carets to distinguish estimators from true parameters of the physical system.  The estimators are functions of measured quantities that are subject to random errors and, consequently, are themselves subject to random errors.  If we repeat the experiment many times, with each repetition containing the same sequence of measurements, we will obtain a distribution of results for each of the estimators.  Nevertheless, estimates of the accuracy of the parameters can be obtained from the curvature of (χ^2)[a] near its minimum.  If the minimum is deep and the curvature with respect to (Subscript[a, μ]) is large, then parameter (Subscript[a, μ]) is determined well by the data but if the curvature is small the uncertainty in that parameter will be large.

The quality of the fit is gauged by the reduced chi-square

\(TraditionalForm\`χ\_red\%2\  = \ \(\(1\/\(N - m\)\) \(χ\^2\)[
            a\&^]\  = \ \(1\/\(N - 
                m\)\) \(∑\+\(i = 1\)\%N\((\(y\_i - y[x\_i; a\&^]\)\/σ\_i)\)\^2\)\)\)

where m is the number of parameters and N - m is the number of degrees of freedom.  If the model is reliable and our estimates of the measurement uncertainties are accurate, we expect to find (Subscript[χ, red]^2) ≈ 1 for a good fit.  However, this statistic is a random variable with its own distribution — the larger the data set the narrower this distribution becomes.

Maximum likelihood methods are quite versatile and can be applied to problems in which each measurement involves several dependent variables and produces several independent variables.  The probability distributions for each measurement can also be much more sophisticated.  It may be useful to apply constraints upon the fitted parameters.  Naturally, the most efficient methods for solving these maximum-likelihood problems depend upon the nature of the data set and the fitting function.  In the next few sections we consider increasingly sophisticated methods, starting with the simplest case of a linear function.

Linear least squares

Two-parameter model, linear wrt x

Suppose that our physical model of the system suggests that y[x] = (Subscript[a, 0]) + (Subscript[a, 1]) x is a linear function of two parameters {(Subscript[a, 0]), (Subscript[a, 1])}.  Note that in interests of streamlining the notation, it is customary to leave the dependence of the fitting function upon the model parameters implicit.  The most important feature of this model is that the fitting function is linear with respect to the properties of the system that we seek to measure, namely the parameters {(Subscript[a, 0]), (Subscript[a, 1])}.  The fact that it is also linear with respect to x simplifies some of the algebra in this section, but is not an essential feature of the method of linear least squares.

Maximizing the probability provides the system of normal equations

\(TraditionalForm\`∂χ\^2\/∂a\_μ == 
      0\ \ ⟹\ \ \ \(∑\+\(i = 1\)\%N\(∂y[
                    x\_i]\/∂a\_μ\) \(y\_i - y[x\_i]\)\/σ\_i\^2\)\  == \ 0\)

where μ ∈ {0, 1} for our two-parameter model.  Thus,

\(TraditionalForm\`\(\(\ \ \)\(∑\+\(i = 1\)\%N\( y\_i - a\_0 - \(a\_1\) x\_i\)\/σ\_i\^2\  == \ 
      0\ \ ⟹\ \ \ S\_y - \(a\_0\) W - \(a\_1\) 
          S\_x\  == \ 0

    \ \ ∑\+\(i = 1\)\%N\( 
              x\_i\) \(y\_i - a\_0 - \(a\_1\) x\_i\)\/σ\_i\^2\  == \ 
      0\ \ ⟹\ \ \ S\_\(x, y\) - \(a\_0\) 
          S\_x - \(a\_1\) S\_\(x, x\)\  == \ 0\)\)\)

where

(Subscript[S, x])  =  (∑_(i = 1)^N) (Subscript[w, i]) 
          (Subscript[x, i])
(Subscript[S, y])  =  (∑_(i = 1)^N) (Subscript[w, i]) 
          (Subscript[y, i])
(Subscript[S, x, x])  =  (∑_(i = 1)^N) (Subscript[w, i]) 
          ((Subscript[x, i])^2)
(Subscript[S, x, y])  =  (∑_(i = 1)^N) 
              (Subscript[w, i]) (Subscript[x, i]) (Subscript[y, i])
W = (∑_(i = 1)^N) (Subscript[w, i])

are expressed in terms of the weights (Subscript[w, i]) = 1/((Subscript[σ, i])^2).  Therefore, we obtain the maximum-likelihood estimators

\(TraditionalForm\`a\&^\_0\  = \ \(\(S\_y\) S\_\(x, x\) - \(S\_x\) S\_\(x, y\)\)\/Δ\ , \ \ \ \ a\&^\_1\  = \ \(W\ S\_\(x, y\) - \(S\_x\) S\_y\)\/Δ\)

where

Δ  =  
      W (Subscript[S, x, x])  -  (Subscript[S, x]^2)

is the discriminant for this pair of equations.  These are simple formulas that can easily be evaluated using a spreadsheet, Mathematica, or other software.  They can even be evaluated by a hand calculator or (gasp!) using pencil and paper if N is not large.

The uncertainty in a derived quantity, such as \(TraditionalForm\`a\&^\_0\), can be estimated using our familiar techniques for error propagation.  We assume that the measurements are independent and write

\(TraditionalForm\`σ[a\&^\_0]\^2\  = \ ∑\+\(i = 1\)\%N\((\(∂a\&^\_0\/∂y\_i\) σ\_i)\)\^2\ , \ \ \ \ σ[a\&^\_1]\^2\  = \ ∑\+\(i = 1\)\%N\((\(∂a\&^\_1\/∂y\_i\) σ\_i)\)\^2\)

where the contribution of the uncertainty in measurement (Subscript[y, i]) is determined by the variation of \(TraditionalForm\`a\&^\_μ\) with respect to (Subscript[y, i]).  Thus, substituting

\(TraditionalForm\`∂a\&^\_0\/∂y\_i\  = \ \(\(Δ\^\(-1\)\)(\(S\_\(x, 
                  x\)\) ∂S\_y\/∂y\_i - \(S\_x\) ∂S\_\(x, y\)\/∂y\_i)\ \  = \ \(w\_i\/Δ\) \((S\_\(x, x\) - \(S\_x\) x\_i)\)\)\)
\(TraditionalForm\`∂a\&^\_1\/∂y\_i\  = \ \(\(Δ\^\(-1\)\)(
          W\ ∂S\_\(x, y\)\/∂y\_i - \(S\_x\) ∂S\_y\/∂y\_i)\ \  = \ \(w\_i\/Δ\) \((W\ x\_i - S\_x)\)\)\)

into

\(TraditionalForm\`σ[a\&^\_0]\^2\  = \ \(Δ\^\(-2\)\ \ \(∑\+\(i = 1\)\%N\(\( w\_i\)(S\_\(x, x\) - \(S\_x\) x\_i)\)\^2\)\  = \ \(Δ\^\(-2\)\ \ \(∑\+\(i = 1\)\%N\( w\_i\)(
                S\_\(x, x\)\^2 - 
                  2 \( S\_x\) \(S\_\(x, x\)\) x\_i + \(S\_x\^2\) 
                    x\_i\^2)\)\  = \ \(\(S\_\(x, x\)\)(W\ S\_\(x, x\) - S\_x\^2)\)\/Δ\^2\)\)\)

\(TraditionalForm\`σ[a\&^\_1]\^2\  = \ \(Δ\^\(-2\)\ \ \(∑\+\(i = 1\)\%N\(\( w\_i\)(W\ x\_i - S\_x)\)\^2\)\  = \ \(Δ\^\(-2\)\ \ \(∑\+\(i = 1\)\%N\( w\_i\)(\(W\^2\) x\_i\^2 - 
                  2  W\ \(S\_x\) x\_i + 
                  S\_x\^2)\)\  = \ \(W(W\ S\_\(x, x\) - S\_x\^2)\)\/Δ\^2\)\)\)

we obtain

\(TraditionalForm\`σ[a\&^\_0]\^2\  = \ 
      S\_\(x, x\)\/Δ\ , \ \ \ \ \ \ σ[a\&^\_1]\^2\  = \ 
      W\/Δ\)

relatively easily.

Often we do not know (Subscript[σ, i]) but believe that the uncertainties are pretty uniform.  Because the weights appear linearly in both numerator and denominator, only relative weights affect the fitted values for the parameter; hence, we can set all (Subscript[w, i]) -> 1 if the weights are uniform.  The measurement uncertainties can then be estimated using the standard deviation of the data from the model, according to

\(TraditionalForm\`σ\_y\^2 = \ \(1\/\(N - 
              2\)\) \(∑\+\(i = 1\)\%N\((y\_i - y\&_\_i)\)\^2\)\)

where the expected or fit values

\(TraditionalForm\`y\&_\_i\  = \ y[x\_i; a\&^]\)

are based upon our best estimates of the model parameters and the number of degrees of freedom is reduced by the number of parameters fit to the data, here  two.  Under these circumstances the uncertainties in fitted parameters become

\(TraditionalForm\`w\_i -> 
        1\ \ \ \ ⟹\ \ \ σ[
            a\&^\_0]\  = \ \(σ\_y\) \@\(S\_\(x, x\)\/Δ\)\ , \ \ \ \ \ \ σ[
        a\&^\_1]\  = \ \(σ\_y\) \@\(N\/Δ\)\)

The quality of the fit should be judged graphically.  If the data scatter fairly symmetrically about the fitted line, the model is probably reliable and (Subscript[σ, y]) provides an estimate of the quality of the data.  However, if the scatter is not uniform the assumption of a linear model becomes dubious.

example: unweighted linear fit

To illustrate this procedure, we construct a data set with uniformly spaced values of the independent variable and values of the dependent variable based upon a linear function with random fluctuations drawn from a normal distribution with prescribed σ.  Pretending that we do not know σ, the weights are initialized to unity.

x = Range[-0.9, 0.9, 0.1];
y = 
        117  -  23  x + 
          Table[Random[
              NormalDistribution[0, 2. ]], {Length[
                x]}];
w = Table[1, {Length[x]}];

Next we evaluate the sums defined above using the listability of Mathematica functions

(Subscript[S, x]) = w . x;
(Subscript[S, y]) = w . y;
(Subscript[S, x, x]) = w . (x^2);
(Subscript[S, x, y]) = w . (x*y);
W = Plus @@ w;
Δ = W (Subscript[S, x, x]) - ((Subscript[S, x])^2);

and evaluate the estimators and their uncertainties.

{(Subscript[a, 0]), 
        (Subscript[a, 1])} = {(((Subscript[S, y]) (Subscript[S, x, x]) - (Subscript[S, x]) (Subscript[S, x, y]))/Δ),  ((W (Subscript[S, x, y]) - (Subscript[S, x]) (Subscript[S, y]))/Δ)} // 
        N

\(\(y\&_ = a\_0 + \(a\_1\) x;\)\)

\(σ\_y\  = \ 
      Sqrt[Plus @@ \(\((y - y\&_)\)\^2\)/\((Length[y] - 2)\)]\)

{(Subscript[σ, 0]), (Subscript[σ, 1])} = (Subscript[σ, y]) 
          Sqrt /@ {((Subscript[S, x, x])/Δ),  ((W )/Δ)} // N

Finally, we compare the fit with the data.

DisplayTogether[
ListPlot[
          Evaluate[Transpose[{x, y}]], 
          PlotStyle -> PointSize[0.02]], 
Plot[
          (Subscript[a, 0]) + (Subscript[a, 1]) x, {x, -1, 1}], Frame -> True, 
        FrameLabel -> {"x", "y"}];

The data do appear to scatter pretty symmetrically about the fit, but a clearer appraisal is provided by plotting the residuals

\(TraditionalForm\`r\_i = \(y\_i - y\&_\_i\)\/σ\_i\)

If the model is reliable, the residuals should be scattered uniformly and if the data are normally distributed approximately two-thirds should be contained within ±1.  The present data conform with these expectations by construction.

\(\(DisplayTogether[
        ListPlot[Evaluate[Transpose[{x, \(y - y\&_\)\/σ\_y}]], 
          PlotRange -> {Automatic, {\(-3\), 3}}, 
          PlotStyle -> PointSize[0.02], Frame -> True, 
          FrameLabel -> {"x", "residuals"}], 
        Plot[{\(-1\), 1}, {x, \(-1\), 1}, 
          PlotStyle -> Dashing[{0.02, 0.02}]]];\)\)

example: linear fit to nonlinear data

Next suppose we construct a data set for a nonlinear model

x = Range[-0.9, 0.9, 0.1];
y = 
        117  -  23  x + 15  (x^2) + 
          Table[Random[
              NormalDistribution[0, 2. ]], {Length[
                x]}];
w = Table[1, {Length[x]}];

but attempt to fit those data using linear least squares

(Subscript[S, x]) = w . x;
(Subscript[S, y]) = w . y;
(Subscript[S, x, x]) = w . (x^2);
(Subscript[S, x, y]) = w . (x*y);
W = Plus @@ w;
Δ = W (Subscript[S, x, x]) - ((Subscript[S, x])^2);

{(Subscript[a, 0]), 
        (Subscript[a, 1])} = {(((Subscript[S, y]) (Subscript[S, x, x]) - (Subscript[S, x]) (Subscript[S, x, y]))/Δ),  ((W (Subscript[S, x, y]) - (Subscript[S, x]) (Subscript[S, y]))/Δ)} // 
        N

\(\(y\&_ = a\_0 + \(a\_1\) x;\)\)

\(σ\_y\  = \ 
      Sqrt[Plus @@ \(\((y - y\&_)\)\^2\)/\((Length[y] - 2)\)]\)

{(Subscript[σ, 0]), (Subscript[σ, 1])} = (Subscript[σ, y]) 
          Sqrt /@ {((Subscript[S, x, x])/Δ),  ((W )/Δ)} // N

The standard deviation is now substantially larger and comparisons between the fit and the data show pronounced systematic deviations, not just random fluctuations.  Although we cannot reject the linear model with absolute certainty without reliable estimates for the measurement uncertainties, these figures cast considerable doubt upon the hypothesis that y is a linear function of x and suggest that the model should be expanded to include an (x^2) contribution (note that parabolic nature of the residuals).

DisplayTogether[
ListPlot[
          Evaluate[Transpose[{x, y}]], 
          PlotStyle -> PointSize[0.02]], 
Plot[
          (Subscript[a, 0]) + (Subscript[a, 1]) x, {x, -1, 1}], Frame -> True, 
        FrameLabel -> {"x", "y"}];

\(\(DisplayTogether[
        ListPlot[Evaluate[Transpose[{x, \(y - y\&_\)\/σ\_y}]], 
          PlotRange -> {Automatic, {\(-3\), 3}}, 
          PlotStyle -> PointSize[0.02], Frame -> True, 
          FrameLabel -> {"x", "residuals"}], 
        Plot[{\(-1\), 1}, {x, \(-1\), 1}, 
          PlotStyle -> Dashing[{0.02, 0.02}]]];\)\)

Basis functions

Next we consider a more general fitting function of the form

y[
        x]  =  (∑_(μ = 1)^m) (Subscript[a, μ]) 
          (Subscript[f, μ])[x]

where a = {(Subscript[a, μ])} is the set of parameters we seek and f = {(Subscript[f, μ])} is a set of basis functions that depend upon the independent variable x but not upon the parameters (Subscript[a, μ]).  Therefore, y[x] is a linear function of m parameters.  The normal equations

\(TraditionalForm\`∂χ\^2\/∂a\_μ == 
      0\ \ ⟹\ \ \ \(∑\+\(i = 1\)\%N\(∂y[
                    x\_i]\/∂a\_μ\) \(y\_i - y[x\_i]\)\/σ\_i\^2\)\  == \ 
      0\ \ ⟹\ \ \(∑\+\(ν = 1\)\%m\( 
                a\_ν\) \(∑\+\(i = 1\)\%N\( w\_i\) f\_μ[x\_i] 
                f\_ν[x\_i]\)\)\ \  == \ ∑\+\(i = 1\)\%N\( 
              w\_i\) f\_μ[x\_i] y\_i\)

can be expressed in the matrix form

M a  ==  b

where

(Subscript[M, μ, ν])  =  (∑_(i = 1)^N) 
              (Subscript[w, i]) (Subscript[f, μ])[(Subscript[x, i])] (Subscript[f, ν])[(Subscript[x, i])]

is the curvature matrix for the experiment, a = {(Subscript[a, μ]), μ = 1, m} is the parameter vector, and

(Subscript[b, μ])  =  (∑_(i = 1)^N) (Subscript[w, i]) (Subscript[y, i]) 
          (Subscript[f, μ])[(Subscript[x, i])]

is the gradient vector.  The terminology for β and M arises from their origin in first and second derivatives of (χ^2).  Notice that the curvature matrix depends upon the values of the independent variable and the precision of the measurements, but not upon the actual results of those measurements.  It represents the sensitivity of the proposed measurements to the parameters of the system.  The gradient vector projects the actual measurements onto the basis functions.  Thus, the maximum-likelihood estimates of the model parameters are given by

\(TraditionalForm\`\(a\&^\) = \(M\^\(-1\)\) b\)

and the problem is reduced to inversion of a symmetric matrix.  One can easily show that the previous results are recovered for a two-parameter linear model.  (Try it!)

The error analysis applies the same technique used for the two-dimensional problem.  It is convenient to identify the error matrix ε = (M^(-1)) with the inverse of the curvature matrix; the reason for this terminology will become apparent momentarily.  Starting with

\(TraditionalForm\`σ[a\&^\_μ]\^2\  = \ ∑\+\(i = 1\)\%N\((\(∂a\&^\_μ\/∂y\_i\) σ\_i)\)\^2\)

and using

\(TraditionalForm\`∂a\&^\_μ\/∂y\_i\  = \ \(∑\+\(ν = 1\)\%m\( ε\_\(μ, ν\)\) ∂b\_ν\/∂y\_i\  = \ \ ∑\+\(ν = 1\)\%m\( ε\_\(μ, ν\)\) \(w\_i\) 
            f\_ν[x\_i]\)\)

we soon obtain

\(TraditionalForm\`σ[a\&^\_μ]\^2\  = \ \(∑\+i\( 
                w\_i\) \(∑\+\(ν, γ\)\((\(ε\_\(μ, ν\)\) f\_ν[x\_i] \(ε\_\(μ, γ\)\) 
                  f\_γ[
                    x\_i])\)\)\  = \ \(∑\+\(ν, γ\)\(ε\_\(μ, ν\)\) \(ε\_\(μ, γ\)\) \(∑\+i\( w\_i\) f\_ν[x\_i] 
                  f\_γ[
                    x\_i]\)\  = \(∑\+\(ν, γ\)\(ε\_\(μ, ν\)\) \(ε\_\(μ, γ\)\) 
                M\_\(γ, ν\)\ \  = \ ε\_\(μ, μ\)\)\)\)\)

where in the last step contraction of γ index yields the unit matrix (Subscript[δ, μ, ν]) and contraction of the ν index isolates the appropriate diagonal element.  Therefore, diagonal elements of the error matrix give the variances in the corresponding parameters.

linfit::usage="linfit[data,f,var] performs linear least-squares fitting.  The data array should be a list of either {x,y,dy} or {x,y} for each measurement, where the fit is unweighted if uncertainties in the dependent variable are omitted.  The list of basis functions is f and var is the symbol used in f for the independent variable.";
linfit::BadDim = "Bad dimensions in data array: `1`";

linfit[data_, 
        f_?VectorQ, var_, opts___Rule] := 
      Module[{x, y, dy, w, M, b, ε, basis, ndata, nparam, ndim, 
          i, j, params, yfit, redchi, errors, 
          stdev}, 

 (* Extract x, y, 
          dy vectors from data array. *) 

{ndata, ndim} = Dimensions[data]; 
If[
          ndata ≤ Length[f] || ndim < 2 || 
            ndim > 3, (Message[linfit::BadDim, {ndata, ndim}]; 
            Return[])]; 
If[
          ndim == 2, 
({x, y} = Transpose[data]; 
            w = Table[1, {ndata}]), 
({x, y, dy} = 
              Transpose[data]; 
            w = 1/(dy^2))]; 

 (* 
          Accumulate the gradient vector and the curvature matrix. *) 

nparam = 
          Length[f]; 
M = 
          Table[0, {i, nparam}, {j, nparam}]; 
b = 
          Table[0, {i, 
              nparam}]; 
Do[
basis = 
            f /. var -> 
                x[[i]]; 
M = 
            M + w[[i]] 
                Outer[Times, basis, basis]; 
b = 
            b + w[[i]] 
                y[[i]] 
                basis, 
{i, 1, 
            ndata}]; 

 (* 
          Evaluate fit parameters,  error matrix,  
          and reduced (χ^2). *) 

ε = 
          PseudoInverse[M]; 
params = ε . 
            b; 
errors = 
          Table[(Sqrt[ε])[[i, i]], {i, 
              nparam}]; 
yfit = (params . f) /. 
            var -> x; 
redchi = 
          w . (((y - yfit))^2)/(ndata - 
                nparam); 
stdev = (Sqrt[redchi]); 

 (* 
          Return results as replacement rules. *) 

If[
          ndim == 3, {Parameters -> params, 
            Uncertainties -> errors, 
            ReducedChiSquare -> redchi}, {Parameters -> params, 
            Uncertainties -> errors*stdev, 
            StandardDeviation -> stdev}]]

?linfit

example

We return to the same data set used in the preceding example based upon a parabolic model.  If we fit these data using a quadratic, we obtain reduced chi-square near unity and recover the model parameters with deviations similar to their estimated uncertainties.  The points scatter normally around the fit with residuals that show no systematic dependence upon x.  Life is good!

f = {1, z, (z^2)};
dy = Table[2. , {Length[y]}];
data = Transpose[{x, y, dy}];
results = linfit[data, f, z]

\(\(y\&_ = \(Parameters . f /. results\) /. z -> x;\)\)

DisplayTogether[
ErrorListPlot[
          data], 
Plot[
          Evaluate[Parameters . f /. results], {z, -1, 1}], 
        Frame -> True, 
        FrameLabel -> {"x", "y"}];

\(\(DisplayTogether[
        ListPlot[Evaluate[Transpose[{x, \(y - y\&_\)\/dy}]], 
          PlotRange -> {Automatic, {\(-3\), 3}}, 
          PlotStyle -> PointSize[0.02], Frame -> True, 
          FrameLabel -> {"x", "residuals"}], 
        Plot[{\(-1\), 1}, {x, \(-1\), 1}, 
          PlotStyle -> Dashing[{0.02, 0.02}]]];\)\)

What happens if we try to fit too many parameters?  The uncertainties in the parameters begin to grow and some of the parameters become unstable without improving the fit.  Clearly, it does not pay to include too many parameters.

linfit[data, {1, z, (z^2), (z^3)}, z]

linfit[data, {1, z, (z^2), (z^3), (z^4)}, z]

Linearizing a model

Often a model that does not appear to be a linear function of its parameters can be transformed to linear form, but one must be careful to transform the uncertainties in a consistent manner.  For example, consider a radioactive decay model

y[t]  =  A Exp[-t/τ]

with parameters {A, τ}, where τ is the lifetime we seek to measure and y is the number of counts in an interval around time t.  To apply the method of linear least squares, we define

z  =  Log[y]  =  
        A  -  t/τ

with uncertainties

\(TraditionalForm\`σ\_z\  = \ \(\(∂z\/∂y\) σ\_y\  = \ σ\_y\/y\)\)

The uncertainties for a counting experiment are generally governed by the Poisson distribution, such that

\(TraditionalForm\`σ\_y = \(\(\@y\)\ \ \ ⟹\ \ σ\_z\  = \ \(\@y\/y = 1\/σ\_y\)\)\)

To illustrate this technique, we generate a set of data using random number based upon the Poisson distribution.  Although (χ^2) is based upon the Gaussian distribution and is not strictly correct for data drawn from the Poisson distribution, it is common to ignore this subtlety and to fit counting data by minimizing (χ^2).  Perhaps we will return to that issue later.  Also note that we apply a slight offset y in order to avoid numerical errors due to Log[0] or 1/0, but that the results are insensitive to the value of the offset if it is small enough.

t = Range[0, 1, 0.02];
p = {(10^3), 0.1};
yfit[t_] := 
        p[[1]] 
          Exp[-t/
              p[[2]]];
y = (10^(-8)) + Random[PoissonDistribution[yfit[#]]] & /@ 
          t;
dy = (Sqrt[y]);  z = Log[y]; dz = Abs[dy/y];
zresults = linfit[Transpose[{t, z, dz}], {1, q}, q]

The parameters returned by linfit represent Log[A] and -1/τ because those are the coefficients of the basis functions.  The desired parameters and their uncertainties are obtained according to

\(TraditionalForm\`A = 
      Exp[a\_0]\ , \ σ\_A\^2\  = \ \(\((\(∂A\/∂a\_0\) σ\_0)\)\^2\ \ ⟹\ \ \ σ\_A\  = \ 
        A\ σ\_\(\_0\)\)\)
\(TraditionalForm\`τ = \(-\(1\/a\_1\)\), \ \ σ\_τ\^2\  = \ \(\((\(∂τ\/∂a\_1\) σ\_1)\)\^2\ \ ⟹\ \ σ\_τ\  = \ σ\_1\/a\_1\^2\)\)

where the parameters are labeled {(Subscript[a, 0]), (Subscript[a, 1])} and their uncertainties {(Subscript[σ, 0]), (Subscript[σ, 1])}.  Thus, for this example we obtain

{pz, dpz} = {Parameters, Uncertainties} /. 
          zresults;
{{Exp[pz[[1]]], -1/
          pz[[2]]}, {Exp[
            pz[[1]]] 
          dpz[[1]], 
        dpz[[2]]/
          pz[[2]]^2}}

and recover the model parameters within their estimated uncertainties.  Note, however, that an unweighted fit

linfit[Transpose[{t, z}], {1, q}, q]

would fail badly.

Below we compare the fit to data either in linear

DisplayTogether[
ErrorListPlot[
          Transpose[{t, y, dy}]], 
Plot[
          Evaluate[Exp[Parameters . {1, tt}] /. zresults], {tt, 0, 1}], 
        PlotRange -> All, Frame -> True, 
        FrameLabel -> {"t", "y"}];

or semilog formats.  Unfortunately, although the semilog format is preferable for this example, Mathematica still lacks a method for plotting error bars on semilog plots.

DisplayTogether[
LogListPlot[
          Evaluate[Transpose[{t, y}]], 
          PlotStyle -> PointSize[0.02]], 
LogPlot[
          Evaluate[Exp[Parameters . {1, tt} /. zresults]], {tt, 0, 1}], 
        PlotRange -> {Automatic, {-1, 3.2}}, Frame -> True, 
        FrameLabel -> {"t", "y"}];

Including uncertainties in the independent variable

Often measurement of the independent variable is also subject to nonnegligible random errors, but most standard software for linear least-squares does not accommodate uncertainties in both variables simultaneously.  For example, the LinearRegression package in Mathematica provides many sophisticated options, but the use of uncertainties in more than one variable is not among them.  Therefore, we need to develop our own fitting procedure.

Suppose that our data set takes the form \(TraditionalForm\`{{\(x\_\(\(1\)\(,\)\)\) δ x\_1, 
          y\_1, δ y\_1}, ⋯, {x\_N, δ x\_N, y\_N, δ y\_N}}\) while the model is given by y[
          x]  =  (∑_(μ = 1)^m) (Subscript[a, μ]) (Subscript[f, μ])[x].  The measurements for both variables

\(TraditionalForm\`x\_i = 
      x\& . \_i + δ x\_i\ , \ \ \ y\_i = 
      y\& . \_i + δ y\_i\)

differ from the true values, indicated by dots, by random fluctuations based upon normal distributions with variances \(TraditionalForm\`{δ x\_i, δ y\_i}\).  If the model parameters were known, the expected (Subscript[y, i]) could be computed from the expected (Subscript[x, i]) according to

\(TraditionalForm\`y\& . \_i\  = \ 
      y[x\& . \_i]\  ≈ \ 
        y[x\_i]\  - \ \(\((∂y\/∂x)\)\_\(x = 
                  x\_i\)\) δ x\_i\)

where the final term is the first-order change due to the fluctuation in x.  The deviation in y then has two contributions

\(TraditionalForm\`y\_i - 
        y\_i[x\_i]\  ≈ \ δ y\_i\  - \ \(\((∂y\/∂x)\)\_\(x = 
                x\_i\)\) δ x\_i\)

where δ (Subscript[x, i]) and δ (Subscript[y, i]) are both random variables.  Therefore, the variance for a measurement is given by

\(TraditionalForm\`σ\_i\^2 ≈ \ \((δ y\_i\ )\)\^2 + \ \((\(\((∂y\/∂x)\)\_\(x = x\_i\)\) δ x\_i)\)\^2\)

Application of this result requires iteration.  First we fit without using the information about δ (Subscript[x, i]) and obtain a set of parameters {(Subscript[a, μ]^((1)))}.  These parameters are used to compute ∂y/∂x and the contributions to (Subscript[σ, i]) from δ (Subscript[x, i]).  Then we refit and obtain new parameters {(Subscript[a, μ]^((2)))}.  This set is then used to re-evaluate {(Subscript[σ, i])} and the procedure is repeated until the changes in all of the parameters become acceptably small or until our patience is exhausted and a maximum number of iterations is reached.

This iterative procedure is implemented in the linfit2 program listed below.  The fits are made with linfit and the rest of the code adds the uncertainties due to the two variables and tests for convergence.

linfit2::usage="linfit2[data,f,var] performs linear least-squares fitting to data with one independent and one dependent variable with uncertainties in both.  The data array should be a list of {x,dx,y,dy} for each measurement.  The list of basis functions is f and var is the symbol used in f for the independent variable.";
linfit2::BadDim = "Bad dimensions in \
data array: `1`";
linfit2::conv = "Warning: failed to converge.  Maximum \
relative change = `1`.";

Options[linfit2] = {MaxNITS -> 10, Conv -> 0.001};

linfit2[data_, 
        f_?VectorQ, var_, opts___Rule] := 
      Module[{x, dx, y, dy, σ, ndata, nparam, ndim, i, j, params, 
          slopes, fp, errors, yfit, redchi, oldparams, nits, mnits, change, 
          conv}, 

 (* Extract x, dx, 
          y, dy vectors from data array. *) 

{ndata, ndim} = 
          Dimensions[data]; 
If[
          ndata ≤ Length[f] || 
            ndim ≠ 4, (Message[linfit2::BadDim, {ndata, ndim}]; 
            Return[])]; 
{x, dx, y, dy} = 
          Transpose[data]; 

 (* 
          Derivative of basis functions. *) 

nparam = 
          Length[f]; 
slopes = 
          Table[0, {ndata}]; 
fp = (∂_var) 
            f; 

 (* 
          Loop using weights that include slope of fit function. *) 

{mnits, 
            conv} = {MaxNITS, Conv} /. {opts} /. 
            Options[linfit2]; 
oldparams = 
          Table[0, {i, 
              nparam}]; 

For[(nits = 
              0; change = 1), (nits < mnits && 
              change > 
                conv), nits++, 
σ = (Sqrt[(dy^2) + (((slopes*dx))^2)]); 
{params, 
              errors} = {Parameters, Uncertainties} /. 
              linfit[Transpose[{x, y, σ}], f, 
                var]; 
yfit = 
            params . f /. var -> x; 
slopes = 
            params . fp /. var -> x; 
change = 
            Max[Abs[(params - oldparams)/
                  params]]; 
oldparams = 
            params]; 
If [mnits > 1 && change > conv, 
          Message[linfit2::conv, 
            NumberForm[change, 
              2]]]; 
redchi = (σ^(-2)) . (((y - yfit))^2)/(ndata - nparam); 

 (* 
          Return results as replacement rules. *) 

{Parameters -> params, 
          Uncertainties -> errors, ReducedChiSquare -> redchi, 
          NITS -> nits}]

?linfit2

example: fit Legendre coefficients to data with error bars in both variables

Consider a model of the form

y[
        x]  =  (∑_μ) (Subscript[a, μ]) (Subscript[P, μ])[x]

where (Subscript[P, μ]) is the Legendre polynomial with degree μ.  The model is linear with respect to the coefficients (Subscript[a, μ]), so data can be fit using the method of linear least squares.  Below we construct and fit a sample data set with Gaussian fluctuations in both x and y.  Note that we assumed that the relative uncertainty in y is uniform.

Clear[z];
p = {3.02, 7.3, -9.2};
f = LegendreP[{0, 1, 2}, z];
x = Range[-0.95, 0.95, 0.05];
dx = Table[0.05, {Length[x]}];
y = p . f /. z -> x;
dy = 0.1 Abs[y];
x = 
        x + (Random[NormalDistribution[0, #]] & /@ 
              dx);
y = 
        y + (Random[NormalDistribution[0, #]] & /@ 
              dy);
data = Transpose[{x, dx, y, dy}];
results = 
      linfit2[data, f, z, MaxNITS -> 10, 
        Conv -> (10^(-4))]

DisplayTogether[
        MultipleListPlot[
          Evaluate[
            Transpose[{x, dx, y, dy}] /. {xx_, dxx_, yy_, 
                  dyy_} -> {{xx, yy}, ErrorBar[dxx, dyy]}], 
          PlotRange -> All, 
          PlotStyle -> PointSize[0.02]], 
Plot[
          Evaluate[Parameters . f /. results], {z, -1, 1}], 
        Frame -> True, 
        FrameLabel -> {"x", "y"}];

Nonlinear least squares

Often the model cannot be transformed into a linear function of its parameters and we must employ more sophisticated algorithms to search for the minimum in (χ^2)[a].  Mathematica provides a very powerful function for this purpose, NonlinearRegress, but a serious experimentalist must always be prepared to develop a custom search procedure whenever the problem at hand does not abide by the limitations of this function or other commercial software.  For example, uncertainties in the independent variable are not considered by NonlinearRegress.  Furthermore, it is more instructive to write our routine than to use a black box, even if that black box provides more sophisticated output.  Therefore, we develop and implement a nonlinear least squares algorithm before discussing the built-in version.

Levenberg-Marquardt algorithm

The basic problem is to locate the minimum in a function of several variables, (χ^2)[a], which can be described as a surface in a multidimensional space.  Many strategies and programs can be found in the literature, but we will focus upon the Levenberg-Marquardt algorithm — it is relatively simple but I have found that it provides a very effective method for a wide variety of physics problems.  This algorithm uses a single parameter, λ, to smoothly combine two distinct strategies.  Imagine that you are a hiker in a mountain range seeking the comfort of a lodge locating in the deepest valley.  If we start our search at location on the (χ^2) surface with a steep slope, the most effective strategy would be to head directly downhill.  On the other hand, if our current location is already near a minimum the slope is small and the gradient search becomes inefficient.  On the other hand, the shape of the surface near a minimum can be well approximated by parabolic function of all variables and we can use that information to locate the minimum accurately.  Iteration is still needed because we can only approximate the derivatives at the minimum by those at our current location, but the number of iterations required to obtain the desired precision is much smaller than needed for the gradient method.  The combination of these methods, gradient where the slope is large and parabolic where it is small, will almost always find a minimum from almost any starting point.  However, we cannot guarantee that this minimum is actually the global minimum that we seek — the (χ^2) surface for a multidimensional search can be quite complicated and many contain many local minima.  If the reduced chi-square for our minimum is substantially greater than unity, we should probably vary the starting conditions and repeat the fit several times to determine whether there is a deeper minimum.  Although we will not implement this extension of the method, the most robust methods often perform fits on a grid of initial conditions in order to maximize the chances of finding the true global minimum.  If we understand the physics of our problem relatively well already and can obtain reasonably good starting conditions by combining insight with visual inspection of the data, the Levenberg-Marquardt algorithm will usually provide the correct minimum without more laborious testing of local minima.  I have found that it is usually fairly tolerant of poor starting conditions — they should be reasonable but need not be especially precise.

Let (Subscript[χ, 0]^2) = (χ^2)[(a^((0)))] represent the value of chi-square at the starting point (a^((0))) = {(Subscript[a, μ]^((0))), μ = 1, m}.  The shape of the surface can be approximated by the Taylor expansion

\(TraditionalForm\`χ\^2\  ≈ \ χ\_0\%2\  + \ ∑\+ν∂χ\_0\%2\/∂a\_ν\ Δ a\_ν\)

where Δ (Subscript[a, ν]) = 
        (Subscript[a, ν]) - (Subscript[a, ν]^((0))) is the displacement in the parameter with index ν relative to its initial value.  If we assume that this approximation is valid for a sufficiently large region to include a minimum, we could locate that minimum by solving the normal equations

\(TraditionalForm\`∂χ\^2\/∂a\_μ == 
      0\ \ ⟹\ \ \ \ \(∂χ\_0\%2\/∂a\_μ\)\  + \ ∑\+ν∂\^2 χ\_0\%2\/\(∂a\_μ ∂a\_ν\)\ Δ a\_ν\  == \ 0\)

It is convenient to define

\(TraditionalForm\`β\_μ\  = \ \(\(-\(1\/2\)\) ∂χ\_0\%2\/∂a\_μ\  = \ ∑\+\(i = 1\)\%N\(\( w\_i\)(
              y - y\&_\_i)\) ∂y\&_\_i\/∂a\_μ\)\)

proportional to the gradient vector and

\(TraditionalForm\`α\_\(μ, ν\)\  = \ \(\(1\/2\) ∂\^2 χ\_0\%2\/\(∂a\_μ ∂a\_ν\)\  = \ ∑\+\(i = 1\)\%N\( 
                w\_i\)(\(∂y\&_\_i\/∂a\_μ\) ∂y\&_\_i\/∂a\_ν - \((y\_i - 
                    y\&_\_i)\) ∂\^2 y\&_\_i\/\(∂a\_μ ∂a\_ν\))\)\)

as the curvature matrix, such that the normal equations take the form

\(TraditionalForm\`β\  == \ α\ Δ a\ \ \ ⟹\ \ Δ a\  = \ ε\ β\)

where ε = (α^(-1)) is the error matrix.  Ideally, the new parameters

a  =  
      (a^((0))) + Δ a

should be closer to the minimum than the initial parameters, but are not likely to be exact because our representation of (χ^2) was only approximate.  Even if the direction of Δ a is good, (χ^2) might actually increase if the step is too large.  Therefore, we should compare (χ^2)[a] with (χ^2)[(a^((0)))] and reduce the step size if (χ^2) increases.  If (χ^2) decreases, we replace (a^((0))) with a and repeat the procedure, treating the new location as the starting point for the next iteration.

The Levenberg-Marquardt algorithm improves upon this procedure by modifying the curvature matrix.  If α happened to be a diagonal matrix, then each Δ (Subscript[a, μ]) would only depend upon the corresponding derivative (Subscript[β, μ]), with a scale factor determined by (Subscript[ε, μ, μ]) -> 
        1/(Subscript[α, μ, μ]).  The search path then becomes similar to the gradient search.  This modification is accomplished by the simple expedient of multiplying the diagonal elements of the curvature matrix by a scale factor

\(TraditionalForm\`α\&~\_\(μ, ν\) = \ \((1 + λ\ δ\_\(μ, ν\))\) α\_\(μ, ν\)\ , \ \ \ ε\&~ = \ α\&~\^\(-1\), \ \ \ \ Δ a\  = \ ε\&~\ β\)

If λ is large the diagonal dominates and the search is of gradient type with relatively small step, while for small λ the search is essentially parabolic.  Assuming that the user will supply reasonably good initial conditions, we start with small λ, say λ = 0.001.  Each time we evaluate Δ a we adjust λ: if (χ^2) decreases we take the step Δ a and reduce λ by a factor of 10, but if (χ^2) increases we skip that step, increase λ by a factor of 10, and try again.  Increasing λ reduces the step size and tends to align its direction with the gradient.

Covariance and correlation matrices

Rigorous analysis of the uncertainties in fitted parameters becomes rather difficult for nonlinear least squares because we do not have an analytical relationship between the fitted parameters and the data with their uncertainties.  Nor can one be certain that the fitted parameters are distributed in a Gaussian fashion.  The distribution of fitted parameters can be simulated by pretending that the results obtained by analyzing an experimental data set are correct and using those parameters, together with a random number generator, to produce many sets of pseudodata containing the same set of measurements as the experimental data but with random fluctuations in the results.  The distribution of fluctuations in the parameters fitted to pseudodata can then by used to analyze the uncertainties in and the correlations between fitted parameters.  However, this type of analysis is too laborious and too sophisticated for the present purposes.  Simpler estimates can be obtained by assuming that the fitting function can be linearized within a small region of the parameter space that is centered upon \(TraditionalForm\`a\&^\).

Being optimistic, we assume that the uncertainties in the fitted parameters are small and that their values are close to the true values.  We can then linearize the fitting function according to

\(TraditionalForm\`y[x; a]\  ≈ \ 
      y[x; a\&^]\  + \ \(\(\(\ \)\(∑\)\)\+μ\) \(\((∂y\/∂a\_μ)\)\_\(a = 
                a\&^\)\) Δ a\_μ\)

where \(TraditionalForm\`Δ a = a - a\&^\).  If we were to fit other experimental data sets from the same parent distribution using this linearized version of the model, we would expect the deviations δ (Subscript[a, μ]) to be normally distributed about 0 and could interpret the variance in the distribution of deviations as the standard error in (Subscript[a, μ]).  However, we already know that the variance in the parameters for linear least-squares is given by

\(TraditionalForm\`σ[\(a\&^\)\_μ]\^2\  = \ ε\_\(μ, μ\)\)

and expect that this relationship remains valid whenever Δ a is small enough to employ the linear approximation to the fitting function.

Similarly, near the minimal (χ^2) we expand

\(TraditionalForm\`\(χ\^2\)[
          a]\  ≈ \ \(χ\^2\)[
            a\&^]\  + \ \(1\/2\) \(∑\+\(μ, ν\)\(\((∂\^2 χ\^2\/\(∂a\_μ ∂a\_ν\))\)\_\(a = 
                      a\&^\)\) Δ \(a\_μ\) Δ a\_ν\)\  = \ \(χ\^2\)[
          a\&^]\  + ∑\+\(μ, ν\)\ Δ \(a\_μ\) \(α\_\(μ, ν\)\) Δ a\_ν\)

where the first-order contribution vanishes.  Remembering that (χ^2) is closely related to the joint probability distribution for the experiment, we write

\(TraditionalForm\`P\  ∝ \ 
      Exp[\(-\(χ\^2\)[a]\)]\  ∝ \ 
      Exp[\(-\(∑\+\(μ, ν\)\ Δ \(a\_μ\) \(α\_\(μ, ν\)\) Δ a\_ν\)\)]\)

in the form of a multivariate Gaussian distribution and identify ε = (α^(-1)) as the covariance matrix.  Thus, the statistical correlation between fitted parameters is given by

\(TraditionalForm\`C\_\(μ, ν\)\  = \ \(⟨Δ \(a\_μ\) Δ a\_ν⟩\/\@\(⟨\((Δ a\_μ)\)\^2⟩ ⟨\((Δ a\_ν)\)\^2⟩\)\ \  = \ ε\_\(μ, ν\)\/\@\(\(ε\_\(μ, μ\)\) ε\_\(ν, ν\)\)\)\)

where we leave the evaluation of the integrals as an exercise for the reader.   The diagonal elements have unit values while nondiagonal elements are limited to the range (-1, 1).  If two parameters are uncorrelated, the product \(TraditionalForm\`Δ \(a\_μ\) Δ a\_ν\) exhibits fluctuations with either positive or negative signs with equal likelihood such that (Subscript[C, μ, ν]) averages to zero.  If a fluctuation in (Subscript[a, μ]) is more likely to accompanied by a fluctuation in (Subscript[a, ν]) that has the same sign, then (Subscript[C, μ, ν]) > 0 represents a positive correlation between those parameters whereas a negative correlation with (Subscript[C, μ, ν]) < 0 indicates that a fluctuation in Δ (Subscript[a, μ]) is more likely to be compensated by a fluctuation in Δ (Subscript[a, ν]) with opposite sign.

Implementation

The program below implements the Levenberg-Marquardt algorithm for problems in which there is a single independent variable with negligible uncertainty and a single dependent variable.  We also assume that derivatives of the fitting function with respect to each parameter can be obtained symbolically and we neglect the contributions of second derivatives to the curvature matrix.  This approximation is reasonable when the fit is good, such that the coefficient of the second derivative is small, fluctuates with random sign and tends to average to zero.  Furthermore, experience has shown that neglecting the second derivatives often improves the performance of the algorithm when the initial conditions are inaccurate or the data contain some points with greater than normal deviations (outliers).  Finally, evaluation of second derivatives can be computationally expensive.  Therefore, it is customary to omit the second derivatives, but we provide an option for their inclusion.  In the interests of simplicity, we use a simple For loop and have not attempted to minimize the number of evaluations of β and α (re-evaluation is not needed when λ is increased).  It should not be difficult to modify this procedure to accommodate multiple independent or dependent variables or to embed this procedure in a loop that accounts for uncertainties in the independent variable(s), but we leave those extensions as an exercise to the reader that needs them.

curfit::usage="curfit[data,f,var,par,start] performs nonlinear least-squares fitting using the Levenberg-Marquardt algorithm. The data array should be a list of {x,y,dy} for each measurement.  The fitting function f is expressed in terms of the independent variable var and a list of parameters par. The initial conditions are given as a list start.";curfit::BadDim = "Bad dimensions in data array: `1`";
curfit::conv = "Warning: failed to \
converge.  Relative change = `1`; λ = `2`";
curfit::start = "Incorrect number of initial values.";

Options[curfit] = {MaxNITS -> 10, Conv -> 0.001, 
            PrintIterations -> False, 
            IncludeSecondDerivatives -> False};

curfit[data_, f_, 
        var_, par_, start_, opts___Rule] := 
      Module[{npar, ndata, ndim, nits, a, da, λ, x, y, dy, w, yfit, 
          derivs, dmat, deriv2, dmat2, β, α, ε, 
          errors, corr, redchi, change, conv, mnits, i, j, oldchi, newchi, 
          printIterations, 
          includeSecondDerivatives}, 
npar  =  
          Length[par]; 

 (* Extract x, 
          y, dy vectors from data array. *) 

{ndata, ndim} = 
          Dimensions[data]; 
If[
          ndata ≤ npar || ndim < 2 || 
            ndim > 3, (Message[curfit::BadDim, {ndata, ndim}]; 
            Return[])]; 
If[
          Length[start] ≠ npar, (Message[curfit::start]; 
            Return[])]; 
If[
          ndim == 2, 
({x, y} = Transpose[data]; 
            w = Table[1, {ndata}]), 
({x, y, dy} = 
              Transpose[data]; 
            w = 1/(dy^2))]; 

 (* 
          Initialize derivatives and search parameters. *) 

{mnits, conv, printIterations, 
            includeSecondDerivatives} = {MaxNITS, Conv, PrintIterations, 
                IncludeSecondDerivatives} /. {opts} /. 
            Options[curfit]; 

derivs = 
          Table[Simplify[
              D[f, par[[i]]]], {i, 
              npar}]; 
If[
          includeSecondDerivatives, 
deriv2 = 
            Table[Simplify[
                D[f, par[[i]], 
                  par[[j]]]], {i, 
                npar}, {j, npar}]]; 
change = 
          1. ; 
λ = 
          0.001; 

 (* 
          Search for minimum. *) 

For[(nits = 0; 
            a = start), (nits ≤ mnits && 
              change > conv), nits++, 
yfit = f /. 
                Thread[par -> a] /. 
              var -> x; 
oldchi = 
            w . (((y - yfit))^2); 
dmat = 
            Table[derivs /. Thread[par -> a] /. 
                var -> x[[i]], {i, 
                ndata}]; 
β = (w (y - 
                      yfit)) . dmat; 
α = 
            Transpose[w dmat] . dmat; 
If[
            includeSecondDerivatives, 
dmat2 = 
              Table[deriv2 /. Thread[par -> a] /. 
                  var -> 
                    x[[i]], {i, 
                  ndata}]; 
α = α - (w (y - yfit)) . dmat2]; 
Do[α[[i, 
                i]] = (1 + λ) α[[i, i]], {i, 
              npar}]; 
ε = 
            PseudoInverse[α]; 
da = ε . β; 
yfit = f /. 
                Thread[par -> (a + da)] /. 
              var -> x; 
newchi = 
            w . (((y - yfit))^2); 
If[
            newchi < oldchi, (λ = λ/10; a = a + da; 
              change = (oldchi - newchi)/oldchi), λ = 
              10  λ]; 
If[printIterations, 
            Print[{oldchi, newchi, λ, change, 
                a + da}]];
]; 
If [
          mnits > 1 && change > conv, 
          Message[curfit::conv, NumberForm[change, 2], 
            NumberForm[λ, 
              2]]]; 

 (* 
          Evaluate proper error matrix. *) 

α = 
          Transpose[w dmat] . dmat; 
If[
          includeSecondDerivatives, 
dmat2 = 
            Table[deriv2 /. Thread[par -> a] /. 
                var -> x[[i]], {i, 
                ndata}]; 
α = α (w (y - 
                        yfit)) . 
                dmat2]; 
ε = 
          PseudoInverse[α]; 
errors = 
          Table[(Sqrt[ε])[[i, i]], {i, npar}]; 
corr = 
          Table[ε[[i, j]]/(errors[[i]] errors[[j]]), {i, npar}, {j, 
              npar}]; 
yfit = f /. 
              Thread[par -> a] /. 
            var -> x; 
redchi = 
          w . (((y - yfit))^2)/(ndata - 
                npar); 

 (* 
          Return results as replacement rules. *) 

{Parameters -> a, 
          Uncertainties -> errors, CorrelationMatrix -> corr, 
          ReducedChiSquare -> redchi, 
          NITS -> nits}
]

Performance tests

Below we construct pseudodata for a radioactive decay experiment

Clear[t, a, τ];
tdata = Range[0, 1, 0.02];
p = {(10^3), 0.1};
yfit[t_] := 
        p[[1]] 
          Exp[-t/
              p[[2]]];
ydata = Random[PoissonDistribution[yfit[#]]] & /@ tdata // 
          N;
dy = (Sqrt[1.  + ydata]);

and attempt to fit, with reasonable but inaccurate initial conditions, including second derivatives.

results = 
      curfit[Transpose[{tdata, ydata, dy}], a Exp[-t/τ], 
        t, {a, τ}, {1500, 0.7}, MaxNITS -> 20, 
        PrintIterations -> True, 
        IncludeSecondDerivatives -> True]

Although this search fails to converge, simply removing the contribution to the curvature matrix due to second-derivatives of the fitting function

results = 
      curfit[Transpose[{tdata, ydata, dy}], a Exp[-t/τ], 
        t, {a, τ}, {1500, 0.7}, MaxNITS -> 20, 
        PrintIterations -> True]

produces a successful fit even though it was necessary to increase λ several times to compensate for the poor initial conditions.  Starting with the gradient search gets the fit going in the right direction, and then it settles into the minimum quickly.  Experience shows that this algorithm is quite robust!  With relatively simple modifications, one can include constraints upon the parameters that prevent the fit from attempting to move in unphysical directions.  For this problem, τ < 0 represents exponential growth and the initial difficulties in the search path arise when poor initial conditions result in parameter increments that venture into unphysical territory.  Adding (χ^2) a penalty for unphysical parameters could mitigate this problem.  We do not investigate those techniques here but simply re-iterate that the use of one's own code instead of generic commercial software provides that flexibility.

example

The scattering of  high-energy electrons by protons or neutrons is governed by electric and magnetic form factors which nonrelativistically are given by Fourier transforms of charge and magnetization densities.  Until recently it was believed that the proton electric form factor, (Subscript[G, Ep]), was approximately equal to the so-called dipole form factor, (Subscript[G, D]), such that

(Subscript[G, Ep])[(Q^2)]  ≈  (Subscript[G, D])[(Q^2)]  =  
      (1/(((1 + (Q^2)/(Λ^2)))^2))

where (Q^2) is the square of the momentum transferred in the electron-proton collision and (Λ^2) = 0.71 (((GeV/c))^2) is the characteristic scale related to the size of the proton.  However, more recent data using a new technique shows that the ratio (Subscript[G, Ep])/(Subscript[G, D]) decreases fairly rapidly for (Q^2) > 1.  A simple parametrization is the rational function of (Q^2)

\(TraditionalForm\`G\_Ep[
        Q\^2]\  = \ \(∑\+\(μ = 0\)\%m\( a\_μ\) Q\^\(2  μ\)\)\/\(∑\+\(ν = 0\)\%\(m + 2\)\(b\_ν\) Q\^\(2  ν\)\)\)

where the symmetry of the problem limits the expansions to even powers.  Recognizing that (Subscript[G, Ep])[0] = 1 because the charge must be one elementary unit, we can constrain (Subscript[a, 0]) = (Subscript[b, 0]) = 1, but if may be useful to vary (Subscript[a, 0]) if the systematic normalization uncertainties are nonnegligible.  Finally, there are theoretical reasons based upon the quark model to require

\(TraditionalForm\`Q\^2 -> ∞\ \ ⟹\ \ G\_Ep ∝ Q\^\(-4\)\)

This is the reason that we include two more terms in the denominator than in the numerator.

We are now ready to analyze some data.  The following cell reads a data file assuming that it is located in a directory contain in $Path; if that is not true, you will need to include the proper directory in the file name.

gepdata = Import["gep.dat"];

Next we try to fit the simplest version of this model

curfit[gepdata, (1/(1 + b1 qq  +  b2 (qq^2))), 
      qq, {b1, b2}, {1.4, 0}, MaxNITS -> 20]

but find the results quite unsatisfactory.  Therefore, we add another term to both numerator and denominator and obtain a satisfactory fit.

gepresults = 
      curfit[gepdata, ((1 + a1 qq)/(1 + b1 qq  +  b2 (qq^2) + b3 (qq^3))), qq, {a1, b1, b2, b3}, {0, 1.4, 0, 0}, MaxNITS -> 20]

It is useful to display the fit in the form of (Subscript[G, Ep])/(Subscript[G, D]).  Clearly, this parametrization describes the data very well.  The fact that (Subscript[χ, red]^2) is significantly below 1 for a relatively large data set combining results from several different experiments suggests that some of the experimentalists overestimated the random uncertainties in their data or included systematic uncertainties in the published error bars.

{x, y, dy} = Transpose[gepdata];
y = y (((1 + (x/0.71)))^2);
dy = dy (((1 + (x/0.71)))^2);
ratiodata = Transpose[{x, y, dy}];
\(\(DisplayTogether[
        Plot[Evaluate[\(\(1 + a1\ qq\)\/\(1 + b1\ qq\  + \ b2\ qq\^2 + 
                      b3\ qq\^3\)\) \((1 + qq\/0.71)\)\^2 /. 
              Thread[{a1, b1, b2, 
                    b3} -> \((Parameters /. gepresults)\)]], {qq, 0, 
            10}, PlotRange -> {{0, 10}, {0, 1.2}}, Frame -> True, 
          FrameLabel -> {\*""(Q^2)"", \*""(Subscript[G, Ep])/(Subscript[G, D])""}], 
ErrorListPlot[
          ratiodata]];\)\)

NonlinearRegress

Mathematica provides several functions for fitting data, but here we discuss only NonlinearRegress because it is the most general and subsumes the other functions as special cases.  The basic syntax and list of options are given in the reference guide, obtained by following the highlighted link.  We illustrate some of the options by returning to an earlier example.  For use in this function, we provide data in the form {{(Subscript[x, 1]), (Subscript[y, 1])}, ⋯, {(Subscript[x, N]), (Subscript[y, N])}} and a separate list of weights \(TraditionalForm\`{σ\_1\^\(-2\), ⋯, σ\_N\^\(-2\)}\).  The fitting function is given explicitly, followed by a list of independent variables (here only 1) and a list of parameters.  The output is returned as a list of replacement rules.

x = Range[-0.9, 0.9, 0.1];
y = 
        117  -  23  x + 15  (x^2) + 
          Table[Random[
              NormalDistribution[0, 2. ]], {Length[
                x]}];
dy = Table[2. , {Length[y]}];
data = Transpose[{x, y}];
w = 1/(dy^2);

Clear[a, b, c, z];
NonlinearRegress[data, a + b z + c (z^2), z, {a, b, c}, 
      Weights -> w]

Obviously, BestFitParameters returns the fitted parameters.  The ParameterCITable returns the standard errors (σ for each parameter) and a confidence interval that should include the true model parameters for 95% of data samples from the parent distribution.  The desired confidence level is an adjustable option.  The ANOVATable is a more advanced feature, standing for analysis of variance, that we do not discuss; nor do we discuss the curvature table.  The correlation matrix assesses correlations between fitted parameters with nondiagonal elements in the range (-1, 1).

The contents of the regression report can be expanded or compressed; the default report above is labelled a summary report.  The quantities that can be included in the regression report are listed below; hopefully, most of the names should be relatively clear.

RegressionReportValues[NonlinearRegress]

Below we use the regression report to compare the fit to the data and to display the residuals.

results = 
      NonlinearRegress[data, a + b z + c (z^2), z, {a, b, c}, 
        Weights -> w, 
        RegressionReport -> {BestFitParameters, ParameterCITable, 
            EstimatedVariance, BestFit, FitResiduals, 
            StandardizedResiduals}]

DisplayTogether[
        ErrorListPlot[Transpose[{x, y, dy}]], 
Plot[
          Evaluate[BestFit /. results], {z, -1, 1}], Frame -> True, 
        FrameLabel -> {"x", "y"}];

DisplayTogether[
        ListPlot[Evaluate[Transpose[{x, StandardizedResiduals /. results}]], 
          PlotRange -> {Automatic, {-3, 3}}, 
          PlotStyle -> PointSize[0.02], Frame -> True, 
          FrameLabel -> {"x", "residuals"}], 
        Plot[{-1, 1}, {x, -1, 1}, 
          PlotStyle -> Dashing[{0.02, 0.02}]]];

Notice that the estimated variance is equivalent to reduced chisquare.

w . ((((FitResiduals /. results))^2)/(Length[y] - 3))

Next we consider the fit to (Subscript[G, Ep]).  In order to provide initial values for the fitting parameters, we use a poorly documented feature of the syntax.  Because the fitting code employs FindMinimum, the fitting parameters can be listed in the form {{a, (Subscript[a, 0])}, {b, (Subscript[b, 0])}, ⋯}.  To use NonlinearRegress we also have to separate the uncertainties and compute the weights.

{x, y, dy} = Transpose[gepdata];
gepresults2 = 
      NonlinearRegress[
        Transpose[{x, 
            y}], ((1 + a1 qq)/(1 + b1 qq  +  b2 (qq^2) + b3 (qq^3))), 
        qq, {{a1, 0}, {b1, 1.4}, {b2, 0}, {b3, 0}}, 
        Weights -> (dy^(-2))]

{x, y, dy} = Transpose[gepdata];
y = y (((1 + (x/0.71)))^2);
dy = dy (((1 + (x/0.71)))^2);
ratiodata = Transpose[{x, y, dy}];

\(\(DisplayTogether[
        Plot[Evaluate[\(\(1 + a1\ qq\)\/\(1 + b1\ qq\  + \ b2\ qq\^2 + 
                      b3\ qq\^3\)\) \((1 + qq\/0.71)\)\^2 /. 
              Thread[{a1, b1, b2, 
                    b3} -> \((Parameters /. gepresults)\)]], {qq, 0, 
            10}, PlotRange -> {{0, 10}, {0, 1.2}}, Frame -> True, 
          FrameLabel -> {\*""(Q^2)"", \*""(Subscript[G, Ep])/(Subscript[G, D])""}], 
Plot[
          Evaluate[\(\(1 + a1\ qq\)\/\(1 + b1\ qq\  + \ b2\ qq\^2 + 
                      b3\ qq\^3\)\) \((1 + qq\/0.71)\)\^2 /. \((BestFitParameters /. gepresults2)\)], {qq, 0, 10}, 
          PlotStyle -> Hue[0.3]], 
ErrorListPlot[
          ratiodata]];\)\)

\(\(DisplayTogether[
        Plot[Evaluate[\(\(1 + a1\ qq\)\/\(1 + b1\ qq\  + \ b2\ qq\^2 + 
                      b3\ qq\^3\)\) \((1 + qq\/0.71)\)\^2 /. 
              Thread[{a1, b1, b2, 
                    b3} -> \((Parameters /. gepresults)\)]], {qq, 0, 
            10}, PlotRange -> {{0, 2}, {0.8, 1.2}}, Frame -> True, 
          FrameLabel -> {\*""(Q^2)"", \*""(Subscript[G, Ep])/(Subscript[G, D])""}], 
Plot[
          Evaluate[\(\(1 + a1\ qq\)\/\(1 + b1\ qq\  + \ b2\ qq\^2 + 
                      b3\ qq\^3\)\) \((1 + qq\/0.71)\)\^2 /. \((BestFitParameters /. gepresults2)\)], {qq, 0, 10}, 
          PlotStyle -> Hue[0.3]], 
ErrorListPlot[
          ratiodata]];\)\)

These results agree very well with those obtained using our custom program; you cannot distinguish the two curves.  For historical reasons, we note that the agreement was not as good with Mathematica versions prior to 5.2 because there was a bug in NonlinearRegress — while searching for the best fit, earlier versions minimized \(TraditionalForm\`∑\+i\(\( w\_i\^2\)(y\_i - y\&_\_i)\)\^2\) instead of \(TraditionalForm\`∑\+i w\_i\ \((y\_i - y\&_\_i)\)\^2\) yet computed the estimated variance based upon the latter instead of the former.  Fortunately, this bug has now been corrected, but there are still advantages to writing programs whose inner workings you understand well and can explore or extend, which is not so easy with commercial software.  For example, none of Mathematica's fitting functions handle uncertainties in the independent variable, which is sometimes important.  No software is perfect — you must be able to develop your own programs whenever faced with either limitations or reliability issues in commercial software!

The Curfit package

Now that we have functions that perform both linear and nonlinear least-squares fits, it is appropriate to package them for use in other notebooks.  Improved versions of both Linfit and Curfit are available in a package provided for this course.  Both functions accommodate data sets with 2 - 4 columns in the forms {x, y}, {x, y, (Subscript[σ, y])}, or {x, (Subscript[σ, x]), y, (Subscript[σ, y])}.  However, because we used many of the same symbols from the prototype versions presented in preceding sections, it is useful to start with a fresh slate.

ClearAll["Global`*"];
Remove["Global`*"]

Needs["Curfit`"]

?Curfit`*

The following subsections illustrate this package using some of the preceding examples.  In all cases we find good agreement with previous results, but the package provides extended features and simpler access to the relevant functions.  Feel free to use or modify this package for your future data analysis needs.

Legendre coefficients — compare Linfit and Curfit for a linear model

First we verify that the same results are obtained using either Linfit or Curfit when the fit function is linear in its parameters.  It is convenient to use the same Legendre expansion as above.  The subsections below demonstrate that this test succeeds whether or not uncertainties are included for the independent variable.  Note that version of Curfit in the package has been generalized to provide that feature.

error bars only for dependent variable

Clear[z];
p = {3.02, 7.3, -9.2};
f = LegendreP[{0, 1, 2}, z];
x = Range[-0.95, 0.95, 0.05];
y = p . f /. z -> x;
dy = 0.1 Abs[y];
y = 
        y + (Random[NormalDistribution[0, #]] & /@ 
              dy);
data = Transpose[{x, y, dy}];

results1 = Linfit[data, f, z]

start = 
        Table[{a[i], Random[]}, {i, 0, 2}];
results2 = 
      Curfit[data, Sum[a[i] LegendreP[i, z], {i, 0, 2}], z, 
        start]

Note that the initial guesses are irrelevant when the model is linear in its parameters; therefore, we just used random numbers.

error bars for both independent and dependent variables

Clear[z];
p = {3.02, 7.3, -9.2};
f = LegendreP[{0, 1, 2}, z];
x = Range[-0.95, 0.95, 0.05];
dx = Table[0.05, {Length[x]}];
y = p . f /. z -> x;
dy = 0.1 Abs[y];
x = 
        x + (Random[NormalDistribution[0, #]] & /@ 
              dx);
y = 
        y + (Random[NormalDistribution[0, #]] & /@ 
              dy);
data = Transpose[{x, dx, y, dy}];

results1 = Linfit[data, f, z]

start = 
        Table[{a[i], Random[]}, {i, 0, 2}];
results2 = 
      Curfit[data, Sum[a[i] LegendreP[i, z], {i, 0, 2}], z, 
        start]

DisplayTogether[
        MultipleListPlot[
          Evaluate[
            Transpose[{x, dx, y, dy}] /. {xx_, dxx_, yy_, 
                  dyy_} -> {{xx, yy}, ErrorBar[dxx, dyy]}], 
          PlotRange -> All, 
          PlotStyle -> PointSize[0.02]], 
Plot[
          Evaluate[Parameters . f /. results1], {z, -1, 1}], 
        Frame -> True, 
        FrameLabel -> {"x", "y"}];

GEp test

Finally, we reproduce the previous fit to (Subscript[G, Ep]) to illustrate the use of Curfit for a nonlinear model.

gepdata = Import["gep.dat"];

pade[
          Q2_] = ((1 + a[1] Q2)/(1 + b[1] Q2  +  b[2] (Q2^2) + b[3] (Q2^3)));
gepresults = 
      Curfit[gepdata, pade[Q2], 
        Q2, {{a[1], 0}, {b[1], 1.4}, {b[2], 0}, {b[3], 0}}]

{x, y, dy} = Transpose[gepdata];
y = y (((1 + (x/0.71)))^2);
dy = dy (((1 + (x/0.71)))^2);
ratiodata = Transpose[{x, y, dy}];
fit[Q2_] = 
        pade[Q2] (((1 + (Q2/0.71)))^2) /. (Parameters /. 
              gepresults);
\(\(DisplayTogether[
        Plot[fit[Q2], {Q2, 0, 10}, PlotRange -> {{0, 10}, {0, 1.2}}, 
          Frame -> True, 
          FrameLabel -> {\*""(Q^2)"", \*""(Subscript[G, Ep])/(Subscript[G, D])""}], 
ErrorListPlot[
          ratiodata]];\)\)