
% Generate synthetic data:
y0 = 1.;
v0 = 20.;
g = 9.81;
t = linspace(0,4,20)';
m = length(t);

% Exact values:
y = y0 + v0*t - 0.5*g*t.^2;

% Add normally distributed noise of size delta:
delta = 0.3;
noise = delta * randn(m,1);
y = y + noise;

% Plot points
figure(1);
clf
plot(t,y,'o')
title('Data points')

ans = input('Hit return to fit... ')


% Create matrix A:
A = ones(m,3);
% fix second and third columns:
A(:,2) = t;
A(:,3) = -0.5 * t.^2;

% Solve using backslash:
c = A \ y;

disp(['Coefficients estimated with lstsq: (y0, v0, g):'])
c


% Solve using QR directly:
[Q,R] = qr(A);
yhat = Q' * y;
c = R \ yhat;

disp(['Coefficients estimated with QR: (y0, v0, g):'])
c

% Plot the function we found with lots of points:
hold on
tt = linspace(0, 4, 1000);
yy = c(1) + c(2)*tt - 0.5*c(3)*tt.^2;
plot(tt, yy, 'r')
title('Data and least squares fit')
