
% Nonparametric regressions:
%
% For data (X_i,Y_i), find a function
%
%    Y_i = mu(X_i) + error_i
%
% First two ways:
%
% (1) Running-line smoother:
%    At each X_i, let mu(X_i) be the value obtained from
%    the linear regression of Y_i on X_i for a window of points
%    centered at X_i
%
% (2) Kernel smoother:
%    Set  
%       mu(X) = Sum(i=1,nn) Y_i K_h(X-X_i) / Sum(i=1,nn) K_h(X-X_i)
%
%    where
%        K_h(x) = (1/h) K_h(x/h)
%
%    is a kernel with ``bandwidth`` h, for example
%    the mean-zero variance-one Gaussian kernel
%

disp('Spending per high-school student and graduation rates');
disp('  for 50 states');
disp('Table 8.6 p380 in Hollander and Wolfe 1999');

xtab = ...
  [  1.  7971  65.5          26.  3998  79.6 ;
     2.  7151  62.3          27.  3943  85.4 ;
     3.  6564  77.4          28.  3919  69.1 ;
     4.  6230  84.9          29.  3858  77.3 ;
     5.  5471  74.4          30.  3840  65.9 ;
     6.  5329  69.8          31.  3794  76.3 ;
     7.  5207  78.7          32.  3786  74.0 ;
     8.  5201  74.1          33.  3744  61.1 ;
     9.  5051  88.3          34.  3691  71.9 ;
    10.  5017  71.7          35.  3623  75.8 ;
    11.  4989  78.4          36.  3608  65.3 ;
    12.  4789  73.0          37.  3519  88.3 ;
    13.  4747  84.9          38.  3434  61.0 ;
    14.  4692  73.6          39.  3408  64.6 ;
    15.  4462  74.7          40.  3368  66.7 ;
    16.  4457  74.1          41.  3249  79.6 ;
    17.  4386  90.9          42.  3138  61.4 ;
    18.  4369  75.6          43.  3093  71.7 ;
    19.  4246  74.4          44.  3068  69.3 ;
    20.  4246  87.3          45.  3011  69.0 ;
    21.  4164  77.1          46.  2989  77.2 ;
    22.  4149  71.6          47.  2718  74.9 ;
    23.  4124  85.8          48.  2667  75.4 ;
    24.  4092  58.0          49.  2548  66.9 ;
    25.  4076  80.2          50.  2454  79.4  ];

disp(' ');
disp('Data:');
fprintf ('  %3d.   %4d   %4.1f     %3d.   %4d   %4.1f\n', xtab'); %'

xx = [ xtab(:,2);  xtab(:,5) ];
yy = [ xtab(:,3);  xtab(:,6) ];

disp(' ');
disp('Data sorted in increasing order by X');
ytab = sortrows([ xx  yy]);
xx = ytab(:,1);  yy=ytab(:,2);

nn=size(ytab,1);
xmean=sum(xx)/nn;  ymean=sum(yy)/nn;
disp(' ');
fprintf('n=%d  X-Y pairs\n',nn);
xmin=min(xx);  xmax=max(xx);
ymin=min(yy);  ymax=max(yy);
fprintf('  Xmin=%d   Xmax=%d   Ymin=%4.1f   Ymax=%4.1f\n',...
  xmin,xmax, ymin,ymax);
fprintf('  Xmean=%6.1f   Ymean=%5.2f\n', xmean,ymean);

% For linear smoother
dw=9;  hw=floor(dw/2);

disp(' ');
disp('Linear smoother with');
fprintf('   Window size in X-values: %d   Half-width: %d\n', dw,hw);

yls = zeros(nn,1);

% Loop over (X_i,Y_i) pairs
% For each (most) points, make a regression line using the hw points
%   to the right and the hw points to the left, for a total of
%   dw=2*hw+1 points.
% Use the Y values at endpoints
% Otherwise, if there are fewer than hw points to the right or left,
%   replace hw by the largest possible positive value
nw=0;  dwsum=0;   % For calculating the average window width
for i=1:nn;
  % Find minimum distance to i=1 or i=nn
  n1=i-1;  n2=nn-i;  hm=min([n1 n2]);
  % Use the Y values at endpoints and go to the next i value
  if hm==0  yls(i)=yy(i);  continue;
  % Use hw normally but not close to the endpoints
  elseif hm>hw hm=hw;  end;
  % Find the regression line determined by these points
  nm = 2*hm + 1;
  xxr = xx(i-hm:i+hm);
  yyr = yy(i-hm:i+hm);
  xm=sum(xxr)/nm;  ym=sum(yyr)/nm;
  b1=0;  b2=0;
  for j=1:nm;  b1=b1+(yy(j)-ym)*(xx(j)-xm);
    b2=b2+(xx(j)-xm)*(xx(j)-xm);  end;
  beta=b1/b2;  mu=ym-beta*xm;
  % Apply at that X value
  yreg = mu + beta*xx(i);
  % To be safe (not really needed)
  if yreg > ymax  yreg=ymax;
  elseif yreg < ymin yreg=ymin;  end;
  % Store linearly-smoothed value
  yls(i) = yreg;
  % Keep track of average window width
  nw=nw+1;  dwsum=dwsum+(max(xxr)-min(xxr));
  end;

avhw=dwsum/(2*nw);
fprintf('   Average window width: %6.1f  Half-width: %6.1f\n', 2*avhw,avhw);
hh=2*avhw;

disp(' ');
fprintf('Kernel smoother for Gaussian kernel with stddev=h=%g:\n', hh);
disp('(Twice the average half-width.)');

% Find kernel smoother
yks = zeros(nn,1);
for i=1:nn;
  xsum=0;  ysum=0;  xc=xx(i);
  for j=1:nn;
    xnorm = (xx(j)-xc)/hh;
    xwt = exp(-xnorm*xnorm/2);
    ysum=ysum+yy(j)*xwt;
    xsum=xsum+xwt;
    end;
  if xsum>0 yks(i)=ysum/xsum;
  else yks(i)=ymean;  end;
  end;

disp(' ');
disp('Data is blue');
disp('Linear smoother is green (based on 9-point window)');
fprintf('Kernel smoother is red (based on h=%g)\n',hh);

plot(xx,yy,xx,yls,xx,yks);


