function ShowFit()
% Fitting data that looks linear two different ways.
% Set up the plot window...
close all
s = 5;
plot([-s s],[0 0],'--k',[0 0],[-s s],'--k')
z = linspace(-s,s);
axis equal manual
hold on
shg
% Generate and display linear data with added noise...
n = 15;
m0 = .6;
b0 =  -2;
x = linspace(-s+1,s-1,n)';
delta = .1;
noise = delta*randn(n,1);
y = (m0*x + b0) + noise;
plot(x,y,'ok',x,y,'*k')
xlabel('    The Data   ','fontsize',16)
pause(2)
% Fitting data in the L2 (or least squares) sense.
% Choose m abd b to minimize sum((y - (m*x+b)).^2), the sum of
% the squares of the discrepancies.
[m,b] = L2(x,y);
plot(z,m*z+b,'b','linewidth',2)
xlabel('     L2  Data Fitting     ','fontsize',16)
pause(2)
% Add an outlier data point...
xlabel('    Add an Outlier   ','fontsize',16)
xi = 4; yi = 4;
plot(xi,yi,'*m','markersize',10)
xlabel('   Recompute the L2 fit    ','fontsize',16)
pause(2)
% Compute the new L2 fit
[mOut,bOut] = L2([x;xi],[y;yi]);
plot(z,mOut*z+bOut,'m','linewidth',2)
xlabel('            L2 Data Fitting                ','fontsize',16)
title(sprintf('m = %5.3f  b = %5.3f        mOut = %5.3f  bOut = %5.3f',m,b,mOut,bOut),'fontsize',16) 
pause(2)

%%%  Fitting dat in the L1 sense.
% Repeat only now we choose m and b to minimize
% sum(abs(y - (m*x+b))), the sum of the absolute
% values of the discrepancies. 
figure
plot([-s s],[0 0],'--k',[0 0],[-s s],'--k')
axis equal manual
hold on
shg
% Generate and display linear data with added noise...
plot(x,y,'ok',x,y,'*k')
[m,b] = L1(x,y);
plot(z,m*z+b,'b','linewidth',2)
plot(xi,yi,'*m','markersize',10)
[mOut,bOut] = L1([x;xi],[y;yi]);
plot(z,mOut*z+bOut,'m','linewidth',2)
xlabel('            L1 Data Fitting                ','fontsize',16)
title(sprintf('m = %5.3f  b = %5.3f        mOut = %5.3f  bOut = %5.3f',m,b,mOut,bOut),'fontsize',16) 


function [m,b] = L2(x,y)
% x and y are column n-vectors. Fits a line y=mx+b
% to the data (x(1),y(1)),...,(x(n),y(n)) in the L2 sense.
% See FVL p.382-3 
n = length(x);
xbar = mean(x);
ybar = mean(y);
sigxx = sum((x-xbar).^2)/n;
sigxy = sum((x-xbar).*(y-ybar))/n;
m = sigxy/sigxx;
b = ybar-m*xbar;

function [m,b] = L1(x,y)
% x and y are column n-vectors. Fits a line y=mx+b
% to the data (x(1),y(1)),...,(x(n),y(n)) in the L1 sense.
% Do not worry about the implementation--beyond the course!
n = length(x);
A = [ones(n,1) x];
b = y;
z0 = A\b;
z = fminsearch(@(x) phi(x,A,b),z0);
b = z(1);
m = z(2);
  
function alfa = phi(x,A,b)
alfa = sum(abs(A*x-b));
