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

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

if nargin==2,
  steps=100;
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.

s=min([std(A) siqr(A)/1.34]);

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

g1=sqrt(2)*s*(2/(7*n))^(1/9);

g2=sqrt(2)*s*(2/(11*n))^(1/13);

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

% Obtain the kernel weights.

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

% Apply 'fftshift' to kw.

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

% Perform the convolution.

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

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

g3=(-6/(sqrt(2*pi)*psi6*n))^(1/7);

g4=(-210/(sqrt(2*pi)*psi10*n))^(1/11);

% Obtain the kernel weights.

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

% Obtain the kernel weights.

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

H=hns(A,k);

hvec=linspace(0.01,1.5*H,steps);

C=(441/(64*pi))^(1/18)*(4*pi)^(-1/5)*psi4^(-2/5)*psi8^(-1/9);

M=A*ones(size(A'));

Y=(M-M');
 
for i=1:steps,

  g=C*n^(-23/45)*hvec(i)^(-2);

  sig1=sqrt(2*hvec(i)^2+2*g^2);

  sig2=sqrt(hvec(i)^2+2*g^2);

  sig3=sqrt(2*g^2);

  term2=sum(sum(feval(k,Y/sig1)/sig1-2*feval(k,Y/sig2)/sig2+feval(k,Y/sig3)/sig3));
  
  score(i)=1/(n*hvec(i)*2*sqrt(pi))+term2/n^2;

end;

[L,I]=min(score);

h=hvec(I);
