% Example 7.1
% This shows the univariate predictor loess procedure.

load environmental
% This has all four variables. We will
% use the ozone as the response and
% temperature as the predictor.
% First do a scatterplot
plot(Temperature,Ozone,'.')
xlabel('Temperature (Fahrenheit)')
ylabel('Ozone (PPB)')

n = length(Temperature); % Find the number of data points
x0 = 78; % Find the estimate at this point
alpha = 2/3;
lambda = 1;		
k = floor(alpha*n);
% First step is to get the neighborhood. 
dist = abs(x0 - Temperature);
[sdist,ind] = sort(dist);
% Get the points in the neighborhood.
Nx = Temperature(ind(1:k));							
Ny = Ozone(ind(1:k));
delxo = sdist(k);  % Maximum distance of neighborhood

% Delete the ones outside the neighborhood.
sdist((k+1):n) = []; 								
% These are the arguments to the weight function.
u = sdist/delxo;
% Get the weights for all points in the neighborhood.
w = (1 - u.^3).^3;

% Now using only those points in the neighborhood,
% do a weighted least squares fit of degree lambda.
% We will follow the procedure in 'polyfit'.
x = Nx(:); y = Ny(:); w = w(:);
% Get the weight matrix
W = diag(w);
% Get the right matrix form for the X: 1, x, x^2.
A = vander(x);
A(:,1:length(x)-lambda-1) = [];
V = A'*W*A;
Y = A'*W*y;
[Q,R] = qr(V,0); 
p = R\(Q'*Y); 
% The following is to fit Matlab's convention
% for polynomials
p = p';				

% This is the polynomial model for the local fit.
% To get the value at that point, use polyval.
yhat0 = polyval(p,x0);

% Now call the loess function and plot the result.
% Get a domain over which to evaluate the curve.
X0 = linspace(min(Temperature),max(Temperature),50);
yhat = loess(Temperature,Ozone,X0,alpha,lambda);
% Plot the results.
plot(Temperature,Ozone,'.',X0,yhat)
xlabel('Temperature (Fahrenheit)'),ylabel('Ozone (PPB)')
