svmnlusex01.m

数据挖掘的新方法－支持向量机书中算法例子

%SVMnLSPex01.m
%(Another expression of Object and constraint Function With Matrix Style)
%Two Dimension SVM Problem, Two Class and Un-separable Situation
%
%Method from Thorsten Joachims:
%"Making Large-Scale SVM Learning Practical", function (4)/(5)/(6)
%	
%  Objective:   min "f(A)=-sum(A)+(A'*Q*A)/2" ,							function (4);
%					 here	Q = diag(Y)*(X*X')*diag(Y) and in matrix Q any component Q(ij)=yi*yj*xi*xj';
%  Subject to:  sum{A'*Y}=0 ,													function (5);
%					 0 <= ai <= C ,												function (6);
%					 here yi={1,-1} indicates a sample xi belongs to either class of two ;
%The optimizing variables is "Lagrange Multipliers": A=[a1,a2,...,am],m is the number of total samples.
%

function mm=SVMnLSex01(mm) %mm -- select sample number

echo off;
close all;
fclose all;

%Generate Two Set of Data each set with m samples
	m1=mm;					% The number of class A which yi = 1 ;
   m2=mm;					% The number of class B which yi = -1;
   n = 2;					% The dimension of sample vector, that is xi = [xi(1), xi(2)];
   al=pi/8;					%
   
   x=rand(m1,2);
	x(1:m1,1)=2*x(1:m1,1);
	x(1:m1,2)=x(1:m1,2)+2.5+sin(x(1:m1,1)*pi);
	X1(1:m1,1)=x(1:m1,1)*cos(al)-x(1:m1,2)*sin(al);
	X1(1:m1,2)=x(1:m1,2)*cos(al)+x(1:m1,1)*sin(al);

	x=rand(m2,2);
	x(1:m2,1)=2*x(1:m2,1);
	x(1:m2,2)=x(1:m2,2)+1+sin(x(1:m2,1)*pi);
	X2(1:m2,1)=x(1:m2,1)*cos(al)-x(1:m2,2)*sin(al);
	X2(1:m2,2)=x(1:m2,2)*cos(al)+x(1:m2,1)*sin(al);
   
   clear x;
   
   X=[X1;X2];							%X1 is Positive Samples Set and X2 is Negative Samples Set;
   m=m1+m2;								%Total number of samples in matrix X(m*n);
     
   figure(1);
  	whitebg(1,[0,0.1,0]);
   plot(X(1:m1,1),X(1:m1,2),'m*');
   hold on
   plot(X(m1+1:m,1),X(m1+1:m,2),'c*');
   
   Xa = min(0,min(X(1:m,1)));		%Minimum of X axes of X
   Xb = max(X(1:m,1));				%Maximum of X axes of X
   Ya = min(0,min(X(1:m,2)));		%Minimum of Y axes of Y
   Yb = max(X(1:m,2));				%Maximum of Y axes of Y

	axis([Xa,Xb,min(X(1:m,2)),max(X(1:m,2))])
   
   asw = 0;
   while (asw~=1)&(asw~=2)
      asw=input('Continue or no? (1/2)');
   end
   if asw==1
   	%
		%
		%Transmiting Samples to Object and Constraint Function, in SVMnLSex01FUN.m,
		fid=fopen('a.dat','w');
      fwrite(fid,[m1,m2,n],'float');
      for k=1:n
         fwrite(fid,X(1:m,k),'float');
      end
      fclose(fid);
		%
		%Keep to MATLAB prescribe,Here
		%x0 is start point to optimizing;
      %Set options(13)=1 indicate the first constraint function contained in SVMnLSPex01FUN.m is a equality,
      %and others are "<=0",Reference to MATLAB function "CONSTR". Of cause, here we have only one constraint;
      %
   	x0 = 0.5*ones(m,1);	%Initial point of optimizing					
      options(13) = 1;
      
	   C = 500;					
      VLB = zeros(m,1);			%	The lower-bound of Lagrange Multiplier ai
      VUB = C*ones(m,1);		%	The up-bound of Lagrange Multiplier ai
%
% About The up-bound of Lagrange Multiplier:  
% A Multiplier a(i) isn't zero indicate a corresponding constranit in the original optimizing problem is active.
% That is the constraint yi*(W*xi +b )-1>=0 is active, on the other words, sample (xi) is a SVs.
% The other Multiplier aj = 0 corresponding to the samples who aren't SV.
%
% But, too small up-bound imposed to Multipliers will cause a Mutiplier who corrrespond un-active constraint 
% cann't descent to zeros in optimizing procedure. 
% However, we decide if a sample is SV just only dependent on the fact that its Multiplier is zero or no,so
% this will cause difficult when distinguish SVs from samples.
     
		[A,options,L] = constr('SVMnLSPex01FUN', x0, options, VLB, VUB);
      %
		%Return parameter: 
		%  A=[A(i)] is the optimizing variables, Lagrange Multipliyer;

	   Y=[ones(m1,1);-ones(m2,1)];	%The classfication value for both Positive and Negative Samples Set
		W = sum((diag(A)*diag(Y))*X);	%Normal vector of the superplane who divides two classes with largest margin.

		[z,I] = max(A(1:m));			%The Support Vectors are corresponding to a(i) <> 0;
      XI = X(I,1:n)';				%Pick out a arbitrary Support Vector to determine the "b" in classifying function
      YI = Y(I);
      
      %yi*(W*xi +b )-1>=0, when yi=1 then k = -W(1)/W(2) and b = yi/W(2)-W*X0/W(2);
      %The classifing function should be X*W + b = 0. It is y = k*x - b, here k and b are,
      b = YI - (diag(A)*Y)'*K(X,XI); 
      Yk = [Ya:0.025:Yb];
      Xd = [Xa:0.025:Xb+0.0125];
      for k = 1:length(Xd),
         Xy = [Xd(k)*ones(1,length(Yk));Yk];
         fx = (diag(A)*Y)'*K(X,Xy) + b;
         [l,p] = min(abs(fx));
         Yd(k) = Yk(p);
      end
      
      plot(Xd,Yd,'y')
         
   	for n=1:m
 	     if A(n) > 1e-6,
            plot(X(n,1),X(n,2),'wo');
 	     end
	   end
      
      delete('a.dat');

end

   function Kernel = K(U,E)									%Kernel function
      Kernel = (U*E+1).^5;