svmmatlab4.txt

％ 支持向量机Matlab工具箱1.0 - One-Class SVM, 一类支持向量机 ％ 使用平台 - Matlab6.5 希望对大家有用

% 支持向量机Matlab工具箱1.0 - One-Class SVM, 一类支持向量机
% 使用平台 - Matlab6.5
% 版权所有：陆振波，海军工程大学
% 电子邮件：[e-mail]
% 个人主页：http://luzhenbo.88uu.com.cn
% 参数文献：Chih-Chung Chang, Chih-Jen Lin. "LIBSVM: a Library for Support Vector Machines"
%
% Support Vector Machine Matlab Toolbox 1.0 - One-Class Support Vector Machine
% Platform : Matlab6.5 / Matlab7.0
% Copyright : LU Zhen-bo, Navy Engineering University, WuHan, HuBei, P.R.China, 430033 
% E-mail : [e-mail] 
% Homepage : http://luzhenbo.88uu.com.cn 
% Reference : Chih-Chung Chang, Chih-Jen Lin. "LIBSVM: a Library for Support Vector Machines"
%
% Solve the quadratic programming problem - "quadprog.m"

clc
clear
close all

% ------------------------------------------------------------%
% 定义核函数及相关参数

nu = 0.15; % nu -> [0,1] 在支持向量数与错分样本数之间进行折衷
% 支持向量机的 nu 参数(取值越小，异常点就越少)

ker = struct('type','linear');
%ker = struct('type','ploy','degree',3,'offset',1);
%ker = struct('type','gauss','width',200);
%ker = struct('type','tanh','gamma',1,'offset',0);

% ker - 核参数(结构体变量)
% the following fields:
% 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)
% degree - Degree d of polynomial kernel (positive scalar).
% offset - Offset c of polynomial and tanh kernel (scalar, negative for tanh).
% width - Width s of Gauss kernel (positive scalar).
% gamma - Slope g of the tanh kernel (positive scalar).

% ------------------------------------------------------------%
% 构造一类训练样本

n = 100
randn('state',1);
x1 = randn(floor(n*0.95),2);
x2 = 4+randn(ceil(n*0.05),2);
X = [x1;x2]; % 训练样本,n×d的矩阵,n为样本个数,d为样本维数

figure(3);
plot(x1(:,1),x1(:,2),'bx',x2(:,1),x2(:,2),'k.');
axis([-5 8 -5 8]); 
hold on;

% ------------------------------------------------------------%
% 训练支持向量机

tic
svm = One_Class_SVM_Train(X,[],nu,ker);
t_train = toc

% svm 支持向量机(结构体变量)
% the following fields:
% ker - 核参数
% x - 训练样本
% y - 训练目标;
% a - 拉格朗日乘子

% ------------------------------------------------------------%
% 寻找支持向量

a = svm.a;
epsilon = 1e-10; % 如果小于此值则认为是0
i_sv = find(a>epsilon); % 支持向量下标
plot(X(i_sv,1),X(i_sv,2),'ro');

% ------------------------------------------------------------%
% 测试输出

[x1,x2] = meshgrid(-4:0.05:6,-4:0.05:6);
[rows,cols] = size(x1);
nt = rows*cols; % 测试样本数
Xt = [reshape(x1,nt,1),reshape(x2,nt,1)];

tic
Yd = One_Class_SVM_Sim(svm,Xt); % 测试输出
t_sim = toc

Yd = reshape(Yd,rows,cols);
contour(x1,x2,Yd,[0 0],'m'); % 分类面
hold off;

function [K] = CalcKernel(ker,x,y)
% Calculate kernel function. 
%
% x: 输入样本,n1×d的矩阵,n1为样本个数,d为样本维数
% y: 输入样本,n2×d的矩阵,n2为样本个数,d为样本维数
%
% ker 核参数(结构体变量)
% the following fields:
% 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)
% degree - Degree d of polynomial kernel (positive scalar).
% offset - Offset c of polynomial and tanh kernel (scalar, negative for tanh).
% width - Width s of Gauss kernel (positive scalar).
% gamma - Slope g of the tanh kernel (positive scalar).
%
% ker = struct('type','linear');
% ker = struct('type','ploy','degree',d,'offset',c);
% ker = struct('type','gauss','width',s);
% ker = struct('type','tanh','gamma',g,'offset',c);
%
% K: 输出核参数,n1×n2的矩阵

%-------------------------------------------------------------%
% 转成列向量
x = x';
y = y';

%-------------------------------------------------------------%

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 svm = One_Class_SVM_Train(X,tmp,nu,ker)

% 输入参数:
% X 训练样本,n×d的矩阵,n为样本个数,d为样本维数
% nu 控制参数
% ker 核参数(结构体变量)
% the following fields:
% 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)
% degree - Degree d of polynomial kernel (positive scalar).
% offset - Offset c of polynomial and tanh kernel (scalar, negative for tanh).
% width - Width s of Gauss kernel (positive scalar).
% gamma - Slope g of the tanh kernel (positive scalar).

% 输出参数:
% svm 支持向量机(结构体变量)
% the following fields:
% ker - 核参数
% x - 训练样本
% y - 训练目标;
% a - 拉格朗日乘子

% ------------------------------------------------------------%
% 解二次优化

n = size(X,1);
H = Calckernel(ker,X,X);

f = zeros(n,1);
for i = 1:n
f(i, = -Calckernel(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);

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

% ------------------------------------------------------------%
% 输出 svm
svm.ker = ker;
svm.x = X;
svm.y = [];
svm.a = a;

function Yd = One_Class_SVM_Sim(svm,Xt)
% 输入参数:
% svm 支持向量机(结构体变量)
% the following fields:
% ker - 核参数
% 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)
% degree - Degree d of polynomial kernel (positive scalar).
% offset - Offset c of polynomial and tanh kernel (scalar, negative for tanh).
% width - Width s of Gauss kernel (positive scalar).
% gamma - Slope g of the tanh kernel (positive scalar).
% x - 训练样本
% y - 训练目标;
% a - 拉格朗日乘子
%
% Xt 测试样本,n×d的矩阵,n为样本个数,d为样本维数

% 输出参数:
% Yd 测试输出,n×1的矩阵,n为样本个数,值为+1或-1

% ------------------------------------------------------------%

ker = svm.ker;
X = svm.x;
a = svm.a;

% ------------------------------------------------------------%
% 求超球半径 R

epsilon = 1e-8; % 如果小于此值则认为是0
i_sv = find(a>epsilon); % 支持向量下标
n_sv = length(i_sv);

tmp1 = zeros(n_sv,1);
for i = 1:n_sv
index = i_sv(i);
tmp1(i) = Calckernel(ker,X(index,,X(index,);
end

tmp2 = 2*Calckernel(ker,X(i_sv,,X)*a; % 列向量
tmp3 = sum(sum(a*a'.*Calckernel(ker,X,X))); 

R_square = tmp1-tmp2+tmp3;
R = mean(sqrt(R_square)); % 超球半径 R (对所有支持向量求平均的结果)

% ------------------------------------------------------------%
% 测试输出

nt = size(Xt,1); % 测试样本数

tmp4 = zeros(nt,1); % 列向量
for i = 1:nt
tmp4(i) = Calckernel(ker,Xt(i,,Xt(i,);
end

tmp5 = 2*Calckernel(ker,Xt,X)*a; % 列向量

Yd = sign(tmp4-tmp5+tmp3-R^2);