1
$\begingroup$

Below is my own approach to implement the Least Squares Regression algorithm in MATLAB. Could you please take a look and tell me if it makes sense; if it does exactly what is supposed to do?

EDIT: Please, pay attention to the commented commands as well.

function [MSE_train_mean,MSE_train_std,MSE_test_mean,MSE_test_std,w_star] = perform_lsr(X,Y,Z)

% 604 days for the training set (29/05/2009 - 18/10/2011)
% 145 days for the evaluation set (19/10/2011 - 16/05/2012)
Xtran = X(1:end-145,:);
Ytran = Y(1:end-145,:);
Ztran = Z(1:end-145,:);
Xeval = X(end-144:end,:);
Yeval = Y(end-144:end,:);
Zeval = Z(end-144:end,:);

% Set number of runs.
runs = 20;

% Set number of folds.
folds = 5;

% Initialize auxiliary matrices.
MSE_train = zeros(runs,folds);
MSE_test = zeros(runs,folds);
w_star = zeros(runs,folds,size(Xtran,2));

% Perform runs.
for r=1:runs,

% Split original dataset into training and test set.
split_assignments = cross_val_split(folds,size(Xtran,1));

% Perform folds.
for f=1:folds,

% Assign explanatory variables (x).
x_train = Xtran((split_assignments(:,f)==0),:);
x_test = Xtran((split_assignments(:,f)==1),:);

% Assign respond variable (y).
y_train = Ytran((split_assignments(:,f)==0),:);
y_test = Ytran((split_assignments(:,f)==1),:);

% Retrieve size of matrices.
[l_train_n,l_train_m] = size(x_train);
[l_test_n,l_test_m] = size(x_test);

% Estimate parameters (w0,w1,...) of the 1st order model from training set.
w_star(r,f,:) = (x_train' * x_train) \ (x_train' * y_train);
w_star_temp = w_star(r,f,:);
w_star_temp = w_star_temp(:);

% Apply the learned weights on both training and test sets and compute the
% corresponding MSEs.
MSE_train(r,f) = (1 / l_train_n) * (w_star_temp' * (x_train') * x_train * w_star_temp - 2 * y_train' * x_train * w_star_temp + y_train' * y_train);
MSE_test(r,f) = (1 / l_test_n) * (w_star_temp' * (x_test') * x_test * w_star_temp - 2 * y_test' * x_test * w_star_temp + y_test' * y_test);

end
end

% % Plot both training and test sets' MSEs as a function of runs and folds.
% figure, mesh(MSE_train), title('Training Set`s MSEs vs Runs and Folds');
% figure, mesh(MSE_test), title('Testing Set`s MSEs vs Runs and Folds');
% 
% % Plot both training and test sets' mean of MSEs accompanied by their
% % corresponding std of MSEs.
% figure, errorbar(1:runs,mean(MSE_train,2),std(MSE_train,0,2));
% figure, errorbar(1:runs,mean(MSE_test,2),std(MSE_test,0,2));

% Average over folds and then over runs.
MSE_train_mean = mean(mean(MSE_train,2),1);
MSE_train_std = std(mean(MSE_train,2),0,1);
MSE_test_mean = mean(mean(MSE_test,2),1);
MSE_test_std = std(mean(MSE_test,2),0,1);
w_star_temp = mean(mean(w_star,2),1);
w_star = w_star_temp(:);

% % Calculate in-sample residuals.
% y_hat_tran = Xtran * w_star;
% e_tran = Ytran - y_hat_tran;
% % Check that they sum up to zero.
% fprintf('In-sample residuals sum up to %i.\n',sum(e_tran));
% % Finally, plot them.
% figure, scatter(1:size(e_tran,1),e_tran);
% 
% % Calculate out-of-sample residuals.
% y_hat_eval = Xeval * w_star;
% e_eval = Yeval - y_hat_eval;
% % Check that they sum up to zero.
% fprintf('Out-of-sample residuals sum up to %i.\n',sum(e_eval));
% % Finally, plot them.
% figure, scatter(1:size(e_eval,1),e_eval);

% Print summary statistics.
fprintf('Over %i runs, %i folds, and regarding the training set, the mean and standard deviation is %.10f and %.10f respectively.\n',runs,folds,MSE_train_mean,MSE_train_std);
fprintf('Over %i runs, %i folds, and regarding the testing set, the mean and standard deviation is %.10f and %.10f respectively.\n',runs,folds,MSE_test_mean,MSE_test_std);

end
$\endgroup$
  • 3
    $\begingroup$ I think this is too general a question. Stackexchange isn't meant to be a place where people check other people's codes but a place where one can ask concrete questions. So if there is a specific part of the algorithm that you have questions on, then this here would be the right place to ask. If all you want is apply the LS algorithm to a data set I suggest you simply apply the built-in commands provided by many different software packages. $\endgroup$ – Wolfgang Bangerth Aug 11 '12 at 11:29
2
$\begingroup$

Chapter 44, "Approximation of Discrete Data", in Advanced Engineering Mathematics by Robert J. Lopez gives the following algorithm for least squares regression:

clear
d = 6;
x = [1, 2, 3, 4, 5, 6, 7]
y = [3, 4, 1, 3, 2, 4, 9]
p = [0:1:d]
for i=1:length(y)
    for j=1:length(p)
        A1(i,j)=x(i)^p(j);
    end
end
A1
A2 = A1'
B0=y
for k=1:length(y)-1
    B0=cat(1,B0,y);
end
B0
B1=B0'
A3 = A2*B1
C1 = A2*A1
C2 = inv(C1)
C3 = C2*A3
C4 = C3'
pol = C4(1,1:length(p))
C5 = A1*C4
yreg = sum(C5,2)

which I have written in Scilab. Should work in Matlab since Scilab is a Matlab clone. "pol" is the coefficients of the polynomial of degree "d" that best fits the data.

Using Mathematica:

(*"d" is the degree of the polynomial to be fitted to x,y pairs*)
d = 6;
x = {1, 2, 3, 4, 5, 6, 7};
y = {3, 4, 1, 3, 2, 4, 9};
p = Range[d + 1] - 1;
A1 = Table[Table[x[[i]]^p[[j]], {j, 1, Length[p]}], {i, 1, Length[y]}];
A2 = Transpose[A1];
B1 = Transpose[Table[y, {i, 1, Length[y]}]];
A3 = A2.B1;
C1 = A2.A1;
C2 = Inverse[C1];
C3 = C2.A3;
C4 = Transpose[C3];
Print["polynomial"]
pol = Total[C4[[1, All]]*(z^p)]
g1 = ListPlot[
   Prepend[Table[{x[[i]], y[[i]]}, {i, 1, Length[x]}], {0, 0}], 
   PlotMarkers -> Automatic];
g2 = Plot[pol, {z, 0, Length[x] + 1}];
Show[g1, g2, PlotRange -> {Min[y] - 1, Max[y] + 1}, ImageSize -> Full]
C5 = A1*C4;
yreg = Total[Transpose[C5]];

to plot the polynomials for different values of the polynomial degree parameter "d" we get:

d=0

polynomial of degree 0

= 26/7 = 3.71429

degree zero polynomial or average of y values of data

d=1

polynomial of degree 1

= 1 + (19 z)/28

= 1. + 0.678571 z

polynomial of degree 1

d=2

polynomial of degree 2

= 46/7 - (85 z)/28 + (13 z^2)/28

= 6.57143 - 3.03571 z + 0.464286 z^2

polynomial of degree 2

d=3

polynomial of degree 3

= 11/7 + (335 z)/126 - (101 z^2)/84 + (5 z^3)/36

= 1.57143 + 2.65873 z - 1.20238 z^2 + 0.138889 z^3

polynomial of degree 3

d=4

polynomial of degree 4

= 2 + (2743 z)/1386 - (7 z^2)/8 + (31 z^3)/396 + z^4/264

= 2. + 1.97908 z - 0.875 z^2 + 0.0782828 z^3 + 0.00378788 z^4

polynomial of degree 4

d=5

polynomial of degree 5

= -(107/7) + (7947 z)/220 - (571 z^2)/24 + (3631 z^3)/528 - (241 z^4)/264 + (11 z^5)/240

= -15.2857 + 36.1227 z - 23.7917 z^2 + 6.87689 z^3 - 0.912879 z^4 + 0.0458333 z^5

polynomial of degree 5

d=6

polynomial of degree 6

= -110 + (5093 z)/20 - (4179 z^2)/20 + (3965 z^3)/48 - (815 z^4)/48 + (419 z^5)/240 - (17 z^6)/240

= -110. + 254.65 z - 208.95 z^2 + 82.6042 z^3 - 16.9792 z^4 + 1.74583 z^5 - 0.0708333 z^6

polynomial of degree 6

Nevermind the dummy point {0,0} in origo. It is just there to make the plots start at zero which I did not know how to do in a better way. {x,y} = {0,0} is not part of the data used in the least squares regression.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Explicit formulation and use of the normal equations. Rather dodgy, to put it mildly. Hey look Ma, I just squared the condition number. $\endgroup$ – Mark L. Stone Sep 13 '15 at 17:32
  • $\begingroup$ Well as I said, this is how Lopez does it. Thanks for the comment anyways! $\endgroup$ – Mats Granvik Sep 14 '15 at 13:07
  • $\begingroup$ Lopez is just one among many authors who apparently don't know what they 're doing in numerical methods/analysis. $\endgroup$ – Mark L. Stone Sep 14 '15 at 13:13

Not the answer you're looking for? Browse other questions tagged or ask your own question.