% 支援向量機SVM分類演算法

clear all
% ------------------------------------------------------------%
% 構造兩類訓練資料集     x2=aa*x1+bb+(-)b1

aa=3;
bb=6;
b1=0.2;       

x1(:,1) = -1:0.1:1;

n = length(x1(:,1));

x1(:,2) = aa.*x1(:,1) + bb + b1 + abs(randn(n,1));
y1 = ones(n,1);

x2(:,1) =  -1:0.1:1;
x2(:,2) = aa.*x2(:,1) + bb - b1 - abs(randn(n,1));

y2 = -ones(n,1);


figure;
plot(x1(:,1),x1(:,2),'bx',x2(:,1),x2(:,2),'k.');
hold on;

X = [x1;x2]; % 訓練樣本
Y = [y1;y2]; % 訓練目標,n×1的矩陣,n為樣本個數,值為+1或-1
% ------------------------------------------------------------%

tic
% 解二次優化方城
n = length(Y);
H = (Y*Y').*(X*X'); % liner kernel 
f = -ones(n,1);
A = [];
b = [];
Aeq = Y';
beq = 0;
lb = zeros(n,1);
ub = 100*ones(n,1);
a0 = zeros(n,1);
options = optimset;
options.LargeScale = 'off';
options.Display = 'off';
[a,fval,eXitflag,output,lambda] = quadprog(H,f,A,b,Aeq,beq,lb,ub,a0,options);
eXitflag

time=toc


   
%先複製上面的程式碼 執行後  再複製下面的程式碼 執行 呵呵 （對Matlab不熟，不知道為什麼一次執行和分開執行看到的結果不同）

% 以下是分類平面：
Y2=a.*Y;
 W(1)=sum(Y2.*(X(:,1)));
 W(2)=sum(Y2.*(X(:,2)));

aLarge=find(a>0.1);
j=aLarge(1);

S(:,1)=Y.*a.*X(:,1);
S(:,2)=Y.*a.*X(:,2);

S2=S*(X(j,:)');

b=Y(j)-sum(S2);


 xx1=x1(:,1);
 xx2=-(W(1)*xx1+b)/W(2);
 plot(xx1,xx2);