function h=hdpi(A,k,l);

% HDPI        L-stage Direct plug-in estimate of smoothing parameter.
%
%             Example: load kddata
%                      hdpi(suicide,'knorm',1);
%
%             HDPI(A,K,L)  
%
%             L can be 0, 1, 2 or 3.
%
%             Christian C. Beardah 1994

if nargin==2,
  l=0;
end;

inc=128;

P=inc*2;

n=length(A);

xmin=min(A);    % Find the minimum value of A.
xmax=max(A);    % Find the maximum value of A.
xrange=xmax-xmin; % Find the range of A.

if k=='knorm',
  mu2=1;
  R=1/(2*sqrt(pi));
  r=4;
end;

if siqr(A)>0,
  s=min([std(A) siqr(A)/1.34]);
  else,
    s=std(A);
end;

% xa holds the x 'axis' vector, defining a grid of x values where 
% the k.d. function will be evaluated and plotted.

ax=xmin-xrange/8;
bx=xmax+xrange/8;

xa=linspace(ax,bx,inc); 

c=zeros(inc,1);
deltax=(bx-ax)/(inc-1);
binx=floor((A-ax)/deltax)+1;

% Obtain the grid counts.

for i=1:n, % Loop over data points in x direction.
  c(binx(i))=c(binx(i))+(xa(binx(i)+1)-A(i))/deltax;
  c(binx(i)+1)=c(binx(i)+1)+(A(i)-xa(binx(i)))/deltax;  
end;

if l==3,

  psi10=-945/(64*sqrt(pi)*s^11);

  [k40,k60,k80]=deriv(0,k);

  g0=(-2*k80/(mu2*psi10*n))^(1/11);

% Obtain the kernel weights.

  [kw4,kw6,kw8]=deriv((bx-ax)*[0:inc-1]/((inc-1)*g0),k);

% Apply 'fftshift' to kw.

  kw=[kw8,0,kw8([inc:-1:2])]';

% Perform the convolution.

  z=real(ifft(fft(c,P).*fft(kw)));
 
  psi8=sum(c.*z(1:inc))/(n^2*g0^9);

end;

if l==2 | l==3,

  if l==2,
    psi8=105/(32*sqrt(pi)*s^9);
  end;

  [k40,k60]=deriv(0,k);

  g1=(-2*k60/(mu2*psi8*n))^(1/9);

% Obtain the kernel weights.

  [kw4,kw6]=deriv((bx-ax)*[0:inc-1]/((inc-1)*g1),k);

% Apply 'fftshift' to kw.

  kw=[kw6,0,kw6([inc:-1:2])]';

% Perform the convolution.

  z=real(ifft(fft(c,P).*fft(kw)));

  psi6=sum(c.*z(1:inc))/(n^2*g1^7);

end;

if l==1 | l==2 | l==3,

  if l==1,
    psi6=-15/(16*sqrt(pi)*s^7);
  end;

  k40=deriv(0,k);

  g2=(-2*k40/(mu2*psi6*n))^(1/7);

% Obtain the kernel weights.

  kw4=deriv((bx-ax)*[0:inc-1]/((inc-1)*g2),k);

% Apply 'fftshift' to kw.

  kw=[kw4,0,kw4([inc:-1:2])]';

% Perform the convolution.

  z=real(ifft(fft(c,P).*fft(kw)));

  psi4=sum(c.*z(1:inc))/(n^2*g2^5);

end;

if l==0,
  psi4=3/(8*sqrt(pi)*s^5);
end;

h=(R/(mu2^2*psi4*n))^(1/5);