点击上方蓝字关注“公众号”
基于matlab的支持向量机分类、回归问题
Part.1
支持向量机(support vector machines)是一种二分类模型，它的目的是寻找一个超平面来对样本进行分割，分割的原则是间隔最大化，最终转化为一个凸二次规划问题来求解。
对于线性可分来说程序一大筐，今天我们主要来讲非线性问题。主要有二分类、一分类和回归问题的处理。
对于这样的问题，可以将训练样本从原始空间映射到一个更高维的空间，使得样本在这个空间中线性可分，如果原始空间维数是有限的，即属性是有限的，那么一定存在一个高维特征空间使样本可分。
令ϕ(x)表示将 x 映射后的特征向量，于是在特征空间中，划分超平面所对应的的模型可表示为：
f(x)=wTϕ(x)+b
于是有最小化函数：
在实际应用中，通常人们会从一些常用的核函数里选择(根据样本数据的不同，选择不同的参数，实际上就得到了不同的核函数)。
Part.2
下面介绍每一个程序的思想。
核函数
首先是编写核函数：
function [K] = kernel(ker,x,y)
% x: 输入样本,d×n1的矩阵,n1为样本个数,d为样本维数
% y: 输入样本,d×n2的矩阵,n2为样本个数,d为样本维数
% ker 核参数(结构体变量)
% K: 输出核参数,n1×n2的矩阵
switch ker.type
case 'linear'
K = x'*y;
case 'ploy'
d = ker.degree;
c = ker.offset;
K = (x'*y+c).^d;
case 'gauss'
s = ker.width;
rows = size(x,2);
cols = size(y,2);
tmp = zeros(rows,cols);
for i = 1:rows
for j = 1:cols
tmp(i,j) = norm(x(:,i)-y(:,j));
end
end
K = exp(-0.5*(tmp/s).^2);
case 'tanh'
g = ker.gamma;
c = ker.offset;
K = tanh(g*x'*y+c);
otherwise
K = 0;
end
决策函数
function Yd = svmSim(svm,Xt)
cathe = 10e+6; % kernel输出的元素个数的上限
nx = size(svm.x,2); % 训练样本数
nt = size(Xt,2); % 测试样本数
block = ceil(nx*nt/cathe); % 分块处理
num = ceil(nt/block); % 每块测试样本数
for i = 1:block
if (i==block)
index = [(i-1)*num+1:nt];
else
index = (i-1)*num+[1:num];
end
Yd(index) = svmSim_block(svm,Xt(:,index));
% 测试输出
end
function Yd = svmSim_block(svm,Xt);
% x - 训练样本
% y - 训练目标;
% a - 拉格朗日乘子
% Xt 测试样本,d×n的矩阵,n为样本个数,d为样本维数
% 输出参数:
% Yd 测试输出,1×n的矩阵,n为样本个数,值为+1或-1
type = svm.type;
ker = svm.ker;
X = svm.x;
Y = svm.y;
a = svm.a;
% 测试输出
epsilon = 1e-8; % 如果小于此值则认为是0
i_sv = find(abs(a)>epsilon); % 支持向量下标
switch type
case 'svc_c',
tmp = (a.*Y)*kernel(ker,X,X(:,i_sv)); % 行向量
b = Y(i_sv)-tmp;
b = mean(b);
tmp = (a.*Y)*kernel(ker,X,Xt);
tmp = tmp+b;
Yd = sign(tmp);
case 'svc_nu',
% 与 'svc_c' 情况相同
tmp = (a.*Y)*kernel(ker,X,X(:,i_sv)); % 行向量
b = Y(i_sv)-tmp;
b = mean(b);
tmp = (a.*Y)*kernel(ker,X,Xt);
Yd = sign(tmp+b);
case 'svm_one_class',
n_sv = length(i_sv);
tmp1 = zeros(n_sv,1);
for i = 1:n_sv
index = i_sv(i);
tmp1(i) = kernel(ker,X(:,index),X(:,index));
end
tmp2 = 2*a*kernel(ker,X,X(:,i_sv)); % 行向量
tmp3 = sum(sum(a'*a.*kernel(ker,X,X)));
R_square = tmp1-tmp2'+tmp3;
R_square = mean(R_square); % 超球半径 R^2 (对所有支持向量求平均的结果)
nt = size(Xt,2); % 测试样本数
tmp4 = zeros(nt,1); % 列向量
for i = 1:nt
tmp4(i) = kernel(ker,Xt(:,i),Xt(:,i));
end
tmp5 = 2*a*kernel(ker,X,Xt); % 行向量
Yd = sign(tmp4-tmp5'+tmp3-R_square);
case 'svr_epsilon',
tmp = a*kernel(ker,X,X(:,i_sv)); % 行向量
b = Y(i_sv)-tmp;
% 符号不一样,决策函数就不一样,实际上是一回事!
%b = Y(i_sv)+tmp;
b = mean(b);
tmp = a*kernel(ker,X,Xt);
% 符号不一样,决策函数就不一样
%tmp = -a*kernel(ker,X,Xt);
Yd = (tmp+b);
case 'svr_nu',
% 与'svr_epsilon' 情况相同
tmp = a*kernel(ker,X,X(:,i_sv)); % 行向量
b = Y(i_sv)-tmp;
% 符号不一样,决策函数就不一样
%b = Y(i_sv)+tmp;
b = mean(b);
tmp = a*kernel(ker,X,Xt);
% 符号不一样,决策函数就不一样
%tmp = -a*kernel(ker,X,Xt);
Yd = (tmp+b);
otherwise,
end
训练函数
function svm = svmTrain(svmType,X,Y,ker,p1,p2)
% 输入参数:
% X 训练样本,d×n的矩阵,n为样本个数,d为样本维数
% Y 训练目标,1×n的矩阵,n为样本个数,值为+1或-1
% ker 核参数(结构体变量)
% 输出参数:
% svm 支持向量机(结构体变量)
% type - 支持向量机类型 {'svc_c','svc_nu','svm_one_class','svr_epsilon','svr_nu'}
% a - 拉格朗日乘子,1×n的矩阵
options = optimset;
options.LargeScale = 'off';
options.Display = 'off';
switch svmType
case 'svc_c',
C = p1;
n = length(Y);
H = (Y'*Y).*kernel(ker,X,X);
f = -ones(n,1);
A = [];
b = [];
Aeq = Y;
beq = 0;
lb = zeros(n,1);
ub = C*ones(n,1);
a0 = zeros(n,1);
[a,fval,eXitflag,output,lambda] = quadprog(H,f,A,b,Aeq,beq,lb,ub,a0,options);
case 'svc_nu',
nu = p1;
n = length(Y);
H = (Y'*Y).*kernel(ker,X,X);
f = zeros(n,1);
A = -ones(1,n);
b = -nu;
Aeq = Y;
beq = 0;
lb = zeros(n,1);
ub = ones(n,1)/n;
a0 = zeros(n,1);
[a,fval,eXitflag,output,lambda] = quadprog(H,f,A,b,Aeq,beq,lb,ub,a0,options);
case 'svm_one_class',
nu = p1;
n = size(X,2);
H = kernel(ker,X,X);
f = zeros(n,1);
for i = 1:n
f(i,:) = -kernel(ker,X(:,i),X(:,i));
end
A = [];
b = [];
Aeq = ones(1,n);
beq = 1;
lb = zeros(n,1);
ub = ones(n,1)/(nu*n);
a0 = zeros(n,1);
[a,fval,eXitflag,output,lambda] = quadprog(H,f,A,b,Aeq,beq,lb,ub,a0,options);
case 'svr_epsilon',
C = p1;
e = p2;
n = length(Y);
Q = kernel(ker,X,X);
H = [Q,-Q;-Q,Q];
f = [e*ones(n,1)-Y';e*ones(n,1)+Y'];
% 符号不一样,决策函数就不一样
%f = [e*ones(n,1)+Y';e*ones(n,1)-Y'];
A = [];
b = [];
Aeq = [ones(1,n),-ones(1,n)];
beq = 0;
lb = zeros(2*n,1);
ub = C*ones(2*n,1);
a0 = zeros(2*n,1);
[a,fval,eXitflag,output,lambda] = quadprog(H,f,A,b,Aeq,beq,lb,ub,a0,options);
a = a(1:n)-a(n+1:end);
case 'svr_nu',
C = p1;
nu = p2;
n = length(Y);
Q = kernel(ker,X,X);
H = [Q,-Q;-Q,Q];
f = [-Y';+Y'];
% 符号不一样,决策函数就不一样
%f = [+Y';-Y'];
A = [];
b = [];
Aeq = [ones(1,n),-ones(1,n);ones(1,2*n)];
beq = [0;C*n*nu];
lb = zeros(2*n,1);
ub = C*ones(2*n,1);
a0 = zeros(2*n,1);
[a,fval,eXitflag,output,lambda] = quadprog(H,f,A,b,Aeq,beq,lb,ub,a0,options);
a = a(1:n)-a(n+1:end);
otherwise,
end
eXitflag
% 输出 svm
svm.type = svmType;
svm.ker = ker;
svm.x = X;
svm.y = Y;
svm.a = a';
Part.3
主程序运行代码。
C二分类
% C-SVC, C二类分类算法
%clc
clear
%close all
% 定义核函数及相关参数
C = 200; % 拉格朗日乘子上界
ker = struct('type','linear');
%ker = struct('type','ploy','degree',3,'offset',1);
%ker = struct('type','gauss','width',1);
%ker = struct('type','tanh','gamma',1,'offset',0);
% 构造两类训练样本
n = 50;
randn('state',6);
x1 = randn(2,n);
y1 = ones(1,n);
x2 = 5+randn(2,n);
y2 = -ones(1,n);
figure(1);
plot(x1(1,:),x1(2,:),'bx',x2(1,:),x2(2,:),'k.');
axis([-3 8 -3 8]);
title('C-SVC')
hold on;
X = [x1,x2];
% 训练样本,d×n的矩阵,n为样本个数,d为样本维数
Y = [y1,y2];
% 训练目标,1×n的矩阵,n为样本个数,值为+1或-1
% 训练支持向量机
tic
svm = svmTrain('svc_c',X,Y,ker,C);
t_train = toc
% 寻找支持向量
a = svm.a;
epsilon = 1e-8;
% 如果小于此值则认为是0
i_sv = find(abs(a)>epsilon);
% 支持向量下标
plot(X(1,i_sv),X(2,i_sv),'ro');
% 测试输出
[x1,x2] = meshgrid(-2:0.1:7,-2:0.1:7);
[rows,cols] = size(x1);
nt = rows*cols;
% 测试样本数
Xt = [reshape(x1,1,nt);reshape(x2,1,nt)];
tic
Yd = svmSim(svm,Xt);
% 测试输出
t_sim = toc
Yd = reshape(Yd,rows,cols);
contour(x1,x2,Yd,[0 0],'m');
% 分类面
hold off;
Nu二分类
% Nu-SVC, Nu二类分类算法
clc
clear
close all
nu = 0.2; % nu -> (0,1] 在支持向量数与错分样本数之间进行折衷
ker = struct('type','linear');
%ker = struct('type','ploy','degree',3,'offset',1);
%ker = struct('type','gauss','width',1);
%ker = struct('type','tanh','gamma',1,'offset',0);
% 构造两类训练样本
n = 50;
randn('state',6);
x1 = randn(2,n);
y1 = ones(1,n);
x2 = 5+randn(2,n);
y2 = -ones(1,n);
figure;
plot(x1(1,:),x1(2,:),'bx',x2(1,:),x2(2,:),'k.');
axis([-3 8 -3 8]);
title('C-SVC')
hold on;
X = [x1,x2];
% 训练样本,d×n的矩阵,n为样本个数,d为样本维数
Y = [y1,y2];
% 训练目标,1×n的矩阵,n为样本个数,值为+1或-1
% 训练支持向量机
tic
svm = svmTrain('svc_nu',X,Y,ker,nu);
t_train = toc
% 寻找支持向量
a = svm.a;
epsilon = 1e-8; % 如果小于此值则认为是0
i_sv = find(abs(a)>epsilon); % 支持向量下标
plot(X(1,i_sv),X(2,i_sv),'ro');
% 测试输出
[x1,x2] = meshgrid(-2:0.1:7,-2:0.1:7);
[rows,cols] = size(x1);
nt = rows*cols; % 测试样本数
Xt = [reshape(x1,1,nt);reshape(x2,1,nt)];
tic
Yd = svmSim(svm,Xt); % 测试输出
t_sim = toc
Yd = reshape(Yd,rows,cols);
contour(x1,x2,Yd,[0 0],'m'); % 分类面
hold off;
一种一类
%One-Class SVM, 一类支持向量机
clc
clear
close all
% 定义核函数及相关参数
nu = 0.2; % nu -> [0,1] 在支持向量数与错分样本数之间进行折衷
% 支持向量机的 nu 参数(取值越小，异常点就越少)
ker = struct('type','linear');
% type - linear : k(x,y) = x'*y
% poly : k(x,y) = (x'*y+c)^d
% gauss : k(x,y) = exp(-0.5*(norm(x-y)/s)^2)
% tanh : k(x,y) = tanh(g*x'*y+c)
% 构造一类训练样本
n = 100
randn('state',1);
x1 = randn(2,floor(n*0.95));
x2 = 5+randn(2,ceil(n*0.05));
X = [x1,x2]; % 训练样本,d×n的矩阵,n为样本个数,d为样本维数
figure;
plot(x1(1,:),x1(2,:),'bx',x2(1,:),x2(2,:),'k.');
axis([-5 8 -5 8]);
title('One-Class SVM');
hold on;
% 训练支持向量机
tic
svm = svmTrain('svm_one_class',X,[],ker,nu);
t_train = toc
% 寻找支持向量
a = svm.a;
epsilon = 1e-10; % 如果小于此值则认为是0
i_sv = find(abs(a)>epsilon); % 支持向量下标
plot(X(1,i_sv),X(2,i_sv),'ro');
% 测试输出
[x1,x2] = meshgrid(-5:0.05:6,-5:0.05:6);
[rows,cols] = size(x1);
nt = rows*cols; % 测试样本数
Xt = [reshape(x1,1,nt);reshape(x2,1,nt)];
tic
Yd = svmSim(svm,Xt); % 测试输出
t_sim = toc
Yd = reshape(Yd,rows,cols);
contour(x1,x2,Yd,[0 0],'m'); % 分类面
hold off;
Epsilon回归算法
% Epsilon-SVR, Epsilon回归算法
clc
clear
%close all
% 定义核函数及相关参数
C = 100; % 拉格朗日乘子上界
e = 0.2; % 不敏感损失函数的参数，Epsilon越大，支持向量越少
% 构造两类训练样本
n = 50;
rand('state',42);
X = linspace(-4,4,n); % 训练样本,d×n的矩阵,n为样本个数,d为样本维数,这里d=1
Ys = (1-X+2*X.^2).*exp(-.5*X.^2);
f = 0.2; % 相对误差
Y = Ys+f*max(abs(Ys))*(2*rand(size(Ys))-1)/2; % 训练目标,1×n的矩阵,n为样本个数,值为期望输出
figure;
plot(X,Ys,'b-',X,Y,'b*');
title('\epsilon-SVR');
hold on;
% 训练支持向量机
tic
svm = svmTrain('svr_epsilon',X,Y,ker,C,e);
t_train = toc
% 寻找支持向量
a = svm.a;
epsilon = 1e-8; % 如果"绝对值"小于此值则认为是0
i_sv = find(abs(a)>epsilon); % 支持向量下标,这里对abs(a)进行判定
plot(X(i_sv),Y(i_sv),'ro');
% 测试输出
tic
Yd = svmSim(svm,X); % 测试输出
t_sim = toc
plot(X,Yd,'r--');
hold off;
Nu回归算法
% Nu-SVR, Nu回归算法
clc
clear
%close all
% 定义核函数及相关参数
C = 100; % 拉格朗日乘子上界
nu = 0.05; % nu -> (0,1] 在支持向量数与拟合精度之间进行折衷
% 构造两类训练样本
n = 50;
rand('state',42);
X = linspace(-4,4,n); % 训练样本,d×n的矩阵,n为样本个数,d为样本维数,这里d=1
Ys = (1-X+2*X.^2).*exp(-.5*X.^2);
f = 0.2; % 相对误差
Y = Ys+f*max(abs(Ys))*(2*rand(size(Ys))-1)/2; % 训练目标,1×n的矩阵,n为样本个数,值为期望输出
figure;
plot(X,Ys,'b-',X,Y,'b*');
title('\nu-SVR');
hold on;
% 训练支持向量机
tic
svm = svmTrain('svr_nu',X,Y,ker,C,nu);
t_train = toc
% 寻找支持向量
a = svm.a;
epsilon = 1e-8; % 如果"绝对值"小于此值则认为是0
i_sv = find(abs(a)>epsilon); % 支持向量下标,这里对abs(a)进行判定
plot(X(i_sv),Y(i_sv),'ro');
% 测试输出
tic
Yd = svmSim(svm,X); % 测试输出
t_sim = toc
plot(X,Yd,'r--');
hold off;
Part.4
仿真结果展示。
C二分类
Nu二分类
一种一类
Epsilon回归算法
Nu回归算法
详细程序工具箱及资料，公众号回复【SVM】可获得。
扫码关注
不迷路
1368
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
