f_svs.m

模式识别工具箱,希望对大家有用!

% F_SVS
%
%  [frac_SV2,alf,R2,Dx,I] = f_svs(sigma,x2,labx,frac_error);
%
% Compute the best fitting sphere around data x, using kerneltype
% gaussian. labx indicates if x is element of the target set (labx=1)
% or if it is an outlier (labx=-1). This improved version requires not
% the real data, but the squared distances between the objects (so, we
% need x2 = distm(x,x)).  The free parameter in the method
% is given in sigma. frac_SV2 is the fraction of the target set which
% becomes support vector (this is required for the function
% minimization routine 'fminbnd').
%
% The method returns in Dx the distance of the trainpoints to the
% centre of the sphere, R2 is the squared radius, alf are the optimized
% alfa's, I are the indices of the support vectors.
%
% frac_error  can also be a vector with length 2, such that different
% errors are allowed for the different classes. The first frac_error
% is used for class y=+1, the second frac_error for class y=-1.
%
% see also svdd_fmin, svdd.
% 19-4-2000  Dxd

% Copyright: D.M.J. Tax, R.P.W. Duin, [e-mail]
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands


function [frac_SV2,alf,R2,Dx,I] = f_svs(sigma,x2,labx,frac_error,fracsv);

if (nargin<5)
  fracsv = 0;
end
if (nargin<4)
  frac_error = [0, 0];
end

num_points = size(x2,1);

if (length(frac_error)<2)
  frac_error(2) = frac_error(1);
end

if ((frac_error(1)>=1.0)|(frac_error(2)>1.0))
  warning('gsvs: frac_error>=1.0: should I reject all data?');
end
if (frac_error(1)<=0)
  C1 = 1.0;
else
  C1 = 1/(num_points*frac_error(1));
end
if (frac_error(2)<=0)
  C2 = 1.0;
else
  C2 = 1/(num_points*frac_error(2));
end

if (size(labx,2)>1)
  error('svs: Please make labx a column vector! (containing 1/-1)');
end

K = exp(-x2/(sigma*sigma));
D = (labx*labx').*K;
f = labx.*diag(D);
% make sure D is positive definite:
i = -30;
while (pd_check(D + (10.0^i)*eye(num_points)) == 0)
  i = i+1;
%  disp('*')
end
i = i+5;
D = D + (10.0^i)*eye(num_points);

% equality constraints:
A = labx';
b = 1.0;

% lower and upper bounds:
lb = zeros(num_points,1);
ub = lb;
ub(find(labx==1)) = C1;
ub(find(labx==-1)) = C2;

% initialization (not sure if this is really necessary):
rand('seed', sum(100*clock));
p = [ 0.5*rand(num_points,1) ];

% other:
% these procedures *maximize* the functional L
if (exist('qld') == 3)
  alf = qld (2.0*D, -f, -A, b, lb, ub, p, 1);
else
  alf = quadprog(2.0*D,-f,[],[],A,b,lb,ub,p);
end

if (isempty(alf))
  warning('No solution for the SVDD could be found!');
  alf = ones(num_points,1)/num_points;
end

%important: change sign for negative examples:
alf = labx.*alf;

% the support vectors and errors:
I = find(abs(alf)>1e-8);

% distance to center of the sphere (ignoring the offset):
Dx = - sum( (ones(num_points,1)*alf').*K, 2);

% threshold squared radius:
borderx = I(find((alf(I) < ub(I))&(alf(I) > 0)));
if (size(borderx,1)<1)  % hark hark
  borderx = I;
end
R2 = mean(Dx(borderx,:));
% although all support vectors should give the same results, sometimes
% they do not. So I take the positive approach and make the rest of
% the SV inside the boundary:
%R2 = max(Dx(borderx,:));

% set all nonl-support-vector alpha to 0
J = find(abs(alf)<1e-8);
alf(J) = 0.0;

% what is the difference between the wanted and the current fraction
% of erroneously classified data?
diff = fracsv + frac_error(1) - length(I)/num_points;
frac_SV2 = diff.*diff;

%fprintf(' s=%f,  #SV=%d\n',sigma,length(I));
%fprintf(' goal=%f,  now=%f\n',fracsv+frac_error(1),length(I)/num_points);

return