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$

closed as too localized by Arnold Neumaier, David Ketcheson, Geoff Oxberry Aug 21 '12 at 20:08

This question is unlikely to help any future visitors; it is only relevant to a small geographic area, a specific moment in time, or an extraordinarily narrow situation that is not generally applicable to the worldwide audience of the internet. For help making this question more broadly applicable, visit the help center. If this question can be reworded to fit the rules in the help center, please edit the question.

  • 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.

$\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.