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 {\tt 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;
Disp  = Cantilever(:,2) * 2.54 / 100;

%% Rename and create model data
x = Force;
y = Disp;
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;
Disp  = Cantilever(:,2) * 2.54 / 100;

%% Rename and create model data
x = Force;
y = Disp;
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)

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;
Disp  = Cantilever(:,2) * 2.54 / 100;

%% Rename and create model data
x = Force;
y = Disp;
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.