function h=hste(A,k)

% HSTE        2-stage Solve the Equation estimate of h.
%
%             Examples: load kddata;
%                       hste(suicide,'knorm');
%
%             HSTE(A,K)
%
%             Christian C. Beardah 1995 

n=length(A);

inc=128;

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

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 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;

psi6=-15/(16*sqrt(pi)*s^7);

psi8=105/(32*sqrt(pi)*s^9);

k40=deriv(0,k);

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

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

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

% Estimate psi6.

% Obtain the kernel weights.

[kw4,kw6]=deriv((bx-ax)*[0:inc-1]/((inc-1)*g2),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*(n-1)*g2^7);

% Now estimate psi4.

% Obtain the kernel weights.

kw4=deriv((bx-ax)*[0:inc-1]/((inc-1)*g1),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*(n-1)*g1^5);

ho=0;

h=hns(A,k);

k40=deriv(0,k);

while abs(ho-h)/abs(h)>0.01,

  temp=h;

  gamma=((2*k40*mu2*psi4*h^5)/(-psi6*R))^(1/7);

% Now estimate psi4.

% Obtain the kernel weights.

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

% Apply 'fftshift' to kw.

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

% Perform the convolution.

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

  p4=sum(c.*z(1:inc))/(n*(n-1)*gamma^5);

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

  ho=temp;

end;  

