nvldtckn.m

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

% A Novelty Detection Algorithm with radius based kernel function, every samples are mapped in a super-sphere;
% Referenced from; David M.J. Tax and P.W. Duin "Data Domain Description using Support Vectors"
% Parameters:
% n: number of samples
% C: punish factor
% s:
% Use a kernel function as K(x,y)=exp{-[(x-y)/s]^2};


function Nothing=NvlDtc(n,s);

Nothing=[];

close all
Arg=nargin;
if Arg==0,
    n=20;
    s=0.5;
end
if Arg==1,
    s=0.5;
end
C=1;

% Produce a sample set centered (0,0) 
X=randn(n,2);                   % Sample Set
figure(1);
whitebg(1,'k');
plot(X(1:n,1),X(1:n,2),'r*');
axis([-3,3,-3,3])

 asw=input('Is the sample distribution OK ? (1/0)');

if asw==1,
% Learning
%--------------------------------------------------------------------
    again=1;
    while again==1,

        % Sending sample data to subroutine NvlDtcKnFUN
        fid=fopen('a.dat','w');
            fwrite(fid,s,'float');
            fwrite(fid,n,'float');
            fwrite(fid,X(1:n,1),'float');
            fwrite(fid,X(1:n,2),'float');
        fclose(fid)

        % Set parameters for the optimizing M-function
        a0=ones(n,1)*C*0.5;             % The initial value of a, Lagrange factor, optimizing variable in Eq.(7) 
        A=[];                           % No ineqation constraint
        B=[];
        Aeq=ones(1,n);                  % Aeq*a=Beq, according to the constraint sum(a)=1;
        Beq=1;
        LB=zeros(n,1);                  % 0<=a(i)<=C, according to the constraint 0<a(1)<C, 0<a(2)<C, ...
        UB=ones(n,1)*C;

        a=fmincon(@NvlDtcKnFUN,a0,A,B,Aeq,Beq,LB,UB); % a,the Lagrange factor, optimizing variables here

        delete a.dat;

        hold on

        % Detecting and plotting the outline of decision area
        I=ones(1,n);
        x1=X(1:n,1);        % Split learning samples into their coordinate value to take the computation facility; 
        x2=X(1:n,2);
        T=1+sum(sum( diag(a)*(  exp(-((x1*I-(x1*I)').^2 + (x2*I-(x2*I)').^2)/s^2)  )*diag(a) )); % First and last items in Eq.(8),note that K(z,z)=1

        % Drawing the decision area
        x=-2.9:0.1:2.9;
        y=-2.9:0.1:2.9;
        for k=1:length(x),
            for m=1:length(y),
                Z=sum(([x(k)*I;y(m)*I]-X').^2)/(s^2);   % A row vector: ([X(i)-X(j)]/s)^2, i,j = 1,2,...
                D=T - 2*(  exp( -Z)*a  );               % Distance of [x,y] from the center of decision area,which is unable presented clearly
                if abs(abs(D))<= 1,
                    plot(x(k),y(m),'b.');               
                end
            end
        end

        % Plot sample set
        plot(X(1:n,1),X(1:n,2),'r*');
        axis([-3,3,-3,3])
  
        hold off
    
     again=input('Do again to the sample set with another s ? (1/0)');
     if again==1,
        s=input('input new s :');
        figure(2);
     end
    
    end % while

end % if