function [h,hvec,score]=hbcv(A,k,steps)

% HBCV        Biased cross-validation estimate of h.
%
%             Examples: load kddata;
%
%             (1) hbcv(suicide,'knorm');
%
%             (2) [h,hvec,score]=hbcv(maddiam,'knorm'); 
%                 plot(hvec,score)     
%
%             HBCV(A,K,STEPS)
%
%             Christian C. Beardah 1995 

n=length(A);

if nargin==2,
  steps=100;
end;

hvec=linspace(0.01,hns(A,k),steps);

for i=1:steps,

  [x,y]=kdfft1(A,'knorm',128,hvec(i));  

  inc=x(2)-x(1);

  y=diff(diff(y)/inc)/inc;
  
  term=inc*sum(y.^2)-0.2116/(n*hvec(i)^5);

  score(i)=1/(n*hvec(i)*2*sqrt(pi))+hvec(i)^4*term/4;

end;

[L,I]=min(score);

h=hvec(I);