实验三 支持向量机

[代码地址]

(https://github.com/tangyuanzong/ML_SVM)

题目：

模式识别中著名的数据集。本实验通过花萼（sepal）和花瓣（petal）的长和宽，建立SVM分类器来判断样本属于山鸢尾（Iris Setosa）、变色鸢尾（Iris Versicolor）还是维吉尼亚鸢尾（Iris Virginica）。请按要求完成实验。

数据集：

文件iris.txt为该实验的数据集，包含150个样本，对应数据集的每行数据。每行数据包含每个样本的四个特征（按顺序分
鸢尾花数据集（Iris data set）是模别为花萼长度、花萼宽度、花瓣长度、花瓣宽度）和样本的类别信息（Iris Setosa、Iris Versicolor、Iris Virginica中的一种）。

实验要求：

（1）建立SVM分类器并用交叉验证法进行分析。
（2）利用PCA降维，将数据转化为二维，然后绘制出分类决策边界。

具体实现：

交叉验证方法的选取：

在机器学习里，通常来说我们不能将全部用于数据训练模型，否则我们将没有数据集对该模型进行验证，从而评估我们的模型的预测效果。为了解决这一问题，有许多验证方法。本次实验采用最简单的一种：

The Validation Set Approach

第一种是最简单的，也是很容易就想到的。我们可以把整个数据集分成两部分，一部分用于训练，一部分用于验证，这也就是我们经常提到的训练集（training set）和测试集（test set）。

例如，如上图所示，我们可以将蓝色部分的数据作为训练集（包含7、22、13等数据），将右侧的数据作为测试集（包含91等），这样通过在蓝色的训练集上训练模型，在测试集上观察不同模型不同参数对应的MSE的大小，就可以合适选择模型和参数了。但是该方法只用了部分数据进行模型的训练我们都知道，当用于模型训练的数据量越大时，训练出来的模型通常效果会越好。所以训练集和测试集的划分意味着我们无法充分利用我们手头已有的数据，所以得到的模型效果也会受到一定的影响。

我们可以将原始数据集的70%用作训练集，剩下 的30%用作测试集。

线性可分支持向量机：

假设给定一个特征空间上的训练数据集

学习的目标是在特征空间中找到一个分离超平面，能将实例分到不同的类。分离超平面对应于方程法向量方程它由法向量w和截距b决定，可用(w,b) 来表示。分离超平面将特征空间划分为两部分，一部分是正类，一部分是负类。法向量指向的一侧为正类，另一侧为负类。

一般地，当训练数据集线性可分时，存在无穷个分离超平面可将两类数据正确分开。感知机利用误分类最小的策略，求得分离超平面，不过这时的解有无穷多个。线性可分支持向量机利用间隔最大化求最优分离超平面，这时，解是唯一的。

给定线性可分训练数据集，通过间隔最大化或等价地求解相应的凸二次规划问题学习得到的分离超平面为：

以及相应的分类决策函数 ：

称为线性可分支持向量机。

考虑如图所示的二维特征空间中的分类问题。图中“。”表示正例，“x ”表示负例。训练数据集线性可分，这时有许多直线能将两类数据正确划分。线性可分支持向量机对应着将两类数据正确划分并且间隔最大的直线，如图所示。

函数间隔和几何间隔

对于给定的训练数据集T和超平面(w,b),定义超平 面(w,b)关于样本点(xi，yi）的函数间隔为

定义超平面(w,b)关于训练数据集T的函数间隔为超平面(w,b)关于T中所有样本点样本点.png的函数间隔之最小值，即

函数间隔可以表示分类预测的正确性及确信度。但是当我们成比例地改变w和b,例如将它们改为2w和2b,超平面并没有改变，但函数间隔却成为原来的2倍。这一事实启示我们，函数间隔不能用作间隔最大化的标准。

于是得到几何间隔

对于给定的训练数据集T和超平面(w，b),定义超平 面(w，b)关于样本点样本点的几何间隔为

所以同样可以得到数据集的几何间隔为：

我们发现函数间隔与集合间隔具有如下关系：

当我们成倍地增加w和b的值时，超平面并不会改变。

间隔最大化

支持向量机学习的基本想法是求解能够正确划分训练数据集并且几何间隔最大的分离超平面。对线性可分的训练数据集而言，线性可分分离超平面有无穷多个（等价于感知机)，但是几何间隔最大的分离超平面是唯一的。这里的间隔最大化又称为硬间隔最大化（与将要讨论的训练数据集近似线性可分时的软间隔最大化相对应)。

间隔最大化的直观解释是：对训练数据集找到几何间隔最大的超平面意味着以充分大的确信度对训练数据进行分类。也就是说，不仅将正负实例点分开，而且对最难分的实例点（离超平面最近的点）也有足够大的确信度将它们分开。这样的超平面应该对未知的新实例有很好的分类预测能力。

1、最大间隔分离超平面

下面考虑如何求得一个几何间隔最大的分离超平面，即最大间隔分离超平面。具体地，这个问题可以表示为下面的约束最优化问题：

考虑几何间隔和函数间隔的关系式几何间隔和函数间隔的关系式,可将这个问题改写为

我们上面说到：成倍地改变w和b的值，超平面不会改变。因此可以改写该问题，在约束条件的两边同时除以，使得最右边为1，并且最大化和最小化是等价的，因此得到新关系式：

这是一个凸二次规划（convex quadratic programming)问题。

支持向量和间隔边界

在线性可分情况下，训练数据集的样本点中与分离超平面距离最近的样本点的实例称为支持向量（support vector)。支持向量是使约束条件式成立的点。

如图所示，在H1和H2上的点就是支持向量。

注意到H1和H2平行，并且没有实例点落在它们中间。在H1和H2之间形成一条长带，分离超平面与它们平行且位于它们中央。长带的宽度，即H1和H2之间的距离称为间隔（margin)。间隔依赖于分离超平面的法向量w。H1和H2称为间隔边界。

在决定分离超平面时只有支持向量起作用，而其他实例点并不起作用。如果移动支持向量将改变所求的解；但是如果在间隔边界以外移动其他实例点，甚至去掉这些点，则解是不会改变的。由于支持向量在确定分离超平面中起着决定性作用，所以将这种分类模型称为支持向量机。支持向量的个数一般很少，所以支持向量机由很少的“重要的”训练样本确定。

学习的对偶算法

原始的最优化问题为：

首先构建拉格朗日函数(Lagrange function)。为此，对每一个不等式约束不等式约束。引进拉格朗日乘子（Lagrauge multiplier)拉格朗日乘子得到拉格朗日函数：

根据拉格朗日对偶性，原始问题的对偶问题是极大极小问题：

最后得到等价的最优化问题（证明过程省略）：

参数的计算：

通过引入拉格朗日函数，我们得到原问题：

的等价问题为：

w的求解：

b的求解：

分类超平面：

分类决策函数：

松弛变量

通常情况是，训练数据中有一些特异点（outlier),将这些特异点除去后，剩下大部分的样本点组成的集合是线性可分的。线性不可分意味着某些样本点不能满足函数间隔大于等于1的约束条件。
为了解决这个问题，可以对每个样本点引进一个松池变量，使函数间隔加上松弛变量大于等于1。这样，约束条件变为：

这样最优化问题变为：

SMO算法

输入：其中

输出：近似解

（1）取初始值

（2）选取优化变量 解析求解两个变量的最优化问题

求得最优解 更新为

（3）若在精度 范围内满足停机条件

其中，

则转(4）;否则令

（4） 取

python实现：

while (iter < maxIter):        alphaPairsChanged = 0   #alpha是否已经进行了优化        for i in range(m):            #   w = alpha * y * x;  f(x_i) = w^T * x_i + b            # 预测的类别            fXi = float(multiply(alphas,labelMat).T*dataMatrix*dataMatrix[i,:].T) + b                Ei = fXi - float(labelMat[i])   #得到误差，如果误差太大，检查是否可能被优化            #必须满足约束            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):                 j = selectJrand(i,m)                fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b                Ej = fXj - float(labelMat[j])                alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy()                                if (labelMat[i] != labelMat[j]):                                                              L = max(0, alphas[j] - alphas[i])                    H = min(C, C + alphas[j] - alphas[i])                else:                    L = max(0, alphas[j] + alphas[i] - C)                    H = min(C, alphas[j] + alphas[i])                if L==H:                    continue                # Eta = -(2 * K12 - K11 - K22)，且Eta非负，此处eta = -Eta则非正                eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T                if eta >= 0:                   continue                alphas[j] -= labelMat[j]*(Ei - Ej)/eta                alphas[j] = clipAlpha(alphas[j],H,L)                  #如果内层循环通过以上方法选择的α_2不能使目标函数有足够的下降，那么放弃α_1                if (abs(alphas[j] - alphaJold) < 0.00001):                     continue                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])                b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T                b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T                if (0 < alphas[i]) and (C > alphas[i]): b = b1                elif (0 < alphas[j]) and (C > alphas[j]): b = b2                else: b = (b1 + b2)/2.0                alphaPairsChanged += 1        if (alphaPairsChanged == 0): iter += 1        else: iter = 0

降维和PCA

常见的降维算法有主成分分析（principal component analysis,PCA）、因子分析（Factor Analysis）和独立成分分析（Independent Component Analysis，ICA），其中PCA是目前应用最为广泛的方法。

在PCA中，数据从原来的坐标系转换到新的坐标系，新坐标系的选择是由数据本身决定的。第一个坐标轴的选择是原始数据中方差最大的方向，从数据角度上来讲，这其实就是最重要的方向，即下图总直线B的方向。第二个坐标轴则是第一个的垂直或者说正交（orthogonal）方向，即下图中直线C的方向。该过程一直重复，重复的次数为原始数据中特征的数目。而这些方向所表示出的数据特征就被称为“主成分”。

那怎么来求出这些主成分呢？由线性代数的知识可以知道，通过数据集的协方差矩阵及其特征值分析，我们就可以求得这些主成分的值。一旦得到协方差矩阵的特征向量，就可以保留最大的N个值。然后可以通过把数据集乘上这N个特征向量转换到新的空间。

在python的numpy包中linalg模块的eig()方法可以用于求特征值和特征向量。
从上面的原理分析中我们可以得出讲数据转化成前N个主成分的伪代码如下：
去除平均值
计算协方差矩阵
计算协方差矩阵的特征值和特征向量
将特征值从大到小排序
保留最上面的N个特征向量
将数据转换到上述N个特征向量构建的新空间中

python实现：

def pca(dataMat, topNfeat):              #pca降维    re = dataMat - mean(dataMat, axis = 0)   #去除平均值     covMat = cov(re,rowvar=0) #计算协防差矩阵      eVals, eVects = linalg.eig(mat(covMat))      eValInd = argsort(eVals)      #从小到大对N个值排序      eValInd = eValInd[: -(topNfeat + 1) : -1]      reVects = eVects[:, eValInd]    lowDataMat = re * reVects   #转换到降维空间    return lowDataMat

多分类的解决方法：

在本次实验中，一共有3种类别：

Iris-setosaIris-versicolorIris-virginica

但是前面介绍的方法适合于求解2分类。

为此，我们可以想到的一种解决方法：通过求解多次来分类。

求解第一次： 将类别  Iris-setosa 标记为 +1，     其它类别标记为-1求解第二次： 将类别  Iris-versicolor 标记为+1，  其它类别标记为-1求解第三次： 将类别  Iris-virginica标记为+1 ，   其它类别标记为-1

我们首先来求解三次（使用全部数据）：

通过PCA得到实验结果：

我们发现并没有得到想要的线性分类的效果：

求解两次（类别 Iris-setosa和 Iris-virginica）：

通过PCA得到实验结果：

通过上面的结果，我们发现求解两次就可以了。

使用部分数据集求解：

（松弛变量=0.8 精度=0.0001 求解40次）

训练数据：

5.0,3.6,1.4,0.2,Iris-setosa5.4,3.9,1.7,0.4,Iris-setosa4.6,3.4,1.4,0.3,Iris-setosa5.0,3.4,1.5,0.2,Iris-setosa4.8,3.4,1.6,0.2,Iris-setosa4.8,3.0,1.4,0.1,Iris-setosa4.3,3.0,1.1,0.1,Iris-setosa5.8,4.0,1.2,0.2,Iris-setosa5.7,4.4,1.5,0.4,Iris-setosa5.4,3.9,1.3,0.4,Iris-setosa5.1,3.5,1.4,0.3,Iris-setosa5.1,3.7,1.5,0.4,Iris-setosa4.6,3.6,1.0,0.2,Iris-setosa5.1,3.3,1.7,0.5,Iris-setosa4.8,3.4,1.9,0.2,Iris-setosa5.0,3.0,1.6,0.2,Iris-setosa5.0,3.4,1.6,0.4,Iris-setosa5.2,3.5,1.5,0.2,Iris-setosa5.2,3.4,1.4,0.2,Iris-setosa4.7,3.2,1.6,0.2,Iris-setosa4.8,3.1,1.6,0.2,Iris-setosa5.4,3.4,1.5,0.4,Iris-setosa5.2,4.1,1.5,0.1,Iris-setosa5.5,4.2,1.4,0.2,Iris-setosa4.9,3.1,1.5,0.1,Iris-setosa4.4,3.0,1.3,0.2,Iris-setosa5.1,3.4,1.5,0.2,Iris-setosa5.0,3.5,1.3,0.3,Iris-setosa4.5,2.3,1.3,0.3,Iris-setosa4.4,3.2,1.3,0.2,Iris-setosa4.8,3.0,1.4,0.3,Iris-setosa5.1,3.8,1.6,0.2,Iris-setosa4.6,3.2,1.4,0.2,Iris-setosa5.3,3.7,1.5,0.2,Iris-setosa5.0,3.3,1.4,0.2,Iris-setosa7.0,3.2,4.7,1.4,Iris-versicolor6.4,3.2,4.5,1.5,Iris-versicolor6.9,3.1,4.9,1.5,Iris-versicolor5.5,2.3,4.0,1.3,Iris-versicolor6.5,2.8,4.6,1.5,Iris-versicolor5.7,2.8,4.5,1.3,Iris-versicolor6.3,3.3,4.7,1.6,Iris-versicolor4.9,2.4,3.3,1.0,Iris-versicolor5.2,2.7,3.9,1.4,Iris-versicolor5.0,2.0,3.5,1.0,Iris-versicolor5.9,3.0,4.2,1.5,Iris-versicolor6.0,2.2,4.0,1.0,Iris-versicolor6.1,2.9,4.7,1.4,Iris-versicolor5.9,3.2,4.8,1.8,Iris-versicolor6.4,2.9,4.3,1.3,Iris-versicolor6.6,3.0,4.4,1.4,Iris-versicolor6.8,2.8,4.8,1.4,Iris-versicolor6.7,3.0,5.0,1.7,Iris-versicolor6.0,2.9,4.5,1.5,Iris-versicolor5.7,2.6,3.5,1.0,Iris-versicolor5.5,2.4,3.8,1.1,Iris-versicolor5.5,2.4,3.7,1.0,Iris-versicolor5.8,2.7,3.9,1.2,Iris-versicolor6.0,2.7,5.1,1.6,Iris-versicolor5.4,3.0,4.5,1.5,Iris-versicolor6.0,3.4,4.5,1.6,Iris-versicolor6.1,3.0,4.6,1.4,Iris-versicolor5.8,2.6,4.0,1.2,Iris-versicolor5.0,2.3,3.3,1.0,Iris-versicolor5.6,2.7,4.2,1.3,Iris-versicolor5.7,3.0,4.2,1.2,Iris-versicolor5.7,2.9,4.2,1.3,Iris-versicolor6.2,2.9,4.3,1.3,Iris-versicolor5.1,2.5,3.0,1.1,Iris-versicolor5.7,2.8,4.1,1.3,Iris-versicolor6.3,3.3,6.0,2.5,Iris-virginica5.8,2.7,5.1,1.9,Iris-virginica7.1,3.0,5.9,2.1,Iris-virginica6.3,2.9,5.6,1.8,Iris-virginica6.5,3.0,5.8,2.2,Iris-virginica7.6,3.0,6.6,2.1,Iris-virginica4.9,2.5,4.5,1.7,Iris-virginica7.3,2.9,6.3,1.8,Iris-virginica6.7,2.5,5.8,1.8,Iris-virginica6.8,3.0,5.5,2.1,Iris-virginica5.7,2.5,5.0,2.0,Iris-virginica5.8,2.8,5.1,2.4,Iris-virginica6.4,3.2,5.3,2.3,Iris-virginica6.5,3.0,5.5,1.8,Iris-virginica7.7,3.8,6.7,2.2,Iris-virginica7.7,2.6,6.9,2.3,Iris-virginica6.0,2.2,5.0,1.5,Iris-virginica6.9,3.2,5.7,2.3,Iris-virginica5.6,2.8,4.9,2.0,Iris-virginica6.4,2.8,5.6,2.1,Iris-virginica7.2,3.0,5.8,1.6,Iris-virginica7.4,2.8,6.1,1.9,Iris-virginica6.1,2.6,5.6,1.4,Iris-virginica7.7,3.0,6.1,2.3,Iris-virginica6.9,3.1,5.4,2.1,Iris-virginica6.7,3.1,5.6,2.4,Iris-virginica6.9,3.1,5.1,2.3,Iris-virginica5.8,2.7,5.1,1.9,Iris-virginica6.8,3.2,5.9,2.3,Iris-virginica6.7,3.3,5.7,2.5,Iris-virginica6.7,3.0,5.2,2.3,Iris-virginica6.3,2.5,5.0,1.9,Iris-virginica6.5,3.0,5.2,2.0,Iris-virginica6.2,3.4,5.4,2.3,Iris-virginica5.9,3.0,5.1,1.8,Iris-virginica

测试数据：

4.7,3.2,1.3,0.2,Iris-setosa4.6,3.1,1.5,0.2,Iris-setosa5.7,3.8,1.7,0.3,Iris-setosa5.1,3.8,1.5,0.3,Iris-setosa5.4,3.4,1.7,0.2,Iris-setosa4.9,3.1,1.5,0.1,Iris-setosa5.0,3.2,1.2,0.2,Iris-setosa5.5,3.5,1.3,0.2,Iris-setosa4.4,2.9,1.4,0.2,Iris-setosa4.9,3.1,1.5,0.1,Iris-setosa5.4,3.7,1.5,0.2,Iris-setosa5.0,3.5,1.6,0.6,Iris-setosa5.1,3.8,1.9,0.4,Iris-setosa5.1,3.5,1.4,0.2,Iris-setosa4.9,3.0,1.4,0.2,Iris-setosa5.6,2.9,3.6,1.3,Iris-versicolor6.7,3.1,4.4,1.4,Iris-versicolor5.6,3.0,4.5,1.5,Iris-versicolor5.8,2.7,4.1,1.0,Iris-versicolor6.2,2.2,4.5,1.5,Iris-versicolor5.6,2.5,3.9,1.1,Iris-versicolor6.3,2.3,4.4,1.3,Iris-versicolor5.6,3.0,4.1,1.3,Iris-versicolor5.5,2.5,4.0,1.3,Iris-versicolor6.1,2.8,4.0,1.3,Iris-versicolor6.3,2.5,4.9,1.5,Iris-versicolor6.1,2.8,4.7,1.2,Iris-versicolor6.7,3.1,4.7,1.5,Iris-versicolor5.5,2.6,4.4,1.2,Iris-versicolor6.6,2.9,4.6,1.3,Iris-versicolor6.3,3.4,5.6,2.4,Iris-virginica6.4,3.1,5.5,1.8,Iris-virginica6.0,3.0,4.8,1.8,Iris-virginica7.7,2.8,6.7,2.0,Iris-virginica6.3,2.7,4.9,1.8,Iris-virginica6.7,3.3,5.7,2.1,Iris-virginica7.2,3.2,6.0,1.8,Iris-virginica6.2,2.8,4.8,1.8,Iris-virginica6.1,3.0,4.9,1.8,Iris-virginica7.2,3.6,6.1,2.5,Iris-virginica6.5,3.2,5.1,2.0,Iris-virginica6.4,2.7,5.3,1.9,Iris-virginica7.9,3.8,6.4,2.0,Iris-virginica6.4,2.8,5.6,2.2,Iris-virginica6.3,2.8,5.1,1.5,Iris-virginica

将测试数据降纬，然后显示：

我们发现：

分类器工作正常，正确率97.7777777778