# Generalized Linear Least Squares Curve Fitting

## Summary

Least squares curve fitting can be generalized to work with any linear model to determine unique parameter estimates and uncertainties in terms of matrix equations. Linear models follow a general form. All other models are classified as nonlinear. While nonlinear models do not have a general analytic solution, some nonlinear models can be transformed into linear models. We give two examples of transforming an exponential model and a transformation of the enzymatic rate equation to straight line models using propagation of error. Polynomials form a class of linear models which we implement in MATLAB.

## Linear and nonlinear models

In least square curve fitting, there is a step where we take derivatives to find the parameter estimates. If those equations are too complicated then we will have a hard time solving them. It is straightforward to solve these equations for so called linear models. Linear models are of the form

One example of a linear model with parameters is a polynomial.

Other linear models might have the replaced with general functions . This gives a lot of different possible functions, but some are unnatural and do not occur in usual data sets.

Linear means the parameter is to the first power out in front. The following model is not linear.

For models of the nonlinear form, the optimization problem is usually quite difficult to solve. We will consider what to do with nonlinear models later.

To perform a curve fit on a model with parameters you need at least different points. A good practice is to have many more measured points than parameters. When models are bloated with a lot of parameters you can usually fit a wide variety of functions. Keep it as simple as possible but not simpler. Use only as few parameters are necessary.

## The generalized least squares problem

Now with our linear model of parameters we would like to derive the best parameter estimates and their uncertainties in the least squares weighted framework. The model is.

where are some generally defined continuous functions.

To minizise , we set each derivative equal to zero.

This gives a linear system of equations. These “normal equations” give the optimal fit parameters which we can write as a vector.

Taking terms on each side of the equal sign.

It is useful to do some linear algebra to solve these equations for the parameters and their uncertainties. The “design matrix,” , is defined as the row and column matrix

Second we define the column vector of rows.

A third matrix, has inverse .

Now the normal equations can be rewritten as

To solve for the fit parameters it is easy to numerically compute the inverse in MATLAB.

## General parameter uncertainties

Now, we can find the uncertainty in the parameters by the propagation of error formula.

But the Kroenecker delta because is the inverse of .

so the sum reduces to

These are the uncertainties in the fitting parameters.

The off diagonal terms of also give the covariances. If two different groups use the least squares theory on a linear model they should arrive at the same results for the fitting parameters and their uncertainties which are derived here. To summarize the results for linear least squares models

## Linearizing nonlinear models

Sometimes data following a nonlinear relationship may be transformed to a linear model. One can then estimate the transformed parameters exactly in conjunction with propagation of error. Here is a well known example from biochemistry for the velocity of a chemical reaction illustrating linearization.

This equation can be rewritten as

We start with as the set of information. is the standard deviation of the mean of the repeated measurement of . The transformed variables are

Propagation of error for a single variable has been used to find . One can then weighted fit to a line to find , , , and . and can then be found also by propagation of error. This transformation generates a linear plot called the Lineweaver-Burke plot.

Another straightforward example of linearization involves exponential decay.

One has the set of information where are the standard deviations of the mean.

Now take the logarithm of both sides.

Now we transform

If the error bars were originally similar for most points on the linear plot they will be stretched out at longer times after transformation. After finding and , you can use error propagation to find and its uncertainty.

Example. Fit some simulated exponentially decaying data with linear least squares by linearizing the model.

We generate some model data which are the average of five points of + normrnd(0,4).

clear;

ti = 0.0:0.1:2.5

N0 = 100;

L = 1;

Ni = 5;

for i = 1:length(ti)

temp = [];

for j = 1:Ni

temp(j) = N0*exp(-L*ti(i))+normrnd(0,4);

end

yi(i) = mean(temp);

ei(i) = std(temp)/sqrt(Ni);

end

figure(1);clf;

errorbar(ti,yi,ei,’.k’);

axis([-0.5 3 0 125]);

xlabel(‘x’);

ylabel(‘N’);

figure(2);clf;

zi = log(yi);

si = ei./yi;

errorbar(ti,zi,si,’.k’);

axis([-0.5 3 0 5])

xlabel(‘t’);

ylabel(‘z’);

hold on;

[A sigA B sigB chi2] = lineABx(ti, zi, si);

XX = 0:0.1:3

YY = A + XX.*B;

plot(XX,YY,’k’);

estN = exp(A)

estsigN = exp(A)*sigA

estL = -B

estsigL = sigB

figure(1);

hold on;

XX = 0:0.01:3

YY = estN.*exp(-estL.*XX);

plot(XX,YY,’k’);

Exponential decay with best fit

Linearization and weighted fit

## Least squares polynomial

Example. Implement a polynomial least squares fitting function in MATLAB.

Here is the MATLAB function we will use.

function [ A sigA ] = genpolyfit( M, xi, yi, sigi )

% [A sigA ] = genpolyfi( M, xi, yi, sigi)

% Least sqaures polynomial fit to degree M-1

%

for k = 1:M

X(k,:) = xi.^(k-1);

end

for j = 1:M

for k = 1:M

D(j,k) = sum(sigi.^(-2).*X(j,:).*X(k,:));

end

end

for k = 1:M

B(k) = sum(X(k,:).*yi.*sigi.^(-2));

end

A = inv(D)*B’;

C = inv(D);

for k = 1:M

sigA(k) = sqrt(C(k,k));

end

Example. Simulate some data that can be fit by a parabola.

We can use the least square polynomial function.

clear;figure(1);clf;

A = 1;

B = 2;

C = 3;

Ni = 5;

xi = -10:1:10;

for i = 1:length(xi)

temp = [];

for j = 1:Ni

temp(j) = A + B*xi(i) + C*xi(i).^2 + normrnd(0,10);

end

yi(i) = mean(temp);

si(i) = std(temp)/sqrt(Ni);

end

figure(1);clf;

errorbar(xi,yi,si,’.k’, ‘markersize’,10);

xlabel(‘x’)

ylabel(‘y’)

hold on;

[A sigA] = genpolyfit(3,xi,yi,si);

x1 = -11:0.1:11;

y1 = A(3).*x1.^2 + A(2).*x1 + A(1);

plot(x1,y1, ‘r’);

axis( [ -11 11 -50 1.1*max(y1)]);’

Curve fit to a parabola

Comments are Disabled