linreg.m

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

% A linear SVM regression example;
% Reference from "A Tutorial on Support Vector Regression",pp4,Eq.9;
% Parameters:
%           n: number of samples (random produced);
%           e: tube width in the algorithm (Fig.1);
%           x1,x2,y1,y2: expectative line of the linear distributed samples;
%
% Another default parameter:
%           variance of the random data has been set as s=(y2-y1)/8;
%           the punish parameter C=5
%
% Note: Sometime it is some wrong;

function [w,b]=LinReg(n,e,C,x1,x2,y1,y2,s)

Close all
%
Arg=nargin;
if Arg==0,
    n=20;
    e=2;
    C=5;
    x1=1;
    x2=5;
    y1=3;
    y2=50;
end
if Arg==1;
    e=2;
    C=5;
    x1=1;
    x2=5;
    y1=3;
    y2=50;
end
if Arg==2;
    C=5;
    x1=1;
    x2=5;
    y1=3;
    y2=50;
end
if Arg==3;
    x1=1;
    x2=5;
    y1=3;
    y2=50;
end
if Arg<8,
    s=(y2-y1)/5;
end


% Produce linear distributed samples:Variance=(y2-y1)/3
if x1>x2,
    a=x2;
    x2=x1;
    x1=a;
end
X=(x1+rand(1,n)*(x2-x1));
Y=randn(1,n)*s;
k=(y2-y1)/(x2-x1);
dy=y1+k*(X-min(X));
Y=Y+dy;

figure(1);
whitebg(1,'k')
plot(X,dy,'w--')
hold on
plot(X,Y,'g*')

% Sorting X and Y
[X,I]=sort(X);
Y1=Y;
for i=1:n,
    Y(i)=Y1(I(i));
end
X=X';
Y=Y';
clear Y1 I x1 x2 y1 y2 Arg dy;

% Setting optimal parameter
a=zeros(2*n,1);             % Original of optimal variables, all optimal variables are a(1:L) and a(L+1:2*L)
A=[];                       % No inequation constraints;
B=[];
Aeq=[ones(1,n),-ones(1,n)]; % Summing all optimal variables is equal to 0; 
Beq=0;
LB=zeros(2*n,1);            % Bound of variables 0<ai<C; 
UB=ones(2*n,1)*C;
a0=a;

%Optimization--------1;
%--------------------------------------
% Transfer necessary parameters to FUN
fid=fopen('a.dat','w');
    fwrite(fid,n,'float');
    fwrite(fid,e,'float');
    fwrite(fid,X,'float');
    fwrite(fid,Y,'float');
fclose(fid)

% Optimizing a (Eq.9)
a=fmincon(@LinRegFUN,a0,A,B,Aeq,Beq,LB,UB);

% Calculatting w (Eq.10) 
w=(a(1:n)-a((n+1):2*n))'*X;

    % Calculating b, (Eq.13)
    na=0;
    Sar=0;
    MoP=-1;
    a=round(a*100000)/100000;
    for i=1:2*n,
        if (a(i)>0.0000) & (a(i)<C), 
            if i>n,
                i=i-n;
                MoP=1;      % Indicate a* if MoP=-1 Indicate a;
            end
            na=na+1;
            plot(X(i),Y(i),'c^');
            if Sar==0;      % Indicate intercept b has not been calculated
                b=Y(i)-w*X(i)+MoP*e;
                Sar=1;  
            end
        end
    end

% The optimal regeression function
y=w*X+b;
yu=w*X+b+e;
yl=w*X+b-e;

plot(X,y,'r')
plot(X,yu,':r')
plot(X,yl,':r')
text(max(X),max(y),['e=',num2str(e)]);

delete a.dat
clear na a i bo Sar MoP yu yl k y yu yl s fid b w;

% Optimization--------2
%--------------------------------------
% Re-calculating with e=5*e;
e=3*e;
fid=fopen('a.dat','w');
    fwrite(fid,n,'float');
    fwrite(fid,e,'float');
    fwrite(fid,X,'float');
    fwrite(fid,Y,'float');
fclose(fid)

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

w=(a(1:n)-a((n+1):2*n))'*X;

    % Calculating b, (Eq.13)
    na=0;
    Sar=0;
    MoP=-1;
    a=round(a*100000)/100000;
    for i=1:2*n,
        if (a(i)>0.0000) & (a(i)<C), 
            if i>n,
                i=i-n;
                MoP=1;      % Indicate a* if MoP=-1 Indicate a;
            end
            na=na+1;
            plot(X(i),Y(i),'ws');
            if Sar==0;      % Indicate intercept b has not been calculated
                b=Y(i)-w*X(i)+MoP*e;
                Sar=1;  
            end
        end
    end
clear na a i bo Sar MoP;

hold on
y=w*X+b;
yu=w*X+b+e;
yl=w*X+b-e;

plot(X,y,'y')
plot(X,yu,':y')
plot(X,yl,':y')

delete a.dat;

text(min(X),min(y),['e=',num2str(e)]);