
% Kernel density estimation:
%
% For a sample X_i and a kernel function K(x) such that
%
%     K(x)=K(-x)     Int(-infty,infty) K(y)dy = 1
%
%        Int(-infty,infty) y^2 K(y) dy = 1
%
% estimate the density f(x)=f_X(x) by
%
%        f(x) = (1/n) Sum(i=1,n) (1/h)K((x-X_i)/h)
%
% where h>0 is the ``bandwidth``
%
% There are theorems stating that the ``minimal risk`` setting
% of the bandwidth  h  is
%
%        h = C/n^{1/5}
%
% In the following, we use  K(x) = (1/sqrt{2pi}) exp(-x^2/2)
%

disp('Estimating the density of X_i for the dataset:');
disp('Spending per high-school student (X_i) and graduation rates (Y_i)');
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 = sort([ xtab(:,2);  xtab(:,5) ]);
nn = size(xx,1);

fprintf('\nX-values (n=%d):\n', nn);
r=0;  k=0;  n2 = ceil(nn/10);
for i=1:n2;
  for j=1:10;
    r=r+1;  if r>nn  break;  end;
    fprintf('   %4d', xx(r));  end;
  if j>1 fprintf('\n');  end;
  if r>nn  break;  end;
  end;

cc=1/sqrt(2*pi);    % Normalizing constant for kernel

harr = [ 10  50  200  500 ];
numk=length(harr);

maxx=10000;

yyk = zeros(numk,maxx);
xxk = 1:maxx;

disp(' ');
disp(' ');
disp('Estimating the density f_X(x) for bandwidths:');
for a=1:numk;
  hh=harr(a);
  fprintf('   h=%g\n', hh);
  for i=1:maxx;
    x=i;  fsum=0;
    for j=1:nn;
      xn=(xx(j)-x)/hh;  fsum=fsum+exp(-xn*xn/2);  end;
    yyk(a,i) = cc*fsum/(nn*hh);
    end;
  end;

fprintf('\nApproximating Int(0,%g) by Sum(1,%g):\n', maxx,maxx);
fprintf('  Int = Int(0,%g) f(x)dx    Intsq = Int(0,%g) f(x)^2 dx\n',...
  maxx,maxx);
fprintf('  Intsqsc is after scaling (0,%g) --> (0,1)\n',maxx);
disp(' ');

% For an mxn matrix A, sum(A,2) is an m-column-vector of n-row-sums
% yyk is numk x maxx
intval=sum(yyk,2);
intsqval = sum(yyk.*yyk,2);

itab = [ harr'  intval  intsqval  intsqval*maxx];  %'
fprintf ('  h=%3d    Int=%6.4f   Intsq=%10.8f    Intsqsc=%6.4f\n', itab');  %'

yy1=yyk(1,:);  yy2=yyk(2,:);
yy3=yyk(3,:);  yy4=yyk(4,:);

disp(' ');
%plot(xxk,yy1,xxk,yy2,xxk,yy3,xxk,yy4);


% Axes
aax=[0 maxx 0 0.0008];
%% harr = [ 10  50  200  500 ];

subplot(2,2,1);
plot(xxk,yy1);
title('h = 10');
axis(aax);

subplot(2,2,2);
plot(xxk,yy2);
title('h = 50');
axis(aax);

subplot(2,2,3);
plot(xxk,yy3);
title('h = 200');
axis(aax);

subplot(2,2,4);
plot(xxk,yy4);
title('h = 500');
axis(aax);







