svm168.m

各种SVM分类算法

function [nsv,al3,bi,T]=svm168(x,y,ker,C,par,tol);
%
%--------------------S200502106  郑黎方------------------------------------
% 支持向量机在模式识别方面的应用.    
%
%     [nsv,al3,bi,T]=svm168(x,y,ker,C,par,tol);
%
% 
% 输入参数:
%
%              x: 输入矩阵，每行代表一个样本，列数表示样本特征数。矩阵行数代表样本数 。
%              y: 由 -1 或 +1组成的列向量. 每行表示相应样本的类别。只有两类。
%              ker: 字符串 'linear', 'poly_h', 'poly_i' and 'rbf'。
%                   用于传递内积函数的类型。
%              C: 惩罚因子。
%              par: 内积函数的参数。如多项式内积表示多项式的次数。查看kernel函数。
%              tol: 算法的迭代误差。默认值是 10^-5。 
%
% 输出:
%
%              nsv: 支持向量个数.
%              al3: alpha值.
%              bi: 偏差值.
%              T: 运行时间.
%  
%
% 例子:
%
%      线性内积函数:(见svmtest.m)
%               clear;clc;close all;
%               disp('输入样本矩阵(每行表示一个测试数据):');
%               x=[randn(1,20)-3 randn(1,20)+3;randn(1,40)]'
%               disp('样本所属类别:');
%               y=[-ones(1,20) ones(1,20)]'
%               [nsv,alpha,bias,T]=svm168(x,y,'linear',10);
%               fprintf('运行时间T=%g\n',T);
%               fprintf('支持向量个数：%g\n',nsv);
%
%               svcplot(x,y,'linear',alpha,bias);
%
%               W=x'*alpha;
%               disp('权系数向量:');
%               W
%               disp('测试数据:');
%               X_test=[randn(1,20)-3 randn(1,20)+3;randn(1,40)]'
%               disp('测试数据类别:');
%               Y_test=[-ones(1,20) ones(1,20)]'
%               disp('测试数据类别判决:');
%               Output=X_test*W+bias
%               disp('错误率:');
%               number_error=sum(abs(Y_test-sign(Output)))/2;
%               error=number_error/20
%
%               
%     径向基(或多项式)内积函数:(testrbf.m)
%               clear;clc;close all;
%               disp('输入样本矩阵(每行表示一个测试数据):');
%               x=[randn(1,20)-3 randn(1,20)+3;randn(1,40)]'
%               disp('样本所属类别:');
%               y=[-ones(1,20) ones(1,20)]'
%               [nsv,alpha,bias,T]=svm168(x,y,'rbf',10,4);
%               fprintf('运行时间T=%g\n',T);
%               fprintf('支持向量个数：%g\n',nsv);
%
%               svcplot(x,y,'linear',alpha,bias);
%
%               disp('测试数据:');
%               X_test=[randn(1,20)-3 randn(1,20)+3;randn(1,40)]'
%               disp('测试数据类别:');
%               Y_test=[-ones(1,20) ones(1,20)]'
%               disp('测试数据类别判决:');
%               Output=kernel('rbf',X_test,x(find(alpha),:),4)*alpha(find(alpha))+bias
%               disp('错误率:');
%               number_error=sum(abs(Y_test-sign(Output)))/2;
%               error=number_error/20
%
%
%--------------------------------------------------------------------------

T=clock;
if (nargin<4 | nargin>6)
   help svm168
else
   if(nargin==4 & ker(1)=='p')
       disp('将多项式内积的次数置为2。');
       par=2;
   elseif(nargin==4 & ker(1)=='r')
       fprintf(1,'径向基函数内积的参数（方差）置为 %1.3f\n',sqrt(size(x,2)));
       par=sqrt(size(x,2));
   elseif(nargin==4)
       par=0;
   end
   if(nargin<6)
      tol=10^-5;
   end         
   K=10^6;
   N=size(x,1);
   ns=-1;
   k=2;
   hacer=1;
   bi=0;
   
   H=kernel(ker,x,x,par);
   
   i1p=1:2:N;
   i1n=2:2:N;
   i1=[i1p i1n]';   
   
   a=zeros(N,1);
   a(i1)=C;
   
   al3=zeros(N,1);
   bi=0;
   i2p=[];
   i2n=[];
   
   Lp(1)=N*C;
   e=y;
   while(hacer)
      al3_a=al3;
      bi_a=bi;
      al3=zeros(N,1);
      al3(i2p)=C;
      al3(i2n)=-C;
      if(length(i1))
         Xi=inv([H(i1,i1)+diag(1./(a(i1))) ones(length(i1),1);ones(1,length(i1)) 0]);
         aux=Xi*[y(i1)+C*H(i1,i2n)*ones(length(i2n),1)-C*H(i1,i2p)*ones(length(i2p),1);C*length(i2n)-C*length(i2p)];
         al3(i1)=aux(1:length(i1));
         bi=aux(length(i1)+1);
      end
      e_a=e;
      
      I=find(al3-al3_a);      
      if(length(I))
         e=e_a-H(:,I)*(al3(I)-al3_a(I))-(bi-bi_a);
      end      
      
      Lp(k)=al3'*H*al3/2+C*sum(e(find(e.*y>0)).*y(find(e.*y>0)));
 
      if((abs(Lp(k)-Lp(k-1))/Lp(k-1))<tol & k>(ns+1))
         nsv=length(i1)+length(i2n)+length(i2p);
         hacer=0;                  
      else                     
          if((Lp(k)-Lp(k-1))/Lp(k-1)>tol)
              LL=Lp(k-1);
              HH=Lp(k);
              posL=0;
              posH=1;
              while(posL==0)
                  jj=(posH+posL)/2;
                  al3b=al3*jj+al3_a*(1-jj);          
                  bib=bi*jj+bi_a*(1-jj);
                  I=find(al3b-al3_a);
                  if(length(I))
                      eb=e_a-H(:,I)*(al3b(I)-al3_a(I))-(bib-bi_a);              
                  end
                  if(HH>LL)
                      HHa=HH;
                      posHa=posH;
                      HH=al3b'*H*al3b/2+C*sum(eb(find(eb.*y>0)).*y(find(eb.*y>0)));
                      posH=jj;
                  else
                      LL=al3b'*H*al3b/2+C*sum(eb(find(eb.*y>0)).*y(find(eb.*y>0)));
                      posL=jj;
                  end
              end
              jj=posH+posL;
              al3b=al3*jj+al3_a*(1-jj);          
              bib=bi*jj+bi_a*(1-jj);
              I=find(al3b-al3_a);
              if(length(I))
                  eb=e_a-H(:,I)*(al3b(I)-al3_a(I))-(bib-bi_a);              
              end
              valor=al3b'*H*al3b/2+C*sum(eb(find(eb.*y>0)).*y(find(eb.*y>0)));
              
              
              posL=posH;
              LL=HH;
              posH=jj;
              HH=valor;
              if(valor>Lp(k-1))
                  while(abs(valor-Lp(k-1))/Lp(k-1)>10^-3)                      
                      jj=(posH+posL)/2;
                      al3b=al3*jj+al3_a*(1-jj);          
                      bib=bi*jj+bi_a*(1-jj);
                      I=find(al3b-al3_a);
                      if(length(I))
                          eb=e_a-H(:,I)*(al3b(I)-al3_a(I))-(bib-bi_a); 
                      end
                      valor=al3b'*H*al3b/2+C*sum(eb(find(eb.*y>0)).*y(find(eb.*y>0)));                                            
                      if(valor>Lp(k-1))                    
                          HH=valor;
                          posH=jj;
                      else
                          LL=valor;
                          posL=jj;
                      end
                  end
              end                                   
              pos=jj;
              
              al3=al3*pos+al3_a*(1-pos);
              bi=bi*pos+bi_a*(1-pos);
              I=find(al3-al3_a);      
              if(length(I))
                  e=e_a-H(:,I)*(al3(I)-al3_a(I))-(bi-bi_a);
              end  
              Lp(k)=al3'*H*al3/2+C*sum(e(find(e.*y>0)).*y(find(e.*y>0)));
 
              if(pos<10^-10)
                  hacer=0;
                  nsv=length(i1)+length(i2n)+length(i2p);
              end
          end
                     
          i1p=find(e>=0 & y==1);
          i1n=find(e<=0 & y==-1);
          
         a=zeros(N,1);
         a(i1p)=C./e(i1p);
         a(i1n)=-C./e(i1n);
         a(find(a>K))=K;

         i2p=find(al3>=(C-C/100) & al3<=(C+C/100) & al3_a>=(C-C/100) & al3_a<=(C+C/100) & e>0 & y==1 & 15*e>14*e_a);
         if(length(setdiff(i1p,i2p)))
             i1p=setdiff(i1p',i2p)';
         else
             [val,pos]=min(e(i2p));
             i2p=setdiff(i2p,i2p(pos));
             i1p=setdiff(i1p',i2p)';
         end
         
         i2n=find(al3>=(-C-C/100) & al3<=(-C+C/100) & al3_a>=(-C-C/100) & al3_a<=(-C+C/100) & e<0 & y==-1 & 15*e<14*e_a);
         if(length(setdiff(i1n,i2n)))
             i1n=setdiff(i1n',i2n)';
         else
             [val,pos]=max(e(i2n));
             i2n=setdiff(i2n,i2n(pos));
             i1n=setdiff(i1n',i2n)';
         end
         
         i1=[i1p;i1n];
         
         k=k+1;
      end   
   end
end
T=etime(clock,T);