支援向量機（SVM）的matlab的實現

支援向量機是一種分類演算法之一，matlab中也有相應的函式來對其進行求解；下面貼一個小例子，這個例子來源於我們實際的專案。

clc;
clear;
N=10;
%下面的資料是我們實際專案中的訓練樣例（樣例中有8個屬性）
correctData=[0,0.2,0.8,0,0,0,2,2];
errorData_ReversePharse=[1,0.8,0.2,1,0,0,2,2];
errorData_CountLoss=[0.2,0.4,0.6,0.2,0,0,1,1];
errorData_X=[0.5,0.5,0.5,1,1,0,0,0];
errorData_Lower=[0.2,0,1,0.2
 
,0,0,0,0];
errorData_Local_X=[0.2,0.2,0.8,0.4,0.4,0,0,0];
errorData_Z=[0.53,0.55,0.45,1,0,1,0,0];
errorData_High=[0.8,1,0,0.8,0,0,0,0];
errorData_CountBefore=[0.4,0.2,0.8,0.4,0,0,2,2];
errorData_Local_X1=[0.3,0.3,0.7,0.4,0.2,0,1,0];
 sampleData=[correctData;errorData_ReversePharse;errorData_CountLoss;errorData_X;errorData_Lower;errorData_Local_X;errorData_Z;errorData_High;errorData_CountBefore;errorData_Local_X1]
 
;%訓練樣例

type1=1;%正確的波形的類別，即我們的第一組波形是正確的波形，類別號用 1 表示
type2=-ones(1,N-2);%不正確的波形的類別，即第2~10組波形都是有故障的波形，類別號用-1表示
groups=[type1 ,type2]';%訓練所需的類別號
j=1;
%由於沒有測試資料，因此我將錯誤的波形資料輪流從訓練樣例中取出作為測試樣例
for i=2:10
   tempData=sampleData;
   tempData(i,:)=[];
    svmStruct = svmtrain(tempData,groups);
    species(j) = svmclassify(svmStruct,sampleData(i
 
,:));
    j=j+1;
end
species

輸出結果如下：

 -1    -1    -1    -1    -1    -1    -1     1    -1

從結果可以看出：只有第九個被誤判，其它的都是正確的。

上面只是用於說明matlab中支援向量機中函式的用法，因為在訓練集中我只用了一個正確的波形和九組有故障的波形作為訓練集，因此這種超平面的選取可能不好。但是，在我們的實際的專案中，我們是用到了許多的訓練集的。

上面是呼叫matlab中的函式實現多維屬性中的支援向量機對其進行分類。

下面是自己實現此功能，不呼叫matlab中的支援向量機的函式。程式碼如下：

%主函式
clear all;
clc;
C = 10;
kertype = 'linear';
%訓練樣本
n = 50;
randn('state',6);%可以保證每次每次產生的隨機數一樣
x1 = randn(2,n);    %2行N列矩陣
y1 = ones(1,n);       %1*N個1
x2 = 5+randn(2,n);   %2*N矩陣
y2 = -ones(1,n);      %1*N個-1

figure;
plot(x1(1,:),x1(2,:),'bx',x2(1,:),x2(2,:),'k.'); 
axis([-3 8 -3 8]);
xlabel('x軸');
ylabel('y軸');
hold on;

X = [x1,x2];        %訓練樣本d*n矩陣，n為樣本個數，d為特徵向量個數，在這裡，X為一個2*100的陣列
Y = [y1,y2];        %訓練目標1*n矩陣，n為樣本個數，值為+1或-1，在這裡，Y為一個1*100的陣列
svm = svmTrain(X,Y,kertype,C);
plot(svm.Xsv(1,:),svm.Xsv(2,:),'ro');

%測試
[x1,x2] = meshgrid(-2:0.05:7,-2:0.05:7);  %x1和x2都是181*181的矩陣
[rows,cols] = size(x1);  
nt = rows*cols;                  
Xt = [reshape(x1,1,nt);reshape(x2,1,nt)];
Yt = ones(1,nt);
result = svmTest(svm, Xt, Yt, kertype);

Yd = reshape(result.Y,rows,cols);
contour(x1,x2,Yd,'m');

訓練集函式如下：

function svm = svmTrain(X,Y,kertype,C)
options = optimset;    % Options是用來控制演算法的選項引數的向量
options.LargeScale = 'off';%LargeScale指大規模搜尋，off表示在規模搜尋模式關閉
options.Display = 'off';%這樣設定意味著沒有輸出

n = length(Y);%陣列Y的長度
H = (Y'*Y).*kernel(X,X,kertype);%呼叫kernel函式，

f = -ones(n,1); %f為1*n個-1,f相當於Quadprog函式中的c
A = [];
b = [];
Aeq = Y; %相當於Quadprog函式中的A1,b1
beq = 0;
lb = zeros(n,1); %相當於Quadprog函式中的LB，UB
ub = C*ones(n,1);
a0 = zeros(n,1);  % a0是解的初始近似值
[a,fval,eXitflag,output,lambda]  = quadprog(H,f,A,b,Aeq,beq,lb,ub,a0,options);

epsilon = 1e-8;                     
sv_label = find(abs(a)>epsilon);  %0<a<a(max)則認為x為支援向量     
svm.a = a(sv_label);
svm.Xsv = X(:,sv_label);
svm.Ysv = Y(sv_label);
svm.svnum = length(sv_label);
%svm.label = sv_label;

核函式如下：

function K = kernel(X,Y,type)
%X 維數*個數
switch type
case 'linear'
    K = X'*Y;
case 'rbf'
    delta = 5;
    delta = delta*delta;
    XX = sum(X'.*X',2);%sum(a,2)程式碼中引數2的意思是將a矩陣a中的按“行”為單位進行求和
    YY = sum(Y'.*Y',2);
    XY = X'*Y;
    K = abs(repmat(XX,[1 size(YY,1)]) + repmat(YY',[size(XX,1) 1]) - 2*XY);
    K = exp(-K./delta);
end

測試函式如下

function result = svmTest(svm, Xt, Yt, kertype)
temp = (svm.a'.*svm.Ysv)*kernel(svm.Xsv,svm.Xsv,kertype);
total_b = svm.Ysv-temp;
b = mean(total_b);
w = (svm.a'.*svm.Ysv)*kernel(svm.Xsv,Xt,kertype);
result.score = w + b;
Y = sign(w+b);
result.Y = Y;
result.accuracy = size(find(Y==Yt))/size(Yt);

要說明的是，上面的程式碼是實現的關於我們樣本只有2個屬性的情況；當我們的樣本用多個屬性時，我們需要修改部分程式碼即可。