clear all;
clf;
% First, generate some random data for our "least squares fit"
m_true = 3;
b_true = 1;
N = 21;
x = linspace(0,1,N);
y = m_true*x(:) + b_true + 1*rand(N,1);
% Plot our "measurements"
plot(x,y,'.','markersize',20);
% ------------------------------------------
% Find the best fit to our model equations
% y = mx + b
% Each of the four methods below solves the
% normal equations
% (A'A)*xhat = A'*R;
% but in different ways, some which are more
% numerically more stable than others.
% ------------------------------------------
% Methods 1, 2 and 3 require that we set up the
% linear system explicitly
A = [x(:) ones(N,1)];
R = y;
% Method 1 : (not recommended)
% Solve normal equations directly
AtA = A'*A;
xhat1 = AtA\(A'*R);
% Method 2 :
% Use the 'backslash operator'. This will
% solve the normal equations "behind the scenes"
xhat2 = A\R;
% Method 3 :
% Functional form of A\R
xhat3 = mldivide(A,R);
% Method 4 : Use 'polyfit'
% Note that this approach does not require that
% we construct the matrix A explicitly :
deg = 1; % Degree of polynomial we want to fit :
xhat4 = polyfit(x(:),y,deg);
xhat4 = xhat4(:);
% Each of the four methods above will produce the same results
% although Method 1 is not recommended. Method 4 is probably
% the most convenient, since it doesn't require that we
% form the linear system explicitly.
format long
disp([xhat1 xhat2 xhat3 xhat4]);
xhat = xhat4; % Choose one.
% ------------------------------------------------
% Plot the best fit line
% ------------------------------------------------
% First we need to get a set of points at which to
% evaluate our model equations :
xf = linspace(0,1,100);
% Method 1 : (not recommended)
% To illustrate what the vector 'xhat' is, we
% show how it can be used to evaluate the model
% equations.
m = xhat(1);
b = xhat(2);
yf = m*xf + b; % (not recommended)
% Method 2 : (recommended)
% Pass the coefficients in xhat directly to 'polyval' to
% evaluate the polynomial :
yf = polyval(xhat,xf);
% Both methods produce essentially the same results, but method
% 2 is much more stable numerically.
% Plot the results
hold on;
plot(xf,yf,'r','linewidth',2);
shg;