%% Robust Fitting
% This example shows how to compare the effects of excluding outliers and robust fitting. The example shows how to exclude outliers at an arbitrary distance greater than 1.5 standard deviations from the model. The steps then compare removing outliers with specifying a robust fit which gives lower weight to outliers.
% Create a baseline sinusoidal signal:

xdata = (0:0.1:2*pi)'; 
y0 = sin(xdata);

%
% Add noise to the signal with nonconstant variance.

% Response-dependent Gaussian noise
gnoise = y0.*randn(size(y0));

% Salt-and-pepper noise
spnoise = zeros(size(y0)); 
p = randperm(length(y0));
sppoints = p(1:round(length(p)/5));
spnoise(sppoints) = 5*sign(y0(sppoints));

ydata = y0 + gnoise + spnoise;

%
% Fit the noisy data with a baseline sinusoidal model, and specify 3 output arguments to get fitting information including residuals.

f = fittype('a*sin(b*x)'); 
[fit1,gof,fitinfo] = fit(xdata,ydata,f,'StartPoint',[1 1]);

%
% Examine the information in the fitinfo structure.

fitinfo

%
% Get the residuals from the fitinfo structure.

residuals = fitinfo.residuals;

%
% Identify "outliers" as points at an arbitrary distance greater than 1.5 standard deviations from the baseline model, and refit the data with the outliers excluded.

I = abs( residuals) > 1.5 * std( residuals );
outliers = excludedata(xdata,ydata,'indices',I);

fit2 = fit(xdata,ydata,f,'StartPoint',[1 1],...
           'Exclude',outliers);

%
% Compare the effect of excluding the outliers with the effect of giving them lower bisquare weight in a robust fit.

fit3 = fit(xdata,ydata,f,'StartPoint',[1 1],'Robust','on');

%
% Plot the data, the outliers, and the results of the fits. Specify an informative legend.

plot(fit1,'r-',xdata,ydata,'k.',outliers,'m*') 
hold on
plot(fit2,'c--')
plot(fit3,'b:')
xlim([0 2*pi])
legend( 'Data', 'Data excluded from second fit', 'Original fit',...
    'Fit with points excluded', 'Robust fit' )
hold off

%
% Plot the residuals for the two fits considering outliers:

figure 
plot(fit2,xdata,ydata,'co','residuals') 
hold on
plot(fit3,xdata,ydata,'bx','residuals')
hold off

% Copyright 2015 The MathWorks, Inc.