function least_abs_dev
% Curve Fitting Demonstration comparing least squares to least absolute
% deviation in the case of non-gaussian noise.
hold off

slope=4*(rand(1)-0.5)
yintercept=2*(rand(1)-0.5)

%Make some data, complete with 'noise'
global datax datay %share these with the badness function (below)
datax = -1:0.1:1;
datay = (slope*datax)+yintercept;

figure(1)
hold off
plot(datax,datay,'r--'); %Plot the data *without* noise
hold on %save current display for overplotting

[M,N] = size(datax);
error = 0.5*(rand(M,N)-0.5).*(max(datay)-min(datay)); %create noise
outliers = rand(M,N) > 0.9 %select outliers
error(outliers) = 50*error(outliers);
datay = datay + error; %add noise to data

plot(datax,datay,'ro'); %Overlay more realistic noised data

%********* least-squares *********************************
%Find least-squares best fit slope and yintercept from noised data
result = amoeba('badness',[0,0],40);
slope=result(1)
yintercept=result(2)
result = amoeba('badness',result,40);
slope=result(1)
yintercept=result(2)

%Generate least-squares best fit curve, plot for comparison
fit_y = (slope*datax)+yintercept;
plot(datax,fit_y,'b-'); %Plot the fit

%********* least-abs-deviation *********************************
%Find least-abs-deviation best fit slope and yintercept from noised data
result = amoeba('badness2',[0,0],40);
slope=result(1)
yintercept=result(2)
result = amoeba('badness2',result,40);
slope=result(1)
yintercept=result(2)

%Generate least-squares best fit curve, plot for comparison
fit_y = (slope*datax)+yintercept;
plot(datax,fit_y,'k-'); %Plot the fit
legend('true', 'data', 'least squares', 'least abs dev')
print('least_abs_dev.pdf','-dpdf')

%end of main program




