function z=kdsmpl1d(A,k,inc,adapt,ss,h);

% KDSMPL1D    Compare sample and population K.d.e.'s.
%
%             KDSMPL1D(A,K,INC,ADAPT,SS,H) plots two univariate Kernel 
%             Density Estimates on the same axes:
%
%             1: The 'population' K.d.e. using the whole data set.
%             2: The 'sample' K.d.e. using a sample of size SS.
%
%             The process is repeated for nine samples (without replacement)
%             of size SS.  The kernel in use is given by the function in K.
%             'KNORM' is the multivariate normal distribution.
%             'KEPAN' is the Epanechnikov kernel.
%
%             The data is stored in the rows of A.
%
%             INC is the number of x increments to use 
%             in plotting the density function.  As a guide, for a 
%             rough, quick plot use, INC=10, for a smooth detailed 
%             plot use INC=100.
%
%             If ADAPT<>0 then the Adaptive Kernel method is used.
%
%             The optional parameter H represents the window width.
%             If omitted, a simplistic 'optimal' value is used.
%
%             Christian C. Beardah 1993

n=length(A);    % Establish dimensions of A,
                % n=number of data points,
               
x=A(:,1);       % Extract the first column of A into x.
xmin=min(x);    % Find the minimum value of x.
xmax=max(x);    % Find the maximum value of x.
xrange=xmax-xmin; % Find the range of x.

% xa holds the x 'axis' vector, defining a grid of x values where 
% the k.d. function will be evaluated and plotted.
% xinc is the actual increment used - needed for the numerical integration routine.

xa=linspace(xmin-xrange/8,xmax+xrange/8,inc); 

xinc=1.25*xrange/(inc-1);

% Calculate kernel density estimate values corresponding to the xa vector
% and store them in pz.

if nargin==5, % If no value of h was specified.
  if adapt==0, % Call either adaptive or non-adaptive routine.
    pz=kdest1d(xa,A,k);
    else,
      pz=akdest1d(xa,A,k);
  end;  
  else,
    if adapt==0, % If h was specified. 
      pz=kdest1d(xa,A,k,h);
      else,
        pz=akdest1d(xa,A,k,h);
    end; 
end;

figure('Name','Repeated sampling and K.d.e. plots');

sumI=0;

% Begin loop around nine different samples of size SS.

for i=1:9,
  
  s=sample(A,ss,0);

% Calculate kernel density estimate values corresponding to the xa vector
% and store them in sz.

  if nargin==5, % If no value of h was specified.
    if adapt==0, % Call either adaptive or non-adaptive routine.
      sz=kdest1d(xa,s,k);
      else,
        sz=akdest1d(xa,s,k);
    end;  
    else,
      if adapt==0, % If h was specified. 
        sz=kdest1d(xa,s,k,h);
        else,
          sz=akdest1d(xa,s,k,h);
      end; 
  end;

% Evaluate an estimate of I(f;g) using basic numerical integration.  

  v=sz.*log(sz./pz);
  I=xinc*sum(v);
  sumI=sumI+I;

  subplot(3,3,i);

  plot(xa,sz,'w'); % Plot a 2D k.e. function on the range defined by sa. 
  hold on;
  
  title(['I(f;g) = ',num2str(I)]);
 
  plot(xa,pz,'w--'); % Plot a 2D k.e. function on the range defined by xa. 

end;

fprintf('Average I value = %g\n',sumI/9)