function h=hbcv2(A,k)

% HBCV2       Biased cross-validation for bi-variate data.
%          
%             Example: load kddata
%                      hbcv2(mask4,'knorm');
% 
%             HBCV2(A,K)
% 
%             Christian C. Beardah 1993-94 

n=length(A);

% Transform A so rows have st. dev.=1:

oldA=A;
A=(oldA-ones(n,1)*mean(oldA))*diag(1./std(oldA));
% A=(oldA-ones(n,1)*mean(oldA))*inv(sqrtm(cov(oldA)));

maxi=20;

tol=1e-2;

gradf=zeros(2,1);

delta=1e-6;

i=0;

vx(1)=hste(A(:,1),k);
vy(1)=hste(A(:,2),k);

x0=vx(1);
y0=vy(1);

fn(1)=bcv2(A,vx(1),vy(1));

i=1; % i is the iteration counter.

while i <= maxi,

  gradf(1)=(bcv2(A,x0+delta,y0)-bcv2(A,x0-delta,y0))/(2*delta);
  gradf(2)=(bcv2(A,x0,y0+delta)-bcv2(A,x0,y0-delta))/(2*delta);

  gradf=gradf/(i*norm(gradf));

  r0=[x0 y0]';

  for j=1:31,
    M(j,:)=(r0+(j-16)*gradf/15)';
  end;

  L=(M(:,1)>0) & (M(:,2)>0);

  L=remzero([M,L]);

  M=L(:,1:2);
  
  [rM,cM]=size(M);
  
  for j=1:rM,
    t(j,1)=bcv2(A,M(j,1),M(j,2));
  end;

  [N,I]=min(t);

  r=M(I,:)';

  vx(i+1)=r(1); % Store iterate.
  vy(i+1)=r(2); 

  fn(i+1)=N;

  if (abs(norm(r-r0)) < tol),

    h=r';

% Transform values of h back:

    h=diag(std(oldA))*h';
%    h=sqrtm(cov(oldA))*h';
    
    return;
  end;

  i=i+1; % Update the iteration counter.

  x0=r(1);
  y0=r(2);

  clear t
 
end; % of while statement.

error('Maximum number of iterations reached - terminating execution.');