【原創】支持向量機非線性回歸通用matlab程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    SVMNR.m
%    Support Vector Machine for Nonlinear Regression
%    支持向量機非線性回歸通用matlab程序
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Alpha1,Alpha2,Alpha,Flag,B]=SVMNR(X,Y,Epsilon,C,D)
%  ChengAihua,PLA Information Engineering University,ZhengZhou,China
%  Email:aihuacheng@gmail.com
%  All rights reserved
%本程序使用支持向量機法，實現對數據的非線性回歸
%輸入參數列表
%  X        輸入樣本,n×l的矩陣，n為變量個數，l為樣本個數
%  Y        輸出樣本,1×l的矩陣，l為樣本個數
%  Epsilon  ε不敏感損失函數的參數，Epsilon越大，支持向量越少
%  C        懲罰系數，C過大或過小，泛化能力變差
%  注意：核函數的設定和修改在函數內部進行，數據預處理在函數外部進行
%輸出參數列表
%  Alpha1   α系數
%  Alpha2   α*系數
%  Alpha    支持向量的加權系數（α－α*）向量
%  Flag     1×l標記，0對應非支持向量，1對應邊界支持向量，2對應標準支持向量
%  B        回歸方程中的常數項

%第一步：構造K矩陣
[n,l]=size(X);
K=zeros(l,l);
for i=1:l
    for j=1:l
        xi=X(:,i);
        xj=X(:,j);
        %K(i,j)=sum(xi.*xj);
        %K(i,j)=(sum(xi.*xj)+1)^20;%注意：核函數在此定義！
        K(i,j)=exp(-(sum((xi-xj).^2)/D));
        %K(i,j)=exp(-(sum((xi-xj))/100));
    end
end
%第二步：構造二次規劃模型參數H,Ft,Aeq,Beq,lb,ub
H=[K,-K;-K,K];
H=(H+H')/2;
Ft=[Epsilon*ones(1,l)-Y,Epsilon*ones(1,l)+Y];
Aeq=[ones(1,l),-ones(1,l)];
Beq=0;
lb=eps.*ones(2*l,1);
ub=C*ones(2*l,1);
%第三步：調用優化工具箱quadprog函數求解二次規劃
OPT=optimset;
OPT.LargeScale='off';
OPT.Display='off';
[Gamma,Obj]=quadprog(H,Ft,[],[],Aeq,Beq,lb,ub,[],OPT);
%[Gamma,Obj]=fmincon('myfun',5.*ones(18,1),[],[],Aeq,Beq,lb,ub);
%第四步：整理輸出參數
Alpha1=(Gamma(1:l,1))';
Alpha2=(Gamma((l+1):end,1))';
Alpha=Alpha1-Alpha2;
Flag=2*ones(1,l);
%第五步：支持向量的分類
Err=0.0000000000001;
for i=1:l
    AA=Alpha1(i);
    BB=Alpha2(i);
    if (abs(AA-0)<=Err)&&(abs(BB-0)<=Err)
        Flag(i)=0;
    end
    if (AA>Err)&&(AA<C-Err)&&(abs(BB-0)<=Err)
        Flag(i)=2;
    end
    if (abs(AA-0)<=Err)&&(BB>Err)&&(BB<C-Err)
        Flag(i)=2;
    end
    if (abs(AA-C)<=Err)&&(abs(BB-0)<=Err)
        Flag(i)=1;
    end
    if (abs(AA-0)<=Err)&&(abs(BB-C)<=Err)
        Flag(i)=1;
    end
end
%第六步：計算B
B=0;
counter=0;
for i=1:l
    AA=Alpha1(i);
    BB=Alpha2(i);
    if (AA>Err)&&(AA<C-Err)&&(abs(BB-0)<=Err)
        %計算支持向量加權值
        SUM=0;
        for j=1:l
            if Flag(j)>0
                %SUM=SUM+Alpha(j)*sum(X(:,j).*X(:,i));
                %SUM=SUM+Alpha(j)*(sum(X(:,j).*X(:,i))+1)^20;
                SUM=SUM+Alpha(j)*exp(-(sum((X(:,j)-X(:,i)).^2)/D));
                %SUM=SUM+Alpha(j)*exp(-(sum((X(:,j)-X(:,i)))/100));
            end
        end
        b=Y(i)-SUM-Epsilon;
        B=B+b;
        counter=counter+1;
    end
    if (abs(AA-0)<=Err)&&(BB>Err)&&(BB<C-Err)
        SUM=0;
        for j=1:l
            if Flag(j)>0
                %SUM=SUM+Alpha(j)*sum(X(:,j).*X(:,i));
                %SUM=SUM+Alpha(j)*(sum(X(:,j).*X(:,i))+1)^20;
                SUM=SUM+Alpha(j)*exp(-(sum((X(:,j)-X(:,i)).^2)/D));
                %SUM=SUM+Alpha(j)*exp(-(sum((X(:,j)-X(:,i)))/100));
            end
        end
        b=Y(i)-SUM+Epsilon;
        B=B+b;
        counter=counter+1;
    end
end
B=B/counter;