Additional
Material

Fitting Parameters of Polynomials to Data

function StructureIdentificationExample(n, TrainingStart, TrainingEnd, ...
    TestStart, TestEnd, NoiseRange)
%Creates training and test data for the specified ranges:
%- One input X, one output Y.
%- Y = X^3 - 100X + 100 + noise
%- Uniformly distributed noise is added
%A polynomial modeling approach is executed for order n; parameters are fit
%to training data using the right matrix division command.
%Results are displayed graphically including model quality information (MSE).

% -- preparations -----------

if nargin<6, NoiseRange = 500, end;
if nargin<5, TestEnd = 25, end;
if nargin<4, TestStart = 6, end;
if nargin<3, TrainingEnd = 5, end;
if nargin<2, TrainingStart = -15, end;
if nargin<1, n = 5, end;

rand('state',12345);

X = TrainingStart:TrainingEnd;
for i=1:length(X)
    Y(i,1)  = X(i)^3 - 100*X(i) + 100;
    Y_(i,1) = rand(1)*NoiseRange - NoiseRange/2;
    Y(i,1)  = Y(i,1) + Y_(i,1);
end
Z = zeros(length(X), n+1);
for j = 1 : (n+1)
    for i=1:length(X)
        Z(i,j) = X(i)^(j-1);
    end
end

Y = Y'; Z = Z';

x = TrainingStart:0.01:TrainingEnd;
z = zeros(length(x), n+1);
for j = 1 : (n+1)
    for i=1:length(x)
        z(i,j) = x(i)^(j-1);
    end
end
z = z';

% -- training -----------

poly = Y/Z;
Yhat = poly*Z; error = Yhat-Y;
mse = sum(error.*error)/length(Y)

y = poly*z;

figure
plot(X,Y,'s', 'MarkerEdgeColor','k', 'MarkerFaceColor','k', 'MarkerSize',5)
hold on;
plot(x,y,'k', 'LineWidth',2);
hold off; grid on;
titlestr = ['Order ' num2str(n) ', MSE (training): ' num2str(mse)];
title(titlestr);

% -- test -----------

Xtest = TestStart:TestEnd;
for i=1:length(Xtest)
    Ytest(i,1)  = Xtest(i)*Xtest(i)*Xtest(i) - 100*Xtest(i) + 100;
    Ytest_(i,1) = rand(1)*NoiseRange - NoiseRange/2;
    Ytest(i,1)  = Ytest(i,1) + Ytest_(i,1);
end

Ztest = zeros(length(Xtest), n+1);
for j = 1 : (n+1)
    for i=1:length(Xtest)
        Ztest(i,j) = Xtest(i)^(j-1);
    end
end

Ytest = Ytest'; Ztest = Ztest';
YhatTest = poly*Ztest;
errorTest = YhatTest-Ytest;
mseTest = sum(errorTest.*errorTest)/length(YhatTest)

x = TrainingStart:0.01:TestEnd;
z = zeros(length(x), n+1);
for j = 1 : (n+1)
    for i=1:length(x)
        z(i,j) = x(i)^(j-1);
    end
end
z = z';
y = poly*z;

figure
plot([X Xtest],[Y Ytest],'s', 'MarkerEdgeColor','k', ...
    'MarkerFaceColor','k', 'MarkerSize',5)
hold on;
plot(x,y,'k', 'LineWidth',2);
hold off; grid on;
titlestr = ['Order ' num2str(n) ', MSE (test): ' num2str(mseTest)];
title(titlestr);