ekozinec2.m

该工具箱是用于统计模式识别的

function [alpha,theta,solution,t,alpha1,alpha2]=ekozinec2(X,I,epsilon,tmax)
% EKOZINEC2 finds epsilon-optimal decision hyperplane, dichotomy.
%  [alpha,theta,solution,t,alpha1,alpha2]=ekozinec2(X,I,epsilon,tmax)
%
% EKOZINEC2 is faster version of EKOZINEC (see help ekozinec).
%
%  Input:
%   X [DxM] contains M training points in D-dimensional the feature 
%      space. X=[x1,x2,..XM] where xi is i-th column vectors.
%   I [1xM] contains class labels of the points in X. A class label 
%      has  to be 1 for the first class and 2 for the second class.
%   epsilon [1x1] desired quality of the solution.
%   tmax [1x1] maximal number of iterations.
%
%  Output:
%   alpha [Dx1] found normal vector of the hyperplane or, if tmax==-1,
%      badly classified point.
%   theta [1x1] found threshold of the hyperplane.
%   solution [1x1], 1 ... solution is found,
%                   0 ... solution is not found.
%   t [1x1] is number of steps the algorithm performed.
%
% See also EKOZINEC, KOZINEC, PERCEPTR, LINSVM.
%

% Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac
% (c) Czech Technical University Prague, http://cmp.felk.cvut.cz
% Written Vojtech Franc (diploma thesis) 02.11.1999, 13.4.2000
% Modifications
% 21-apr-2001, V.Franc, created

% Process input arguments
if nargin < 5,
   t=0;
   alpha1=0;
   alpha2=0;
end
if nargin < 4,
   tmax=inf;
end


X1=X(:,find(I == 1));
X2=X(:,find(I == 2));

n1 = size(X1,2);
n2 = size(X2,2);

alpha1=X1(:,1);
alpha2=X2(:,2);

solution=0;
t = 0;
while solution == 0 & tmax > t,

  t = t + 1;
  solution = 1;


  da1 = (alpha1 - alpha2);
  nda1 = norm( da1 );
  
  XA2 = X1 - repmat( alpha2, 1, n1 );
  
  [minXDA1, inx] = min( XA2' * da1 /nda1 );
  
  if minXDA1 <= (nda1 - epsilon/2 ),
    xt = X1(:, inx );

    k = min(1, (da1'*(alpha1-xt)) / ((alpha1-xt)'*(alpha1-xt)));
    
    alpha1 = alpha1 * (1 - k) + xt * k;
    
    
    solution = 0;
  else
    da2 = (alpha2 - alpha1 );
    nda2 = norm( da2 );
  
    XA1 = X2 - repmat( alpha1, 1, n2 );
  
    [minXDA2, inx ] = min( XA1' * da2  / nda2 );
  
    if minXDA2 <= (nda2 - epsilon/2 ),
      xt = X2(:, inx );

      k = min(1, (da2'*(alpha2-xt)) / ((alpha2-xt)'*(alpha2-xt)));

      alpha2 = alpha2 * (1 - k) + xt * k;
      
      solution = 0;
    end
  end
    
end

theta = 0.5 * ( alpha1'*alpha1 - alpha2'*alpha2 );
alpha = alpha1 - alpha2;

return;