%
% 
% epage 432
% 
% Written by:
% -- 
% John L. Weatherwax                2008-02-20
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

clear all; 
close all; 
clc; 

addpath('../../Code/CSTool');
load environ; 

[swind,inds] = sort(wind); 
sozone = ozone(inds);
x=swind; y=sozone; 

fd=figure; rd=plot( wind, ozone, '.k' ); hold on; 
title( 'the environ data with a polynomial model' ); 

% lets assume a simple polynomial model ... maybe a cubic ... 
[pcoef,s] = polyfit(x,y,3); 
  
% evaluate this polynomial model:
yp = polyval(pcoef,x); 
  
figure(fd); pm=plot(x,yp,'-xr' ); 
xlabel( 'independent variable x' ); ylabel( 'y response' ); 

% compute the residuals;
res_eps = y - yp;

%far=figure; plot( x, abs(res_eps), 'k.' ); hold on; 
%title( 'x v.s. abs(residual)' ); 
far=figure; plot( x, res_eps, 'k.' ); hold on; 
title( 'x v.s. residual' ); 
xlabel( 'independent variable x' ); ylabel( 'y response' ); 

% smooth the "function" abs(residual) using loess: 
x_eval = x;
alpha  = 1/3;      % <- the fraction of the data to use for {\em local} fits ... 
deg    = 1;        % <- the degree of the local fit 
%y_loess = csloess(x,abs(res_eps),x_eval,alpha,deg);
y_loess = csloess(x,res_eps,x_eval,alpha,deg);

figure(far); plot( x, y_loess, 'rx-' ); 

% add in this loess correction to the main model: 
figure(fd); lc=plot(x,yp+y_loess, 'bo-' ); axis([2,21,0,120]); grid on; 
legend( [rd,pm,lc], {'raw environ data','polynomial model','loess corrected'} ); 
saveas( gcf, ['../../WriteUp/Graphics/Chapter10/prob_10_4_loess_correction'], 'epsc' ); 

clear functions;