% compute a linear regression to recover a linear model from noisy data
x = transpose(1:100);
y = 10 + 3 * x + 10 * zquant(rand(size(x)));
xx = [ones(size(x)) x];
b = inv(transpose(xx) * xx) * transpose(xx) * y;
b2 = xx \ y;
yp = xx * b;
e = sse(y,yp);
t = sst(y);
% r2 is coefficient of determination
r2 = (t-e)/t;
% plot data and yp model
figure(1);
plot(x,y,'o',x,yp);
figure(2);
plot(yp,y-yp,'o');

% compute a curvilinear regression for noisy data
x = transpose(1:.1:10);
y = 10 + (200 ./ x) + 8 * zquant(rand(size(x)));
xx = [ones(size(x)) (1./x)];
b = xx \ y;
yp = xx * b;
e = sse(y,yp);
t = sst(y);
% r2 is coefficient of determination
r2 = (t-e)/t;
figure(3);
plot(x,y,'o',x,yp);
figure(4);
plot(yp,y-yp,'o');

% compute a multiple linear regression for noisy data
x1 = transpose(1:100);
x2 = transpose((1:100) * 2 - 50);
y = 10 + 3 * x1 + 5 * x2 + 10 * zquant(rand(size(x1)));
xx = [ones(size(x1)) x1 x2];
b = xx \ y;
yp = xx * b;
e = sse(y,yp);
t = sst(y);
% r2 is coefficient of determination
r2 = (t-e)/t;
figure(5);
plot3(x1,x2,y,x1,x2,yp)
figure(6);
plot(yp,y-yp,'o');
