function h=hstt(A,k)

% HSTT        Scott-Tapia-Thompson estimate of h.
%
%             Examples: load kddata;
%                       hstt(suicide,'knorm');
%
%             HSTT(A,K)
%
%             Christian C. Beardah 1995 

n=length(A);

inc=512;

if k=='knorm',
  R=0.28209;
  mu2=1;
  elseif k=='kepan',
    R=0.6;
    mu2=0.2;
    elseif k=='kbiwe',
      R=0.71429;
      mu2=0.14286;
      elseif k=='ktria',
        R=2/3;
        mu2=1/6;
        elseif k=='krect',
          R=0.5;
          mu2=1/3;
          elseif k=='klapl',
            R=0.25;
            mu2=2;
            elseif k=='ktriw',
              R=350/429;
              mu2=1/9;
              elseif k=='klogi',
                R=1/6;
                mu2=pi^2/3;
end;

count=1;
maxit=20;

ho=0;

h=hns(A,k);

% This iteration can cycle (eg. Suicide data, Epanechnikov kernel).
% Don't allow more than maxit iterations.

while abs(ho-h)>0.01 & count<maxit,

  temp=h;

  [xa,y]=kdfft1(A,k,inc,h,0);

  delta=xa(2)-xa(1);

% Estimate psi4=R(f'') using simple finite differences and quadrature.

  for i=2:(inc-1),
    z(i)=((y(i+1)-2*y(i)+y(i-1))/delta^2)^2;
  end;

  psi4=delta*sum(z(2:(inc-1)));

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

  ho=temp;

  count=count+1;

end;