% Plot of some biweight functions
deltav = -2:0.05:2;
C1 = deltav.^2;  % Squared error
C2 = abs(deltav);  % L1
b = 0.5;  C3 = b^2*log(1+deltav.^2/b^2);  % Cauchy
fh = figure;
plot(deltav,C1,'-',deltav,C2,'--',deltav,C3,'-.');
legend('Squared error','L1','Cauchy');
exportfig(gcf,'robustlsqline_fig02.eps','width',10,'height',2.5,'fontmode','fixed', 'fontsize',8,'Color','cmyk');

% y = m*x + c
np = 15;  m0 = 0.3;  c0 = 1;
nstddev = 0.2;
x0 = 0:10;  y0 = m0*x0 + c0;  % actual line

% Observed data
x = 10*rand(1,np);
y = m0*x + c0 + nstddev*randn(size(x));

% Make last point outlier
x(end) = 9.5;  y(end) = 0.5;

% Plot points (the last point is the outlier)
fh = figure;
plot(x0,y0,'b--');
axis([0 10 0 5]);
hold on;  sh = scatter(x(1:end-1),y(1:end-1),'r+');  hold off;
hv = [];
firstpass = 1;
while 1
  if ~isempty(hv), delete(hv); hv = []; end
  hold on;  sh = scatter(x(end),y(end),'go');  hold off;  hv = [hv; sh];
  
  % Best fit least squares
  A = [x' ones(size(x'))];
  theta = inv(A'*A)*A'*y';
  m = theta(1);  c = theta(2);
  
  % Draw best fit line
  x0 = 0:10;
  y0 = m*x0 + c;
  hold on;  ph = plot(x0,y0,'g-');  hold off;  hv = [hv; ph];
  title('Least squares fit line with outlier');

  if firstpass
    exportfig(gcf,'robustlsqline_fig01.eps','width',10,'height',6.5,'fontmode','fixed', 'fontsize',8,'Color','cmyk');
    firstpass = 0;
  end    
  
  % Get next point
  [xp,yp,b] = ginput(1);
  if b==3, break; end
  x(end) = xp;  y(end) = yp;
end

