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

% KDSSMPL     Evaluates I(f;g) for various partitions.
%
%             KDSSMPL(A,K,ADAPT,SS,H) evaluates I(f;g) for different
%             values of the plotting increment. 
%
%             SS is the sample size.
% 
%             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.
%
%             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

figure;

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.

s=sample(A,ss,0);

for inc=1:9,

% 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.

  inc=2^(inc);

  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==4, % 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;

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

  if nargin==4, % 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);

  fprintf('I value = %g\n',I)

  plot(inc,I,'wo');
  if inc==2,
    hold on;
  end;

end;