MATLAB:Fitting

This document contains examples of polynomial fitting, general linear regression, and nonlinear regression. In each section, there will be example code that may come in useful for later courses. The example code is based on the existence of a file in the same directory called Cantilever.dat that contains two columns of data - the first is an amount of mass (in kg) placed at the end of a beam and the second is a displacement, measured in inches, at the end of the beam. For EGR 53, this file is:

0.000000        0.005211
0.113510002     0.158707
0.227279999     0.31399
0.340790009     0.474619
0.455809998     0.636769
0.569320007     0.77989
0.683630005     0.936634
0.797140015     0.999986

Polynomial Fitting

Example Code

%% Initialize the workspace
clear; format short e
figure(1); clf

%% Load and manipulate the data
load Cantilever.dat

Force = Cantilever(:,1) * 9.81;
Displ = Cantilever(:,2) * 2.54 / 100;

%% Rename and create model data
x = Force;
y = Displ;
xmodel = linspace(min(x), max(x), 100);

%% Determine the polynomial coefficients
N = 1;
P = polyfit(x, y, N)

%% Generate estimates and model
yhat   = polyval(P, x);
ymodel = polyval(P, xmodel);

%% Calculate statistics
% Compute sum of the squares of the data residuals
St = sum(( y - mean(y) ).^2)

% Compute sum of the squares of the estimate residuals
Sr = sum(( y - yhat ).^2)

% Compute the coefficient of determination
r2 = (St - Sr) / St

%% Generate plots
plot(x,      y,      'k*',...
     x,      yhat,   'ko',...
     xmodel, ymodel, 'k-');
 xlabel('Independent Value')
 ylabel('Dependent Value')
 title('Dependent vs. Independent and Model')
 legend('Data', 'Estimates', 'Model', 0)

General Linear Regression

Example Code

%% Initialize the workspace
clear; format short e
figure(1); clf

%% Load and manipulate the data
load Cantilever.dat

Force = Cantilever(:,1) * 9.81;
Displ = Cantilever(:,2) * 2.54 / 100;

%% Rename and create model data
x = Force;
y = Displ;
xmodel = linspace(min(x), max(x), 100);

%% Define model equation and A matrix
model = 'linear'
switch model
    case 'linear'
        yeqn = @(coefs, x) coefs(1)*x.^1 + coefs(2)*x.^0;
        A    =                     [x.^1            x.^0];
    case 'quadratic'
        yeqn = @(coefs, x) coefs(1)*x.^2 + coefs(2)*x.^1 + coefs(3)*x.^0;
        A    =                     [x.^2            x.^1            x.^0];
    case 'line through origin'
        yeqn = @(coefs, x) coefs(1)*x.^1;
        A    =                     [x.^1];
    case 'trig'
        yeqn = @(coefs, x) coefs(1)*cos(x) + coefs(2)*sin(x);
        A    =                     [cos(x)            sin(x)];
    case 'trig with offset'
        yeqn = @(coefs, x) coefs(1)*cos(x) + coefs(2)*sin(x) + coefs(3)*x.^0;
        A    =                     [cos(x)            sin(x)            x.^0];
    otherwise
        error('Don''t know the model...')
end

%% Determine the function coefficients
MyCoefs = A\y

%% Generate estimates and model
yhat   = yeqn(MyCoefs, x);
ymodel = yeqn(MyCoefs, xmodel);

%% Calculate statistics
% Compute sum of the squares of the data residuals
St = sum(( y - mean(y) ).^2)

% Compute sum of the squares of the estimate residuals
Sr = sum(( y - yhat ).^2)

% Compute the coefficient of determination
r2 = (St - Sr) / St

%% Generate plots
plot(x,      y,      'k*',...
     x,      yhat,   'ko',...
     xmodel, ymodel, 'k-');
 xlabel('Independent Value')
 ylabel('Dependent Value')
 title('Dependent vs. Independent and Model')
 legend('Data', 'Estimates', 'Model', 0)

Linearized Models

There are three primary nonlinear models which, through a transformation of variables, may yield a linear relationship between the transformed variables: the exponential model, the power-law model, and the saturation-growth model. In the subsections before, these are addressed individually by showing the modelling equation, describing a transformation mechanism, and showing the transformed equation in the form of

\( \eta=P(1)*\xi^1+P(2)*\xi^0\,\! \)

where \(\xi\) is the transformed version of the independent variable \(x\) and \(\eta\) is the transformed version of the dependent variable \(y\). In each of the three cases below, then, the polyfit command can be used, with the transformed variables, to find the coefficients of the straight-line fit.

Exponential Model

Exponential models come from relationships where the derivative of the dependent variable is proportional to the value of the dependent variable itself. That is:

\(\begin{align} \frac{dy}{dt}&=kt \end{align}\)

Integrating, this yields the following model:

\(y=y_0e^{kx}\,\!\)

where \(y_0\) is the value of the dependent variable when the independent variable is 0 and \(k\) is the growth rate. Transform by taking the logarithm of each side. Using the natural logarithm is most useful:

\(\begin{align} y&=y_0e^{kx}\\ \ln\left(y\right)&=\ln\left(y_0e^{kx}\right)\\ \ln\left(y\right)&=k*x^1+\ln\left(y_0\right)*x^0\,\! \end{align}\)

meaning the transformed measurement \(\ln\left(y\right)\) can be fit to a straight line function of \(x\). That is, using the transformations:

\(\begin{align} \xi&=x & \eta&=\ln\left(y\right) \end{align}\)

The relationship between the slope and intercept for the transformed variables and the constants in the model equation are:

\(\begin{align} P(1)&=k & P(2)&=\ln\left(y_0\right)\\ k&=P(1) & y_0&=e^{P(2)} \end{align}\)

Power Law Model

Power Law models work for several different phenomena in nature and are typified by a relationship where the dependent variable is proportional to some power of the independent variable. That is:

\(y=ax^{k}\,\!\)

where \(a\) is the constant of proportionality and \(k\) is the scaling exponent. Transform by taking the logarithm of each side. Unlike the transformation of the exponential model, there is generally no mathematical advantage to using one logarithm over another. There is, however, a graphical advantage if the base-10 logarithm is used since MATLAB's semilogy, semilogx, and loglog plots use the base-10 logarithm.

\(\begin{align} y&=ax^k\\ \log_{10}\left(y\right)&=\log_{10}\left(ax^k\right)\\ \log_{10}\left(y\right)&=k*\left(\log_{10}\left(x\right)\right)^1+\log_{10}\left(a\right)*\left(\log_{10}\left(x\right)\right)^0\,\! \end{align}\)

meaning the transformed measurement \(\log_{10}\left(y\right)\) can be fit to a straight line function of \(\log_{10}\left(x\right)\). That is, using the transformations:

\(\begin{align} \xi&=\log_{10}\left(x\right) & \eta&=\log_{10}\left(y\right) \end{align}\)

The relationship between the slope and intercept for the transformed variables and the constants in the model equation are:

\(\begin{align} P(1)&=k & P(2)&=\log_{10}\left(a\right)\\ k&=P(1) & a&=10^{P(2)} \end{align}\)

Nonlinear Regression

Example Code

In the code below, note that the use of the variable fSSR was taken from Section 14.5 of the Chapra book.

%% Initialize the workspace
clear; format short e
figure(1); clf

%% Load and manipulate the data
load Cantilever.dat

Force = Cantilever(:,1) * 9.81;
Displ = Cantilever(:,2) * 2.54 / 100;

%% Rename and create model data
x = Force;
y = Displ;
xmodel = linspace(min(x), max(x), 100);

%% Define model equation and initial guesses
model = 'linear'
switch model
    case 'linear'
        yeqn = @(coefs, x) coefs(1)*x.^1 + coefs(2)*x.^0;
        InitGuess = [0 0]
    case 'trig with offset'
        yeqn = @(coefs, x) coefs(1)*cos(x*coefs(3)) + coefs(2)*sin(x*coefs(3)) + coefs(4);
        InitGuess = [0 0 100 0]
    otherwise
        error('Don''t know the model...')
end

%% Determine the function coefficients
fSSR = @(coefs, x, y) sum(( y - yeqn(coefs, x) ).^2)
[MyCoefs, Sr] = fminsearch(@(MyCoefsDummy) fSSR(MyCoefsDummy, x, y), InitGuess)

%% Generate estimates and model
yhat   = yeqn(MyCoefs, x);
ymodel = yeqn(MyCoefs, xmodel);

%% Calculate statistics
% Compute sum of the squares of the data residuals
St = sum(( y - mean(y) ).^2)

% Compute sum of the squares of the estimate residuals
Sr = sum(( y - yhat ).^2)

% Compute the coefficient of determination
r2 = (St - Sr) / St

%% Generate plots
plot(x,      y,      'k*',...
     x,      yhat,   'ko',...
     xmodel, ymodel, 'k-');
 xlabel('Independent Value')
 ylabel('Dependent Value')
 title('Dependent vs. Independent and Model')
 legend('Data', 'Estimates', 'Model', 0)

Questions

Post your questions by editing the discussion page of this article. Edit the page, then scroll to the bottom and add a question by putting in the characters *{{Q}}, followed by your question and finally your signature (with four tildes, i.e. ~~~~). Using the {{Q}} will automatically put the page in the category of pages with questions - other editors hoping to help out can then go to that category page to see where the questions are. See the page for Template:Q for details and examples.