 AMath 584, Autumn Quarter, 2011 Applied Linear Algebra and Introductory Numerical Analysis

Least squares¶

In this example exact data points on a parabola are generated and then normally-distributed (Gaussian) random noise is added. You can change the amount of noise by adjusting the parameter delta.

The least squares solution is then computed and plotted along with the data.

This simulates an experiment in which is ball is thrown upwards and released at y = y0 with velocity v0. The height as a function of time should then be y(t) = y0 + v0*t - 0.5*g*t^2. From noisy measurements we might want to estimate g, the third component of the coefficient vector c below. Python version:¶

[lstsq1.py]

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 from pylab import * # Generate synthetic data: y0 = 1. v0 = 20. g = 9.81 t = linspace(0,4,20) m = len(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) y = y + noise # Plot points figure(1) clf() plot(t,y,'o') title('Data points') ans = raw_input("Hit return to fit... ") # Create matrix A: A = ones((m,3)) # fix second and third columns: A[:,1] = t A[:,2] = -0.5 * t**2 # Solve using lstsq: output = lstsq(A,y) c = output # see documentation for other elements of output print "\nCoefficients estimated with lstsq: (y0, v0, g):" print c # Solve using QR directly: Q,R = qr(A) yhat = dot(Q.T, y) c = solve(R, yhat) print "\nCoefficients estimated with QR: (y0, v0, g):" print c # Plot the function we found with lots of points: tt = linspace(0, 4, 1000) yy = c + c*tt - 0.5*c*tt**2 plot(tt, yy, 'r') title('Data and least squares fit')

Matlab version:¶

[lstsq1.m]

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 % 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')