function [xa,ya,f,h,hl]=kdplot(A,k,inc,adapt,h,g);

% KDPLOT      Kernel Density Estimate plot for 1D or 2D data.
%
%             Example: load kddata
%                      kdplot(suicide,'knorm',100,0,0);
%                      kdplot(suicide,'kbiwe',200,0,40);
%                      [x,y,z,h]=kdplot(suicide,'knorm',100,1,0);
%
%             KDPLOT(A,K,INC,ADAPT,H,G) plots the Kernel Density 
%             Estimate with kernel given by the function in K.
%             'KNORM' is the multivariate normal distribution.
%             'KEPAN' is the Epanechnikov kernel.
%             'KEPV2' and 'KEPV3' are described on page 76 of Silverman (1986) 
%                  and are only valid for bi-variate data. 
%             'KTRIA' is the Triangular kernel.
%             'KRECT' is the rectangular kernel.
%             'KBIWE' is the Bi-weight kernel. 
%
%             The data is stored in the rows of A.  For 1D data A 
%             will have one column, for 2D data A will have 2 columns.
%
%             INC is the number of x and y 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 parameter H represents the window width:
%             (1) If a single non-zero value of H is given then smoothing is the 
%             same in directions.
%             (2) If H=[H1, H2] then smoothing is by H1 in the x direction and H2
%             in the y direction.  If H=0 is specified, automatic 'optimal' choices 
%             are made.
%             (3) If H=[H1, H2, H3] then the smoothing can be in orientations
%             different from those of the co-ordinate directions.  Note: the size
%             of H3 must be less than H1*H2.   
%
%             Christian C. Beardah 1994

if nargin==5,
  g=1;
end;

sym='g';

[n d]=size(A);  % Establish dimensions of A,
                % n=number of data points,
                % d=dimension of the data, either 1 or 2.

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.

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

if d==2,        % If data is bi-variate.
  y=A(:,2);     % Extract the second column of A into y.
  ymin=min(y);
  ymax=max(y);
  yrange=ymax-ymin;

% Establish an 'axis' vector for the y values.

  ya=linspace(ymin-yrange/4,ymax+yrange/4,inc);

% X and Y are matrices which define the 2D grid of points at which
% the k.d. function is evaluated.

  [X,Y]=meshgrid(xa,ya);

end;

if d==1,         % If data is 1D.
   ya=[];
   if adapt==0, % Call either adaptive or non-adaptive routine.
    [f,h]=kdest1d(xa,A,k,h);
    else,
      [f,h,hl]=akdest1d(xa,A,k,h);
  end;  
  else,        % If data is 2D. 
    if adapt==0,
     [f,h]=kdest2d(X,Y,A,k,h);
      else,
        [f,h]=akdest2d(X,Y,A,k,h);
    end;  
end;

if d==1 & g~=0,       % If data is 1D.
  plot(xa,f,sym); % Plot a 2D k.e. function on the range defined by xa. 
  ylabel('Rel. Freq.');
  title(['h = ',num2str(h)]);
  elseif d==2 & g~=0,        % If data is 2D.
    mesh(xa,ya,f);          % Output a 3D mesh plot.
%    surf(xa,ya,f);          % Output a 3D surface plot.
    xlabel('x');
    ylabel('y');
    zlabel('Rel. Freq.');
    if length(h)==3,
      title(['h = ',num2str(h(1)),' , ',num2str(h(2)),' , ',num2str(h(3))]);
      elseif length(h)==2,
        title(['h = ',num2str(h(1)),' , ',num2str(h(2))]);
        else,
          title(['h = ',num2str(h)]);
    end;
end;
