lsqd

Contents

sample data set

clear
close all
N = 20;
x = linspace(0,10,N);
a = 1.5;
b = 2.3;
sz = size(x);
y = a*x+b+0.2*randn(sz);
plot(x,y,'ko');
axis([0 10 0 20])
grid
save pset1 x y N

MATLAB regression

clear
close all
load pset1
c = polyfit(x,y,1);
fprintf('slope %g  intercept %g\n',c(1),c(2));
xf = [0 10];
yf = polyval(c,xf);
plot(x,y,'o',xf,yf,'k');
axis([0 10 0 20]);
grid
yfit = polyval(c,x);
v = var(y-yfit);
fprintf('rms fit: %g\n',sqrt(v));

slope 1.50013  intercept 2.2716
rms fit: 0.227891

MATLAB variance

Vyy1 = var(yfit-y);
fprintf('Vyy (using var) %g\n',Vyy1);

Vyy2 = norm(yfit-y)^2/(N-1);
fprintf('Vyy (using norm) %g\n',Vyy2);

sigma = sum((yfit-y).^2)/(N-1);
fprintf('Vyy (using summation) %g\n',sigma);

yrms = sqrt(sigma);
fill([0 0 10 10],[yf(1)-yrms,yf(1)+yrms,yf(2)+yrms,yf(2)-yrms],'y','EdgeColor','none');
grid
hold on;
plot(x,y,'o',xf,yf,'k');
axis([0 10 0 20]);
hold off

Vyy (using var) 0.0519343
Vyy (using norm) 0.0519343
Vyy (using summation) 0.0519343

MATLAB examples of simple linear fits

clear
load pset1

simple method

c = polyfit(x,y,1);
fprintf('slope %g intercept %g\n',c(1),c(2));

slope 1.50013 intercept 2.2716

using variances

Sxx = var(x);
Syy = var(y);
Sxy = (var(x+y)-var(x-y))/4;
m = Sxy/Sxx;
fprintf('(error in y) slope %g angle %g\n',m,atan(m)*180/pi);

(error in y) slope 1.50013 angle 56.3122

errors in x and y

theta = angle(Sxx-Syy+1j*2*Sxy)/2;
fprintf('(errors in x and y) slope %g angle %g\n', tan(theta),theta*180/pi);

(errors in x and y) slope 1.5026 angle 56.3558

Formulas for Minimum Deviation

Deviations in y

$$y-ybar = m (x-xbar)$$

$$b = ybar - m * xbar$$

b = mean(y) - m*mean(x);
fprintf('slope %g intercept %g\n',m,b);

d = y - m*x - b;
n = length(d)-1;
fprintf('std(d) = %g\n',std(d));
fprintf('sqrt(dot(d,d)/n) = %g\n',sqrt(dot(d,d)/n));
D = Sxx*Syy-Sxy*Sxy;
fprintf('sqrt(D/Sxx) = %g\n',sqrt(D/Sxx));

slope 1.50013 intercept 2.2716
std(d) = 0.227891
sqrt(dot(d,d)/n) = 0.227891
sqrt(D/Sxx) = 0.227891

Deviations in x and y

v = cos(theta)*(y-mean(y))-sin(theta)*(x-mean(x));
fprintf('std(v) = %g\n',std(v));
fprintf('sqrt(D/(Sxx+Syy)) = %g\n',sqrt(D/(Sxx+Syy)));

std(v) = 0.126332
sqrt(D/(Sxx+Syy)) = 0.1263

Vector norms

xc = x-mean(x);
v1 = norm(xc,2)^2/N;

fprintf('norm^2 %g\n', v1);

v2 = var(x);

fprintf('var %g\n',v2);

norm^2 9.21053
var 9.69529

“Least-Squares” via Vector Norms

see pslope.m 

pv = 1:8;
for p=pv
    m(p) = pslope(x,y,p);
end

fprintf('%10s%10s\n','p','m');
disp([pv' m']);

         p         m
    1.0000    1.4821
    2.0000    1.5001
    3.0000    1.5053
    4.0000    1.5088
    5.0000    1.5117
    6.0000    1.5140
    7.0000    1.5158
    8.0000    1.5171

Infinity norm

m = pslope(x,y,inf);

fprintf('slope (using infinity norm) %g\n',m);

xc = x - mean(x);
yc = y - mean(y);
Sxx = max(abs(xc))^2;

fprintf('Sxx %g\n',Sxx);

Sxy = (max(abs(xc+yc))^2-max(abs(xc-yc))^2)/4;

fprintf('Sxy %g\n',Sxy);

slope = Sxy/Sxx;

fprintf('infinity slope %g\n',slope);

slope (using infinity norm) 1.54003
Sxx 25
Sxy 38.5008
infinity slope 1.54003


Published with MATLAB® 7.14