nvldtc.m

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

% A Novelty Detection Algorithm with linear kernel function, The samples aren't mapped and computation is simply preceeded in the input space.  
% Referenced from; David M.J. Tax and P.W. Duin "Data Domain Description using Support Vectors"
% Parameters:
% n: number of samples
% C: punish factor
% Use a kernel function as K(x,y)=x*y';


function a=NvlDtc(n,C);

close all
Arg=nargin;
if Arg==0,
    n=25;
    C=0.27;
end
if Arg==1,
    C=0.27;
end


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

fid=fopen('a.dat','w');
    fwrite(fid,n,'float');
    fwrite(fid,X(1:n,1),'float');
    fwrite(fid,X(1:n,2),'float');
fclose(fid)

a0=ones(n,1)*C*0.5;             % The start of a when it is optimized
A=[];                           % No ineqation constraint
B=[];
Aeq=ones(1,n);                  % Aeq*a=Beq, the constraint sum(a)=1;
Beq=1;
LB=zeros(n,1);                  % 0<=a(i)<=C
UB=ones(n,1)*C;


a=fmincon(@NvlDtcFUN,a0,A,B,Aeq,Beq,LB,UB);

delete a.dat;

hold on

X0=[0,0];
nirec=0;
ni=0;
for i=1:n,
    X0=X0+a(i)*X(i,1:2);                    % The center of decision area
    if (abs(a(i))>0.00000001) & (a(i)<C),   % 0<a(i)<C according to a X(i,1:2) is support vector
        plot(X(i,1),X(i,2),'g*');
        plot(X(i,1),X(i,2),'gs');
        if nirec==0;
            ni=i;                           % Record one of the support vectors
        end
    end
    if abs(a(i))==C,                        % a(i)==C means the X(i,1:2) is out of decision area
        plot(X(i,1),X(i,2),'g*');
        plot(X(i,1),X(i,2),'wv');
    end
end

plot(X0(1),X0(2),'ko');         % The centre of decision area 
plot(X0(1),X0(2),'k+');         %  

RV=X(ni,1:2)-X0;                % A Radius vector of the decision area

R=RV*RV';                       % Radius of the decision area

% Detecting and plotting the outline of decision area
x=-4.9:0.1:4.9;
y=-4.9:0.1:4.9;
    for k=1:length(x),
        for m=1:length(y),
            Z=[x(k),y(m)]-X0;
            RZ=Z*Z';
            if abs(RZ-R)<=0.1,
               plot(x(k),y(m),'m.');
            end
         end
    end

hold off