实 验 目 的: 

1.实现支持向量机分类器，并利用半月形数据集进行分类任务。

2.训练模型并评估其性能，包括支持向量数量、分类准确率等指标。

3.对分类结果进行可视化分析，以便深入理解模型的分类能力。

实 验 原 理：

1.数据生成： 在实验中，首先通过生成函数生成一个半月形状的数据集。这个数据集包括两个类别的数据点，它们在空间中呈半月形分布。

2.支持向量机训练： 接下来，利用生成的数据集，训练一个支持向量机分类器。在训练过程中，使用了RBF核函数作为核函数，通过最大化分类间隔和最小化分类误差来找到最优的超平面。

3.模型评估： 在训练完成后，通过测试数据集对训练好的支持向量机模型进行评估。通过计算分类误差率和统计支持向量的数量等指标来评估模型的性能。

4.结果可视化： 最后，通过绘制训练结果和测试结果的图表，将训练数据、测试数据、决策边界等信息可视化展示出来，以直观地展示分类效果和模型性能。

实 验 设 备 与 器 件：

PC机、 MATLAB软件。

实 验 内 容 与 参 数 分 析

第一步

1. rad 和 width 分别定义了半月形的中心半径和宽度，这些参数影响了半月形的形状。
2. num_tr 和 num_te 分别定义了训练数据集和测试数据集的大小，即样本数量。
3. num_samp 定义了总的样本数量，是训练数据集和测试数据集的大小之和。
4. 通过 fprintf 函数打印输出了一些信息，包括数据生成的过程中使用的参数值，例如半月形数据集的生成、总数据点数、半月形的半径和宽度等。
5. 使用 rand('seed',seed) 设置了随机数生成器的种子，以确保每次生成的数据都是一致的，这里使用了 seed 变量作为种子值。
6. 调用 halfmoon() 函数生成半月形数据集，其中 halfmoon() 函数的输入参数包括半月形的中心半径、宽度、两个半月形之间的距离以及总的样本数量，输出包括原始数据 data 和随机打乱后的数据 data_shuffled。

使用 warning off 关闭了警告信息的显示，可能是因为在数据生成过程中可能会出现一些不必要的警告，此处关闭警告信息可能是为了提高代码的整洁性和可读性。第二步

参数定义：

1.C 是支持向量机的正则化参数。它的值越大，模型对错误分类的惩罚越严重，模型的复杂度越高。在这里，C 被设置为25，表示在训练模型时允许一定程度的错误分类。

2.p1 是RBF核函数的参数，用于控制核函数的宽度。RBF核函数的宽度决定了数据在特征空间中的分布范围。在这里，p1 被设置为1，这是一个常用的默认值。

3.epsilon 是支持向量检测的阈值。当拉格朗日乘子的绝对值小于该阈值时，被认为是支持向量。这个值通常很小，表示对支持向量的要求非常严格。

4.err 是一个计数器，用于记录分类错误的输出数量。在这里，初始化为0，之后会在模型测试过程中用于统计分类错误的数量。

5.b0 是隐式偏差，对于RBF核函数来说，通常设置为0。这是因为RBF核函数本身具有偏差项，所以不需要额外的偏差参数。

通过这些参数的设置，可以控制支持向量机模型的复杂度、核函数的宽度以及对错误分类的容忍程度，从而影响模型的性能和泛化能力。

第三步

这段代码的目的是对数据进行归一化处理，主要步骤包括：

1.计算数据的均值，并将数据以均值为中心化。

2.计算每个样本的前两行的绝对值，并找到最大值。

3.将数据除以每个样本前两行的最大绝对值，以进行归一化处理。

这种归一化方法可以确保数据在各个特征上具有相似的尺度，有助于提高模型的收敛速度和性能。

第四步

这段代码主要是将归一化后的数据按照一定的规则划分成训练集和测试集。具体步骤如下：

1.X_tr = data_norm(1:2,1:num_tr)';：从归一化后的数据中提取前 num_tr 个样本的前两列作为训练集的特征，存储在矩阵 X_tr 中。

2.Y_tr = data_norm(3,1:num_tr)';：从归一化后的数据中提取前 num_tr 个样本的第三列作为训练集的标签，存储在向量 Y_tr 中。

3.X_te = data_norm(1:2,1+num_tr:num_tr+num_te)';：从归一化后的数据中提取接下来的 num_te 个样本的前两列作为测试集的特征，存储在矩阵 X_te 中。

4.Y_te = data_norm(3,1+num_tr:num_tr+num_te)';：从归一化后的数据中提取接下来的 num_te 个样本的第三列作为测试集的标签，存储在向量 Y_te 中。

这种划分方式是常见的机器学习实验中常用的训练集和测试集划分方法，以便对模型进行训练和评估。这个过程是神经网络训练的基本步骤，通过这个过程，MLP 能够学习到从输入到输出的映射关系

第五步

1.初始化一个大小为 num_tr × num_tr 的矩阵 H，并将其填充为零。

2.使用两个嵌套的循环遍历训练数据集中的每对样本。

3.将计算得到的 RBF 核值与对应的标签 Y_tr(i) 和 Y_tr(j) 相乘，结果存储在矩阵 H 中。

最后，为了避免优化问题中的数值问题，向核矩阵 H 中添加了一小部分人为噪声。

第六步

1.初始化了一个大小为 1×num_tr1  的 α向量为零。

2.进行了 100 次迭代（由外层的 for 循环控制），在每次迭代中：

3.对训练数据集中的每一对样本，计算了一个新的核矩阵 K

4.根据计算得到的 K 和当前的 α 更新每个样本对应的 α。

5.更新过程中，对 α 做了一个截断操作，确保其在合理范围内。

6.在迭代结束后，对 α 做了一个常数减去操作（减去了 0.05）。

最后将得到的 α转置，并输出了它的大小以及核矩阵 H 的大小。

第七步

 1.svi = find(alpha > epsilon);：找到对应的 α 值大于某个阈值 ϵ 的索引，即找到支持向量的索引。

 2. nsv = length(svi);：计算支持向量的数量，即支持向量的个数。

 3.打印输出训练样本的数量 num_tr。

 4.计算模型的权重向量的范数的平方 w2，即 ∣w0∣2=αTHα，然后打印输出。

5.计算所有支持向量的 α 值之和 sum(alpha)，然后打印输出。

 6.打印输出支持向量的数量和占比。

参数调试和运行截图

距离参数（dist）:-3,  上界参数（C）:200,  RBF核函数的σ参数（p1）:1,  测试点数：24002.距离参数（dist）:-6,  上界参数（C）:200,  RBF核函数的σ参数（p1）:1,  测试点数：11003..距离参数（dist）:-5,  上界参数（C）:400,  RBF核函数的σ参数（p1）:1,  测试点数：4600 4.距离参数（dist）:-5,  上界参数（C）:200,  RBF核函数的σ参数（p1）:3,  测试点数：2200

5.距离参数（dist）:-1,  上界参数（C）200,  RBF核函数的σ参数（p1）:2,  测试点数：3400

实验结果与分析

1. 数据生成

生成了一个半月形数据集，其中包括300个训练数据点和1000个测试数据点。半月形数据的参数设置如下：

中心半径：10

半月形宽度：6

两半月之间的距离：-5

2. 数据预处理

为了提高模型的收敛性和准确性，对数据进行了去均值和归一化处理。具体步骤如下：

计算数据的均值，并将每个样本减去均值。

计算数据的绝对值最大值，并将每个样本除以该最大值。

3. 核矩阵构建

使用径向基函数（RBF）核函数来构建核矩阵。RBF核函数的参数设置为σ=1。通过遍历训练数据点，计算每对数据点之间的核函数值，构建核矩阵H。

4. 优化问题求解

在优化问题的求解过程中，我们使用了拉格朗日乘子法，通过迭代更新α值，逐步逼近最优解。在迭代过程中，设置了上界参数C=200，确保拉格朗日乘子α值在0和C之间。

5. 支持向量计算

通过迭代更新后，计算出支持向量的个数。支持向量的个数为53，占训练数据的17.7%。支持向量的存在对于SVM的决策边界有重要影响。

6. 模型训练结果

通过计算得到的拉格朗日乘子和核矩阵，可以得到模型的权重和偏置，并计算出决策边界。训练过程中，得到了以下结果：

样本点：300

参数C：200

|w0|^2：1.898

边界宽度：1.454

拉格朗日乘子和：848.973

7. 模型测试结果

在测试过程中，使用训练好的模型对1000个测试样本进行分类。分类结果显示，错误分类的样本点为78个，错误率为7.80%。这表明模型在测试数据上的表现良好，能够有效区分两类数据。

8. 图像展示

通过绘制训练和测试数据的分类结果图像，直观展示了SVM的分类效果。在图像中，可以清楚地看到决策边界、支持向量和各类数据点的分布情况。

实验总结

本实验通过使用支持向量机对半月形数据进行分类，验证了SVM在处理非线性可分数据时的有效性。通过构造RBF核矩阵，模型能够在高维特征空间中找到最优的分类超平面。数据预处理和归一化提高了模型的训练效果。实验中识别出的支持向量数量和测试错误率展示了模型的性能。

实验代码

% function err = NN_SVM(dist,C)%[nsv, alpha, b0, data] = NN_SVM(dist,C) %(X,Y,C)

%NN_SVM: Support Vector Machine for Classification

%

% Usage: [nsv alpha bias] = NN_SVM(X,Y,ker,C)

%

% Parameters: X - Training inputs

% Y - Training targets

% ker - kernel function

% C - upper bound (non-separable case)

% nsv - number of support vectors

% alpha - Lagrange Multipliers

% b0 - bias term

%

% Original code: Steve Gunn (srg@ecs.soton.ac.uk)

% Modified by Yanbo Xue (yxue@soma.mcmaster.ca) for non-commercial use

% March 28, 2007

%

%%--------------------- BEGIN OF MAIN ROUTINE -----------------------------

clc;clear;

dist=-5;

C=200;

%%============== Step 0A: Generate halfmoon data ==========================

rad = 10; % central radius of the half moon

width = 6; % width of the half moon

%dist = -4; % distance between two half moons

num_tr = 300; % # of training datasets

num_te = 1000; % # of testing datasets

num_samp = num_tr + num_te;% number of samples

fprintf('Support Vector Machine for Classification\n');

fprintf('_________________________________________\n');

fprintf('Generating halfmoon data ...\n');

fprintf(' ------------------------------------\n');

fprintf(' Points generated: %d\n',num_samp);

fprintf(' Halfmoon radius : %2.1f\n',rad);

fprintf(' Halfmoon width : %2.1f\n',width);

fprintf(' Distance : %2.1f\n',dist);

fprintf(' ------------------------------------\n');

seed=2e5;

rand('seed',seed);

[data, data_shuffled] = halfmoon(rad,width,dist,num_samp);

warning off;

%%================== Step 0B: Initialization of SVM =======================

%C = 25; % upper bound (non-separable case)

p1 = 1; % sigma for RBF, width

epsilon= 1e-6; % threshold for Support Vector Detection

err = 0; % a counter to denote the number of error outputs

b0 = 0; % Implicit bias, 0 for RBF kernal

%%==== Step 1: Preprocess the input data, remove mean and normalize =======

mean0 = [mean(data(1:2,:)')';0]; % mean of the original data

for i = 1:num_samp,

data_norm(:,i) = data_shuffled(:,i) - mean0;%每个样本与均值的偏差

end

abs(data_norm(1:2,:)');

max(abs(data_norm(1:2,:)'))

max0 = [max(abs(data_norm(1:2,:)'))';1];% max of the original data,求出max(|xi-E(x)|)

% max0

for i = 1:num_samp,

data_norm(:,i) = data_norm(:,i)./max0;

end

%%========== Step 2: Fetch Training and Testing Data ======================

X_tr = data_norm(1:2,1:num_tr)';

Y_tr = data_norm(3,1:num_tr)';

X_te = data_norm(1:2,1+num_tr:num_tr+num_te)';

Y_te = data_norm(3,1+num_tr:num_tr+num_te)';

%==================== Step 3: Construct the Kernel matrix =================

fprintf('Constructing RBF kernal matrix ...\n');

H = zeros(num_tr,num_tr);

for i=1:num_tr

for j=1:num_tr

H(i,j) = Y_tr(i)*Y_tr(j)*(exp(-(X_tr(i,:)-X_tr(j,:))*(X_tr(i,:)-X_tr(j,:))'/(2*p1^2)));

end

end

c = -ones(num_tr,1);

H = H + 1E-10*randn(size(H)); % add artifical noise to avoid Hessian problem.

%============ Step 4: Use QP to solve the Optimization Problem ============

% Set up the parameters for the Optimization problem

% fprintf('Optimizing using QP ...\n');

% vlb = zeros(num_tr,1); % Set the bounds: alphas >= 0

% vub = C*ones(num_tr,1);% alphas <= C

% x0 = zeros(num_tr,1); % The starting point is [0 0 0 0]

% neqcstr = 0; % Set the number of equality constraints (1 or 0)

% A = [];

% b = [];

% % Running QP to solve the optimization problem

% % The QP is implemented using C++ and then exported to a .dll file

% st = cputime;

% [alpha lambda how] = new_qp(H, c, A, b, vlb, vub, x0, neqcstr);

% de=size(alpha)

% ita=1;

alpha=zeros(1,num_tr);

for t=1:100

for i=1:num_tr

for j=1:num_tr

K1(i,j)=exp(-(X_tr(i,:)-X_tr(j,:))*(X_tr(i,:)-X_tr(j,:))'/(2*p1^2));

K(i,j)=alpha(j)*Y_tr(j)*K1(i,j);

end

alpha(i)=alpha(i)+(1-Y_tr(i)*sum(K(i,:)))/K1(i,i);

if alpha(i)<0

alpha(i)=0;

else

if alpha(i)>C

alpha(i)=C;

end

end

end

alpha(i)=alpha(i)-0.05;

end

alpha=alpha';

a1=size(alpha)

a2=size(H)

%======== Step 5: Compute the number of Support Vectors ===================

svi = find( alpha > epsilon);

nsv = length(svi);%支持向量的个数

nsv

fprintf(' ------------------------------------\n');

fprintf(' Points trained: %d\n',num_tr);

% fprintf(' Execution time: %4.1f seconds\n',cputime - st);

% fprintf(' Status : %s\n',how);

fprintf(' C : %f\n',C);

w2 = alpha'*H*alpha;%P176，公式6.30

fprintf(' |w0|^2 : %f\n',w2);

fprintf(' Margin : %f\n',2/sqrt(w2));

fprintf(' Sum alpha : %f\n',sum(alpha));

fprintf(' Support Vectors: %d (%3.1f%%)\n',nsv,100*nsv/num_tr);

fprintf(' ------------------------------------\n');

% %============= Step 6: Plot training result ===============================

% figure;

% svcplot_train(X_tr,Y_tr,alpha,b0,epsilon,p1,mean0,max0);

% xlabel('x');ylabel('y');

% title(['Classification using SVM with dist = ',num2str(dist),', radius = ',...

% num2str(rad), ' and width = ',num2str(width), '(Training)']);

%%=============== Step 7: Test the trained SVM ============================

fprintf('Testing the trained SVM ... \n');

H_te = zeros(num_te,num_tr);

for i=1:num_te

for j=1:num_tr

H_te(i,j) = Y_tr(j)*(exp(-(X_te(i,:)-X_tr(j,:))*(X_te(i,:)-X_tr(j,:))'/(2*p1^2)));%H_te

end

end

Y_pred = sign(H_te*alpha + b0);%177，公式6.34. %RBF，b0=0,对全部输入样本x进行分类

size(Y_pred)

% Calculate testing error rate

for i = 1:num_te,

if abs(Y_pred(i)-Y_te(i)) > 1E-6,%将第i个样本的经计算后的分类值与其实际值进行比较

% if abs(Y_pred(i)-Y_te(i)) > 1

err = err + 1;

end

end

fprintf(' ------------------------------------\n');

fprintf(' Points tested : %d\n',num_te);

fprintf(' Error points : %d (%5.2f%%)\n',err,(err/num_te)*100);

fprintf(' ------------------------------------\n');

% %%===================== Step 8: Plot testing result =======================

% figure;

% svcplot_test(X_tr,Y_tr,X_te,Y_pred,alpha,b0,epsilon,p1,mean0,max0);

% xlabel('x');ylabel('y');

% title(['Classification using SVM with d = ',num2str(dist), ', r = ',...

% num2str(rad), ' and width = ',num2str(width), '( C = ',num2str(C),' Testing)']);

%

% fprintf('Mission accomplished!\n');

% fprintf('_________________________________________\n');

%================== Final: Output data for further process ================

% data.X_tr = X_tr; % training input

% data.Y_tr = Y_tr; % training output

% data.X_te = X_te; % testing input

% data.Y_te = Y_te; % ideal testing output

% data.Y_pred=Y_pred;% predicted testing output

% data.mean = mean0;% mean value of the data (training + testing)

% data.max = max0; % max value of the data (training + testing)

% data.dist = dist; % parameter of halfmoon data: distance

% data.rad = rad; % parameter of halfmoon data: radius

% data.width= width;% parameter of halfmoon data: width

% data.err = err; % number of error outputs

%%-------------------- END OF MAIN ROUTINE --------------------------------

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% %%

%% Two subroutines to plot the training and testing result %%

%% %%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% function [h] = svcplot_train(X_tr,Y_tr,alpha,bias,epsilon,p1,mean0,max0)

%svcplot_train: Support Vector Machine Plotting Routine for Training data

%

% Usage: svcplot_train(X_tr,Y_tr,alpha,bias,epsilon,p1,mean0,max0)

%

% Parameters: X_tr - Training inputs

% Y_tr - Training outputs

% alpha - Lagrange Multipliers

% bias - Bias term

% epsilon - Threshold for support vector detection

% p1 - Sigma for RBF, spread width

% mean0 - Mean value of simulation data

% max0 - Max value of simulation data

hold on;

bias=0;

%%================== Plot Training Points =================================

for i = 1:size(Y_tr)

XX = X_tr(i,1)*max0(1) + mean0(1); % scale x back to original

YY = X_tr(i,2)*max0(2) + mean0(2); % scale y back to original

if (Y_tr(i) == 1)

h(1) = plot(XX,YY,'r+','LineWidth',1); % Class A

end

if (Y_tr(i) == -1)

h(2) = plot(XX,YY,'kx','LineWidth',1); % Class B

end

if ((abs(alpha(i)) > epsilon)&&Y_tr(i)==1)

plot(XX,YY,'ro','Linewidth',2); % Support Vector

end

if ((abs(alpha(i)) > epsilon)&&Y_tr(i)==-1)

plot(XX,YY,'ko','Linewidth',2); % Support Vector

end

end

%%================= Plot Decision Boundary ================================

xmin = min(X_tr(:,1));

xmax = max(X_tr(:,1));

ymin = min(X_tr(:,2));

ymax = max(X_tr(:,2));

[x,y]= meshgrid(xmin:(xmax-xmin)/50:xmax,ymin:(ymax-ymin)/50:ymax);

z = bias*ones(size(x));

%wh = waitbar(0,'Plotting...');

for x1 = 1 : size(x,1)

for y1 = 1 : size(x,2)

input(1) = x(x1,y1);

input(2) = y(x1,y1);

for i = 1 : length(Y_tr)

if (abs(alpha(i)) > epsilon)

z(x1,y1) = z(x1,y1) + Y_tr(i)*alpha(i)*(exp(-(input-X_tr(i,:))*(input-X_tr(i,:))'/(2*p1^2)));

end

end

end

%waitbar((x1)/size(x,1),wh)

%set(wh,'name',['Progress = ' sprintf('%2.1f',(x1)/size(x,1)*100) '%']);

end

%close(wh);

% Scale data points back to original value

x = x*max0(1) + mean0(1);

y = y*max0(2) + mean0(2);

xmin = xmin*max0(1) + mean0(1); % scale xmin

xmax = xmax*max0(1) + mean0(1); % scale xmax

ymin = ymin*max0(2) + mean0(2); % scale ymin

ymax = ymax*max0(2) + mean0(2); % scale ymax

set(gca,'XLim',[xmin xmax],'YLim',[ymin ymax]);

% Plot Boundary contour

contour(x,y,z,[0 0],'k','Linewidth',1);

contour(x,y,z,[-1 -1],'k:','Linewidth',1);

contour(x,y,z,[1 1],'k:','Linewidth',1);

hold off

axis on

%set(gca,'XTickLabelMode','auto');

%set(gca,'YTickLabelMode','auto');

% return

figure

% function svcplot_test(X_tr,Y_tr,X_te,Y_pred,alpha,bias,epsilon,p1,mean0,max0)

%svcplot_test: Support Vector Machine Plotting Routine for Testing data

%

%Usage: svcplot_test(X_tr,Y_tr,X_te,Y_pred,alpha,bias,epsilon,p1,mean0,max0)

%

% Parameters: X_tr - Training inputs

% Y_tr - Training outputs

% X_te - Testing inputs

% Y_pred - Testing outputs

% alpha - Lagrange Multipliers

% bias - Bias term

% epsilon - Threshold for support vector detection

% p1 - Sigma for RBF, spread width

% mean0 - Mean value of simulation data

% max0 - Max value of simulation data

%%================= Colormaping the figure here ===========================

%%=== In order to avoid the display problem of eps file in LaTeX. =========

hold on;

xmin = min(X_te(:,1));

xmax = max(X_te(:,1));

ymin = min(X_te(:,2));

ymax = max(X_te(:,2));

[x,y]= meshgrid(xmin:(xmax-xmin)/50:xmax,ymin:(ymax-ymin)/50:ymax);

z = bias*ones(size(x));

%wh = waitbar(0,'Plotting testing result...');

for x1 = 1 : size(x,1)

for y1 = 1 : size(x,2)

input(1) = x(x1,y1);, input(2) = y(x1,y1);

for i = 1 : length(Y_tr)

if (abs(alpha(i)) > epsilon)

z(x1,y1) = z(x1,y1) + Y_tr(i)*alpha(i)*(exp(-(input-X_tr(i,:))*(input-X_tr(i,:))'/(2*p1^2)));

end

end

end

%waitbar((x1)/size(x,1),wh)

%set(wh,'name',['Progress = ' sprintf('%2.1f',(x1)/size(x,1)*100) '%']);

end

%close(wh);

x = x*max0(1) + mean0(1); % scale x back to original

y = y*max0(2) + mean0(2); % scale y back to original

xmin = xmin*max0(1) + mean0(1); % scale xmin

xmax = xmax*max0(1) + mean0(1); % scale xmax

ymin = ymin*max0(2) + mean0(2); % scale ymin

ymax = ymax*max0(2) + mean0(2); % scale ymax

sp = pcolor(x,y,z);

%load red_black_colmap;

% colormap(red_black);

shading flat;

set(gca,'XLim',[xmin xmax],'YLim',[ymin ymax]);

num_te = size(X_te,1);

for i = 1:num_te,

XX = X_te(i,1)*max0(1) + mean0(1); % scale x back to original

YY = X_te(i,2)*max0(2) + mean0(2); % scale y back to original

if Y_pred(i) == 1,

plot(XX,YY,'r+');

else

plot(XX,YY,'kx');

end

end

%%======================= Plot decision boundary ==========================

contour(x,y,z,[0 0],'k','Linewidth',1);

contour(x,y,z,[-1 -1],'k:','Linewidth',1);

contour(x,y,z,[1 1],'k:','Linewidth',1);

% set(gca,'XTickLabelMode','auto');

% set(gca,'YTickLabelMode','auto');

hold off

axis on

% return