对“视觉机器学习20讲配套仿真代码”的研究心得---SVM方法

来源：互联网发布：淘宝云客服考试编辑：程序博客网时间：2024/06/07 01:56

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%功能：演示SVM算法在计算机视觉中的应用
%基于SVM实现特征分类；
%环境：Win7，Matlab2012b
%Modi: NUDT-VAP
%时间：2013-09-23
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% kernel function
C = 200;
ker = struct('type','linear');

&&&&&&&&&&&&&&&

struct在matlab中是用来建立结构体数组的。通常有两种用法：
s = struct('field1',{},'field2',{},...)
这是建立一个空的结构体，field1，field2是域名称
s = struct('field1',values1,'field2',values2,...)
这是在建立结构体的时候就赋值，values1,values2就是针对field1,field2所赋的值。
例如：s = struct('type',{'big','little'},'color',{'red'},'x',{3 4})
这个结构体有三个域'type','color','x'，根据英文意义可以看出分别是类型，颜色和值的意思。类型可以取的值有'big','little'，也就是类型有两种，分别是大和小。颜色只有'red'红色，而值x可以是3或者4
matlab帮助里解释的很详细，你自己可以去check一下

&&&&&&&&&&&&&&&&&

% test sample----------------------------------------------------------测试案例
n = 50;
randn('state',6);
x1 = randn(2,n);
y1 = ones(1,n);
x2 = 5+randn(2,n);
y2 = -ones(1,n);

&&&&&&&&&&&&&&&&&&&&&&&&&

randn：产生正态分布的随机数或矩阵的函数

randn：产生均值为0，方差σ^2 = 1，标准差σ = 1的正态分布的随机数或矩阵的函数。

用法：

Y = randn(n)：返回一个n*n的随机项的矩阵。如果n不是个数量，将返回错误信息。

Y = randn(m,n) 或 Y = randn([m n])：返回一个m*n的随机项矩阵。

Y = randn(m,n,p,...) 或 Y = randn([m n p...])：产生随机数组。

Y = randn(size(A))：返回一个和A有同样维数大小的随机数组。

randn

返回一个每次都变化的数量。

s = randn('state')

举例

Example 1. R = randn(3,4) 将生成矩阵

R =

1.1650 0.3516 0.0591 0.8717

0.6268 -0.6965 1.7971 -1.4462

0.0751 1.6961 0.2641 -0.7012

For a histogram of the randn distribution, see hist.

Example 2. 产生一个随机分布的指定均值和方差的矩阵：将randn产生的结果乘以标准差，然后加上期望均值即可。例如，产生均值为0.6，方差为0.1的一个5*5的随机数方式如下：

x = .6 + sqrt(0.1) * randn(5)

x =

0.8713 0.4735 0.8114 0.0927 0.7672

0.9966 0.8182 0.9766 0.6814 0.6694

0.0960 0.8579 0.2197 0.2659 0.3085

0.1443 0.8251 0.5937 1.0475 -0.0864

0.7806 1.0080 0.5504 0.3454 0.5813

其他类似函数：rand, randperm, sprand, sprandn

均匀分布的随机数或矩阵
语法
Y = rand(n)
Y = rand(m,n)
Y = rand([m n])
Y = rand(m,n,p,...)
Y = rand([m n p...])
Y = rand(size(A))
rand
s = rand('state')
描述
rand函数产生由在(0, 1)之间均匀分布的随机数组成的数组
Y = rand(n) 返回一个n x n的随机矩阵如果n不是数量，则返回错误信息
Y = rand(m,n) 或 Y = rand([m n]) 返回一个m x n的随机矩阵
Y = rand(m,n,p,...) 或Y = rand([m n p...]) 产生随机数组
Y = rand(size(A)) 返回一个和A有相同尺寸的随机矩阵
1，rand(3)*-2 rand（3）是一个3*3的随机矩阵（数值范围在0~1之间）然后就是每个数乘上-2
2 ，用matlab随机产生60个1到365之间的正数1+fix（365*rand（1，60））；
3，用rand函数随机取100个从-1到1的数x1，x2，...，x = rand(1,100) * 2 - 1

&&&&&&&&&&&&&&&&&&&&&&&&&&

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];
Y = [y1,y2];

&&&&&&&&&&&&&&&&&&&&&&&&&&&

matlab中的 figure 命令，能够创建一个用来显示图形输出的一个窗口对象。每一个这样的窗口都有一些属性，例如窗口的尺寸、位置，等等。下面一一介绍它们。

一、概述

总的来说，figure 的使用语法包括：

figure

figure('PropertyName',propertyvalue,...)

figure(h)

h = figure(...)

第一种用法最简单，它创建一个窗口，其各种属性都是使用默认设置。例如它创建的窗口立即成为当前窗口，并显示在其它窗口之上。直到新的窗口被创建或者其它窗口被唤醒（called）。

第二种用法，figure('PropertyName',propertyvalue,...)，则可以指定某些属性。例如“Name”属性，则可以指定该窗口的标题：

第三种用法，figure(h)，则根据参数 h 的不同而表现不同。如果 h 是先前一个窗口的句柄（handle），则figure(h)相当于唤醒了该窗口，使得该窗口为当前窗口；如果 h 并不是某个窗口的句柄，但它是一个整数，则 figure(h)
创建一个句柄为 h 的新窗口。

最后一种用法，h
= figure(...)，相当于前面三种用法中在创建一个窗口的同时，返回该窗口的句柄。

二、用法举例

1、在同一个figure中显示多个contour，可通过“hold on” 来实现：

a = ones(500,500);

figure('name','叠加显示图像');

b = a;

b(10:500-10,10:500-10) = 0;

contour(b,[0 0],'r');
hold on;

c = a;

c(100:400,100:400) = 0;

contour(c,[0 0],'b');

hold off;

2、figure 默认的坐标原点在左下角，这会导致有些图形图像的显示出现上下翻转，这时可通过：set(gca,'ydir','reverse','xaxislocation','top'); % 将坐标原点设在左上角

3、指定坐标轴显示的坐标范围：xlim([xmin xmax])

功能

二维曲线绘图

语法

plot(Y)
plot(X1,Y1,...)
plot(X1,Y1,LineSpec,...)
plot(...,'PropertyName',PropertyValue,...)
plot(axes_handle,...)
h = plot(...)
hlines = plot('v6',...)

描述

plot(Y)如果Y是m×n的数组，以1:m为X横坐标，Y中的每一列元素为Y坐标，绘制n条曲线；如果Y是n×1或者1×n的向量，则以1:n为横坐标，Y为坐标表绘制1条曲线；如果Y是复数，则plot(Y)等效于plot(real(Y),imag(Y))；其它使用情况下，忽略坐标数据中的虚部。

plot(X1,Y1,...)如果X和Y都是数组，按列取坐标数据绘图，此时它们必须具有相同的尺寸；如果X和Y其中一个是向量另一个为数组，X和Y中尺寸相等的方向对应绘制多条曲线；如果X和Y其中一个是标量另一个为向量，那么将绘制垂直X或者Y轴离散的点。

备注

m×n

按列取坐标数据绘制n条曲线

X和Y必须具有相同的尺寸

1×n或n×1

m×n或n×m

自动匹配尺寸相同方向绘制m条曲线

任意四种组合，效果一样

m×n或n×m

1×n或n×1

同上

1×1

1×n或n×1

绘制垂直X轴离散点

Y可以任意向量

1×n或n×1

1×1

绘制垂直Y轴离散点

X可以任意向量

plot(X1,Y1,LineSpec,...)通过参数LineSpec指定曲线的曲线属性，它包括线型、标记符和颜色。plot函数支持同时绘制任意组图形

plot(X1,Y1,LineSpec1,X2,Y2,LineSpec2,...)

此时完全等效于

plot(X1,Y1,LineSpec1,...)

hlod all

plot(X2,Y2,LineSpec2,...)

MATLAB中提供的线型属性有：

线型

说明

标记符

说明

颜色

说明

实线(默认)

加号符

红色

双划线

空心圆

绿色

虚线

星号

蓝色

点划线

实心圆

青绿色

叉号符

洋红色

正方形

黄色

菱形

黑色

上三角形

白色

下三角形

右三角形

左三角形

五角星

六边形

需要说明的是，LineSpec中设置曲线线型、标识符和颜色三项属性时，控制符的顺序不受限制并可以省略或者部分省略。也就是说'r-.*'、'-.r*'、'*-.r'等形式是等效的，都表示使用红色点划线连接各个节点，各节点使用“*”标识。

plot(...,'PropertyName',PropertyValue,...)设置由plot创建的所有曲线句柄对象的属性，Line对象属性和属性值参见附录，具体设置参考下面的实例，当然可以使用set/get进行设置。

plot(axes_handle,...)指定坐标系，也就是在axes_handle坐标系中绘图，在没有指定时默认为gca。

h = plot(...)返回由plot创建的所有曲线句柄对象的句柄。每条曲线对应一个句柄，如果有n条曲线，则h为n×1的数组。

注意

在同时绘制多条曲线时，如果没有指定曲线属性，plot按顺序循环使用当前坐标系中ColorOrder和LineStyleOrder两个属性。

默认情况，MATLAB在每次调用plot函数时将ColorOrder和LineStyleOrder自动重置为DefaultAxesColorOrder和DefaultAxesLineStyleOrder。Default**属性我们可以自定义，有效期至MATLAB关闭，Matlab下次启动时将Default**属性重置为厂家设置(Factory)

set(0,'DefaultAxesColorOrder',’r|g|b|k’,...

      'DefaultAxesLineStyleOrder','-|-.|--|:')

使用hold all命令可以阻止调用plot函数时自动重置ColorOrder和LineStyleOrder属性，而是循环使用。注意hold on只是使多次绘制的图形叠加（相当于NextPlot），但不能阻止属性重置。

另外我们可以通过下面四个属性设置标识符的颜色和大小

LineWidth——指定线宽

MarkerEdgeColor——指定标识符的边缘颜色

MarkerFaceColor——指定标识符填充颜色

MarkerSize——指定标识符的大小

注意上面四个属性是针对当前坐标系中所有曲线的

实例

% by dynamic

% see also

% 2009.8.20

X=1:10;

% 两个都是数组，必须具有相同的尺寸

X1=[X;X;X]';%10×3

Y1=rand(10,3)+1;%10×3

% 其中一个为向量，另一个为数组，自动匹配尺寸相等方向

X2=1:0.1:10;%1×91

Y2=[sin(X2);cos(X2)]';%91×2

% 其中一个是标量，另一为矢量，绘制垂直坐标轴的离散点

X3=1:10;

Y3=-0.5;

fh=figure('numbertitle','off','name','PLOT Usability Demo');%创建figure对象

ah=axes;%创建axes对象

h=plot(...%返回所有曲线句柄

ah,...%指定坐标系，可以省略，此时默认gca

X1,Y1,...%坐标数据

'-.^',...%曲线属性，可以省略或部分省略，此时自动选择

X2,Y2,...

'm-',...

X3,Y3,...

'o',...%注意此组数据设置线型和颜色无效，因为默认绘制离散点

'LineWidth',2,...%线宽

'MarkerEdgeColor','k',...%标识符边缘颜色

'MarkerFaceColor','r',...%标识符填充颜色

'MarkerSize',8)%标识符大小

【原创】Matlab中plot函数全功能解析

在plot写完后 加上axis([0 1 2 3])则表明图线的x轴范围为0~1y轴范围为2~3如   plot(x,y)     axis([0 1 2 3])     grid on

显示希腊字符，matlab是支持latex语法的，所以希腊字母可以这样表示：title(['\alpha ','\beta ','\gamma ','\theta ','\delta ','\Gamma ','\Theta ','\Delta '],'Fontsize',15);
增加上标下标：title(['y=x^{我是上标}_{我是下标}',' \alpha'],'Fontsize',20);
将程序中变量值显示在标题中：
x=100;
title(['x=',num2str(x)],'Fontsize',20);
title颜色：
title('我是蓝色','Fontsize',20,'Color','b');
更改title位置：title(['我在（50，0.7）'],'position',[50,0.6],'Fontsize',15);
更改title字号，更改title字体：
1、title('字体黑体，字号20','Fontname','黑体','Fontsize',20);
2、title(['\fontname{宋体}\fontsize{15}{我是宋体15号}','\fontname{黑体}\fontsize{20}{我是黑体20号}'...
sprintf('\n'),'\fontname{Times new roman}{I am English}']);
title换行(多行显示)：title(['我是第一行',sprintf('\n'),'我是第二行'],'Fontsize',15);

上面的几行代码是整个SVM界面的核心。

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

% train SVM--------------------------- 训练SVM，这里的svmTrain是一个自己编写的函数，在MATLAB中有相应的接口。
tic
svm = svmTrain('svc_c',X,Y,ker,C);----------------------------- 在接下来的关于svmTrain这个函数的讲解中，大家会了解这个函数的真实意图。
t_train = toc

% find sustain vector ----------------------------   找到支持向量
a = svm.a;
epsilon = 1e-8; ---------------------- 希腊语的第五个字母   1e-8 ==0.00000001

i_sv = find(abs(a)>epsilon); --------------这里的意思就是说，找出 abs(a)>epsilon的值。然后赋值给i_sv ，为我们下一步画图做准备。

plot(X(1,i_sv),X(2,i_sv),'ro');----------------

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

PLOT 后面的 ro 代表红色，圆环。

Various line types, plot symbols and colors may be obtained with
plot(X,Y,S) where S is a character string made from one element
from any or all the following 3 columns:

b blue . point - solid
g green o circle : dotted
r red x x-mark -. dashdot
c cyan + plus -- dashed
m magenta * star (none) no line
y yellow s square
k black d diamond
w white v triangle (down)
^ triangle (up)
< triangle (left)
> triangle (right)
p pentagram
h hexagram

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

% test output 测试输出

[x1,x2] = meshgrid(-2:0.1:7,-2:0.1:7);

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

语法

[X,Y] = meshgrid(x,y)

解释：输出X的每一行的数值都是复制的x的值；输出Y的每一列的数值都是复制的y的值。

[X,Y]=meshgrid(x)与[X,Y]=meshgrid(x,x)是等同的

[X,Y,Z]=meshgrid(x,y,z)生成三维数组，可用来计算三变量的函数和绘制三维立体图

相关函数： plot3、mesh、surf、automesh、ndgrid

程序示例

编辑

示例一：

x=-3:1:3;y=-2:1:2;

[X,Y]= meshgrid(x,y);

这里meshgrid（x，y）的作用是分别产生以向量x为行，向量y为列的两个大小相同的矩阵，其中x的行是从-3开始到3，每间隔1记下一个数据，并把这些数据集成矩阵X；同理y的列则是从-2到2，每间隔1记下一个数据，并集成矩阵Y。即

-3 -2 -1 0 1 2 3

Y =

-2 -2 -2 -2 -2 -2 -2

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

0 0 0 0 0 0 0

1 1 1 1 1 1 1

2 2 2 2 2 2 2

示例二：

function main

close all; clear; clc;

M1; M2;

end

function M1

x = rand(3, 4);

y = rand(2, 3);

size_of_x = size(x)

size_of_y = size(y)

[X, Y] = meshgrid(x, y);

size_of_X = size(X)

size_of_Y = size(Y)

end

function M2

x = rand(3, 4, 2);

y = rand(2, 3);

size_of_x = size(x)

size_of_y = size(y)

[X, Y] = meshgrid(x, y);

size_of_X = size(X)

size_of_Y = size(Y)

end

输出结果：

size_of_x =

3 4

size_of_y =

2 3

size_of_X =

6 12

size_of_Y =

6 12

size_of_x =

3 4 2

size_of_y =

2 3

size_of_X =

6 24

size_of_Y =

6 24

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

[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;

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

function Yd = svmSim(svm,Xt)

% parameters
cathe = 10e+6;
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);

type = svm.type;
ker = svm.ker;
X = svm.x;
Y = svm.y;
a = svm.a;

% test output
epsilon = 1e-8;
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',
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);

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',

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)
% SVM Classification:
% svm = svmTrain('svc_c',x,y,ker,C);
% svm = svmTrain('svc_nu',x,y,ker,nu);
%
% One-Class SVM:
% svm = svmTrain('svm_one_class',x,[],ker,nu);
%
% SVM Regression:
% svm = svmTrain('svr_epsilon',x,y,ker,C,e);
% svm = svmTrain('svr_nu',x,y,ker,C,nu);

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

% output svm

svm.type = svmType;
svm.ker = ker;
svm.x = X;
svm.y = Y;
svm.a = a';

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

function [K] = kernel(ker,x,y) 四种基本核：线性核、多项式核、径向基核和S形核，还有faruto提到的几种自定义核
% Calculate kernel function.
【模式识别】SVM核函数 http://blog.csdn.net/xiaowei_cqu/article/details/35993729

SVM入门（七）为何需要核函数

http://www.blogjava.net/zhenandaci/archive/2009/03/06/258288.html

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

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

0 0