机器学习实战笔记(Python实现)-05-支持向量机(SVM)
---------------------------------------------------------------------------------------
本系列文章为《机器学习实战》学习笔记，内容整理自书本，网络以及自己的理解，如有错误欢迎指正。
源码在Python3.5上测试均通过，代码及数据 --> https://github.com/Wellat/MLaction
---------------------------------------------------------------------------------------
1、支持向量机概述
1.1 原理简述
所谓支持向量机，顾名思义，分为两个部分了解，一什么是支持向量(简单来说，就是支持（或支撑）平面上把两类类别划分开来的超平面的向量点)，二这里的“机（machine,机器）”是一个算法。在机器学习领域，常把一些算法看做是一个机器，如分类机(当然，也叫做分类器)，而支持向量机本身便是一种监督式学习的方法，它广泛的应用于统计分类以及回归分析中。
用一个二维空间里仅有两类样本的分类问题来举个小例子。 如图所示
正方形和圆圈是要区分的两个类别，在二维平面中它们的样本如上图所示。中间的直线就是一个分类函数，它可以将两类样本完全分开。在这种情况下边缘实心的几个数据点就叫做支持向量，这也是支持向量机这个分类算法名字的来源。在SVM中，我们寻找一条最优的分界线使得它到两边的边界都最大（最大化支持向量到分隔线(面)的距离）。
一般的，如果一个线性函数能够将样本完全正确的分开，就称这些数据是线性可分的，否则称为非线性可分的。什么叫线性函数呢？在一维空间里就是一个点，在二维空间里就是一条直线，三维空间里就是一个平面，可以如此想象下去，如果不关注空间的维数，这种线性函数还有一个统一的名称——超平面。
在实际中，我们经常会遇到线性不可分的样例，此时，我们的常用做法是把样例特征通过某种核函数映射到高维空间中去，使其线性可分。下图即是映射前后的结果，将坐标轴经过适当的旋转，就可以很明显地看出，数据是可以通过一个平面来分开的。
核函数的价值在于它虽然也是将特征进行从低维到高维的转换，但核函数厉害在它事先在低维上进行计算，而将实质上的分类效果表现在了高维上，避免了直接在高维空间中的复杂计算。
1.2 特点
将SVM算法和前面介绍的Logistic回归和决策树算法作对比，如下图所示，我们能直观看到SVM作为非线性分类器的优势。
优点：泛化错误率低，计算开销不大，结果易解释
缺点：对参数调节和核函数的选择敏感，原始分类器不加修改仅适用于处理二类问题
适用数据类型：数值型和标称型数据，类别标签为+1和-1
2、寻找最大分类间隔
2.1 关于线性分类器
上图中红蓝两类数据点可以用线性函数 g(x)=w*x+b 区分开，关于这个线性函数要注意一下三点
- 式中的 x 不是二维坐标系中的横轴，而是样本的向量表示，例如一个样本点的坐标是(3,8)，则x=(3,8)，而不是x=3
- 这个形式并不局限于二维的情况，在 n 维空间中仍然可以使用这个表达式，只是式中的 w 成为了 n 维向量
- g(x) 不是中间那条直线的表达式，中间那条直线的表达式是g(x)=0，即 w*x+b = 0，也把这个函数叫做分类面
易知中间那条分界线并不是唯一的，把它稍微旋转或平移一下，仍然可以达到上面说的效果。那就牵涉到一个问题，对同一个问题，哪一个函数更好？通常的衡量指标叫做分类间隔。
2.2 分类间隔
在监督学习中，每一个样本由一个特征向量和一个类别标签组成，如下：
Di=(xi,yi)
xi 就是特征向量，就是yi 分类标记。
在二元的线性分类中， 这个表示分类的标记只有两个值，+1和-1。有了这种表示法，我们就可以定义一个样本点到某个超平面的间隔（函数间隔）：
注意到如果某个样本属于该类别的话，那么wi*x+b > 0（这是因为我们所选的g(x)=wx+b就是通过大于0还是小于0来判断分类），而yi也大于0；若不属于该类别的话，那么wi*x+b < 0 ，而yi也小于0，这意味着yi(w*xi+b)总是大于0，而它的值就等于|wxi+b|, 也即|g(xi)|.
现在把w和b进行归一化，即用w/||w||和b/||w||分别代替原来的w和b，那么间隔就可以写成如下形式，叫做几何间隔，几何间隔所表示的正是点到超平面的欧式距离。
其中，||w|| 叫做向量w的范数，范数是对向量长度的一种度量。我们常说的向量长度其实指的是它的2-范数，范数最一般的表示形式为p-范数，可以写成如下形式：
好，到这现在的目标就是找出分离器定义中的w和b。而前面我们知道SVM依据最大化支持向量到分隔线(面)的距离，为此我们必须找到具有最小间隔的数据点，而这些数据点也就是前面提到的支持向量。一旦找到具有最小间隔的数据点，我们就需要对该间隔最大化。这就可以写作：
直接求解上述问题相当困难，可以将它转换成为另一种更容易求解的形式。考察上式中大括号内的部分，我们可以固定其中一个因子而最大化其他因子。如果令所有支持向量 
的都为1，那么就可以通过求 ||w|| 的最小值来得到最终解。但是， 并非所有数据点的 
都等于1，只有那些离分隔超平面最近的点得到的值才为1。而离超平面越远的数据点，其值也就越大。
上述问题是一个带约束条件的优化问题，这里的约束条件是
大于等于1，对于这类问题可以引入拉格朗日乘子法求解。由于这里的约束条件都是基于数据点的，因此我们就可以将超平面写成数据点的形式。于是，优化目标函数最后可以写成：
其约束条件为：
这里我们有个假设即数据100%线性可分，但是实际上并不是所有数据都能达到该要求，因此我们可以通过引入所谓松弛变量——C，来允许有些数据点可以处于分隔面的错误一侧，这时约束条件变为：
这里的常数C用于控制 “最大化间隔” 和 “保证大部分点的函数间隔小于1.0” 这两个目标的权重。在优化算法的实现代码中，常数C是一个参数，因此我们就可以通过调节该参数得到不同的结果。一旦求出了所有的alpha，那么分隔超平面就可以通过这些alpha来表达。
3、SMO算法
本节对前面的两个式子进行优化，一个是最小化的目标函数，一个是遵循的约束条件。优化的方法有很多种，但本章只关注其中最流行的一种实现，即序列最小优化(SMO)算法。SMO算法的目标是求出一系列alpha和b，一旦求出了这些alpha，就很容易计算出权重向量w并得到分隔超平面（2维平面中就是直线），再能够将之用于数据分类。
SMO算法的工作原理是：每次循环中选择两个alpha进行优化处理。一旦找到一对合适的alpha,那么就增大其中一个同时减小另一个。这里所谓的“合适” 就是指两个alpha必须要符合一定的条件，即这两个alpha必须要在间隔边界之外，且这两个alpha还没有进行过区间化处理或者不在边界上。
3.1 简化版SMO算法处理小规模数据集
SMO算法的完整实现需要大量代码。在接下来的第一个例子中，我们将会对算法进行简化处理，以便了解算法的基本工作思路，之后再基于简化版给出完整版。
该算法伪代码大致如下：
创建一个alpha向量并将其初始化为O向量
当迭代次数小于最大迭代次数时（外循环）
对数据集中的每个数据向量（内循环）：
如果该数据向量可以被优化：
随机选择另外一个数据向量
同时优化这两个向量
如果两个向量都不能被优化，退出内循环
如果所有向量都没被优化，增加迭代数目，继续下一次循环
简化版SMO算法：
1 def loadDataSet(fileName): 2 dataMat = []; labelMat = [] 3 fr = open(fileName) 4 for line in fr.readlines(): 5 lineArr = line.strip().split('\t') 6 dataMat.append([float(lineArr[0]), float(lineArr[1])]) 7 labelMat.append(float(lineArr[2])) 8 return dataMat,labelMat 9 10 def selectJrand(i,m): 11 ''' 12 在0-m中随机返回一个不等于i的数 13 ''' 14 j=i 15 while (j==i): 16 j = int(random.uniform(0,m)) 17 return j 18 19 def clipAlpha(aj,H,L): 20 #小于L或大于H的aj将被调整为L或H 21 if aj > H: 22 aj = H 23 if L > aj: 24 aj = L 25 return aj 26 27 def smoSimple(dataMatIn, classLabels, C, toler, maxIter): 28 ''' 29 5个输入参数分别为：数据集、类别标签、常数C、容错率、循环次数 30 ''' 31 dataMatrix = mat(dataMatIn); labelMat = mat(classLabels).transpose() 32 b = 0; m,n = shape(dataMatrix) 33 alphas = mat(zeros((m,1))) 34 iter = 0 35 while (iter < maxIter): 36 #用于记录alpha是否已经优化 37 alphaPairsChanged = 0 38 for i in range(m): 39 #fXi为预测类别 40 fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b 41 #Ei预测值与实际的误差 42 Ei = fXi - float(labelMat[i]) 43 #如果Ei大于容错，且0<alpha<C 进入优化 44 #由于后面alpha小于0或大于C时将被调整为0或C，所以一旦在该if语句中它们等于这两个值的话， 45 #那么它们就已经在“边界”上了，因而不再能够减小或增大，因此也就不值得再对它们进行优化了。 46 if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)): 47 #随机选择第二个alpha 48 j = selectJrand(i,m) 49 fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b 50 Ej = fXj - float(labelMat[j]) 51 alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy(); 52 #计算L和H。保证alpha在0和C之间 53 if (labelMat[i] != labelMat[j]): 54 L = max(0, alphas[j] - alphas[i]) 55 H = min(C, C + alphas[j] - alphas[i]) 56 else: 57 L = max(0, alphas[j] + alphas[i] - C) 58 H = min(C, alphas[j] + alphas[i]) 59 #如果L==H，本次循环结束 60 if L==H: print("L==H"); continue 61 #eta是alpha[j]的最优修改量 62 eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T 63 if eta >= 0: print("eta>=0"); continue 64 alphas[j] -= labelMat[j]*(Ei - Ej)/eta 65 alphas[j] = clipAlpha(alphas[j],H,L) 66 #alpha[j]轻微改变，退出for循环 67 if (abs(alphas[j] - alphaJold) < 0.00001): print("j not moving enough"); continue 68 #对alphas[i]进行修改，修改量与alpha[j]相同，但是方向相反 69 alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j]) 70 #对alpha[i]和alpha[j]进行优化之后，设置一个常数项b 71 b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T 72 b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T 73 if (0 < alphas[i]) and (C > alphas[i]): b = b1 74 elif (0 < alphas[j]) and (C > alphas[j]): b = b2 75 else: b = (b1 + b2)/2.0 76 alphaPairsChanged += 1 77 print("iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)) 78 if (alphaPairsChanged == 0): iter += 1 79 else: iter = 0 80 print("iteration number: %d" % iter) 81 return b,alphas
测试：
1 if __name__ == "__main__": 2 dataMat,labelMat = loadDataSet('testSet.txt') 3 b,alphas = smoSimple(dataMat,labelMat,0.6,0.001,40) 4 #支持向量的位置 5 for i in range(100): 6 if alphas[i]>0.0: print(dataMat[i],labelMat[i])
3.2 完整Platt SMO算法加速优化
完整版SMO算法通过一个外循环和内循环分别选择两个alpha值，算法中许多计算过程要反复调用，这里将它们单独提取出方便使用，具体如下：
1 def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)): #full Platt SMO 2 ''' 3 主函数--外循环 4 ''' 5 #用构建的数据结构容纳所有数据 6 oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler) 7 iter = 0 8 entireSet = True; alphaPairsChanged = 0 9 #第一个alpha值的选择会在两种方式之间进行交替 10 #一种方式是在所有数据集上进行单遍扫描，另一种方式则是在非边界alpha中实现单遍扫描 11 #这里非边界指的是那些不等于边界0或C的alpha值 12 #同时，这里会跳过那些已知的不会改变的alpha值 13 while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)): 14 alphaPairsChanged = 0 15 if entireSet: 16 #在数据集上遍历任意可能的alpha 17 for i in range(oS.m): 18 alphaPairsChanged += innerL(i,oS) 19 print("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)) 20 iter += 1 21 else:#遍历非边界alpha值 22 nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0] 23 for i in nonBoundIs: 24 alphaPairsChanged += innerL(i,oS) 25 print("non-bound, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)) 26 iter += 1 27 if entireSet: entireSet = False #控制循环的flag 28 elif (alphaPairsChanged == 0): entireSet = True 29 print("iteration number: %d" % iter) 30 return oS.b,oS.alphas
其调用到的子函数说明如下：
1 class optStruct: 2 '''数据结构保存数据''' 3 def __init__(self,dataMatIn, classLabels, C, toler): #初始化参数 4 self.X = dataMatIn 5 self.labelMat = classLabels 6 self.C = C 7 self.tol = toler 8 self.m = shape(dataMatIn)[0] 9 self.alphas = mat(zeros((self.m,1))) 10 self.b = 0 11 self.eCache = mat(zeros((self.m,2))) #第一列为有效标志flag 12 13 def innerL(i, oS): 14 ''' 15 内循环函数 16 此处代码集合和smoSimple()函数一模一样，主要变化有 17 1.数据用参数oS传递 18 2.用selectJ函数替代selectJrand来选择第二个alpha的值 19 3.在alpha值改变时更新Ecache 20 21 ''' 22 Ei = calcEk(oS, i) 23 if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)): 24 j,Ej = selectJ(i, oS, Ei) #此处和简化版不同 25 alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy(); 26 if (oS.labelMat[i] != oS.labelMat[j]): 27 L = max(0, oS.alphas[j] - oS.alphas[i]) 28 H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i]) 29 else: 30 L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C) 31 H = min(oS.C, oS.alphas[j] + oS.alphas[i]) 32 if L==H: print("L==H"); return 0 33 eta = 2.0 * oS.X[i,:]*oS.X[j,:].T - oS.X[i,:]*oS.X[i,:].T - oS.X[j,:]*oS.X[j,:].T 34 # eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j] #changed for kernel 35 if eta >= 0: print("eta>=0"); return 0 36 oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta 37 oS.alphas[j] = clipAlpha(oS.alphas[j],H,L) 38 updateEk(oS, j) #更新误差缓存 39 if (abs(oS.alphas[j] - alphaJold) < 0.00001): print("j not moving enough"); return 0 40 oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j]) 41 updateEk(oS, i) #更新误差缓存 42 b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[i,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[i,:]*oS.X[j,:].T43 b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[j,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[j,:]*oS.X[j,:].T44 if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1 45 elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2 46 else: oS.b = (b1 + b2)/2.0 47 #alpha值发生变动返回1，否则返回0 48 return 1 49 else: return 0 50 51 def selectJ(i, oS, Ei): 52 ''' 53 选择第二个alpha，即内循环的alpha值 54 依据最大步长(max(abs(Ei - Ej)))选择 55 ''' 56 maxK = -1; maxDeltaE = 0; Ej = 0 57 #存储Ei，第一位为有效标记 58 oS.eCache[i] = [1,Ei] 59 #构建一个非零表，返回非零E值所对应的index 60 validEcacheList = nonzero(oS.eCache[:,0].A)[0] 61 if (len(validEcacheList)) > 1: 62 for k in validEcacheList: 63 if k == i: continue 64 Ek = calcEk(oS, k) 65 deltaE = abs(Ei - Ek) 66 if (deltaE > maxDeltaE): 67 maxK = k; maxDeltaE = deltaE; Ej = Ek 68 return maxK, Ej 69 else:#第一次循环随机选一个alpha值 70 j = selectJrand(i, oS.m) 71 Ej = calcEk(oS, j) 72 return j, Ej 73 74 def calcEk(oS, k): 75 '''计算预测误差''' 76 fXk = float(multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T)) + oS.b
77 Ek = fXk - float(oS.labelMat[k]) 78 return Ek 79 80 def updateEk(oS, k): 81 '''计算误差值，并存入缓存中''' 82 Ek = calcEk(oS, k) 83 oS.eCache[k] = [1,Ek]
3.3 分类
花了大量时间来计算alpha值和b，现在利用它们来完成分类任务。依据y=w*x+b。
以下函数计算w
1 def calcWs(alphas,dataArr,classLabels): 2 ''' 3 实际计算出的alpha值为0，而非零alpha所对应的也就是支持向量 4 本计算函数遍历所有alpha，但最终起作用的只有支持向量 5 ''' 6 X = mat(dataArr); labelMat = mat(classLabels).transpose() 7 m,n = shape(X) 8 w = zeros((n,1)) 9 for i in range(m): 10 w += multiply(alphas[i]*labelMat[i],X[i,:].T) 11 return w
现在可以采用以下方式测试分类：
以上得出的值大于0，那么其属于1类；如果该值小于等于0，那么则属于-1类。则上例被分为-1类，可以通过如下命令确认分类结果的正确性：
对所有数据的分类结果：
4、在复杂数据上应用核函数
上节两类数据点可以通过一条直线划分开，对于两类分布在一个圆的内部和外部的数据点，如下图，我们就要使用核函数将数据转换成已与分类器理解的形式了。本节介绍一种较为流行的核函数——径向基核函数（RBF）。
4.1 径向基核函数
径向基函数是一个采用向量作为自变量的函数，能够基于向量距离运算输岀一个标量。这个距离可以是从<0,0>向量或者其他向量开始计算的距离。接下来，我们将会使用到径向基函数的高斯版本，其具体公式为：
其中，
是用户定义的用于确定到达率（reach) 或者说函数值跌落到0的速度参数。
用Python实现，添加代码如下：
1 def kernelTrans(X, A, kTup): 2 ''' 3 核函数 4 X(m,n):支持向量集 5 A(1,n):待变换的向量 6 kTup:含两个参数--①所用核函数的类型 ②速度参数sigma 7 输出K(m,1) 8 ''' 9 m,n = shape(X) 10 K = mat(zeros((m,1))) 11 if kTup[0]=='lin': K = X * A.T #线性核 12 elif kTup[0]=='rbf': 13 for j in range(m): 14 deltaRow = X[j,:] - A 15 K[j] = deltaRow*deltaRow.T 16 K = exp(K/(-1*kTup[1]**2)) # 17 else: raise NameError('Error：That Kernel is not recognized！') 18 return K
为顺利使用核函数，还需对上节代码稍作调整，具体变动如下：
optStruct类，引入新变量kTup：
1 class optStruct: 2 '''数据结构保存数据''' 3 def __init__(self,dataMatIn, classLabels, C, toler, kTup): #初始化参数 4 self.X = dataMatIn 5 self.labelMat = classLabels 6 self.C = C 7 self.tol = toler 8 self.m = shape(dataMatIn)[0] 9 self.alphas = mat(zeros((self.m,1))) 10 self.b = 0 11 self.eCache = mat(zeros((self.m,2))) #第一列为有效标志flag 12 self.K = mat(zeros((self.m,self.m))) #for kernel 13 for i in range(self.m): 14 self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)
innerL()函数修改如下：
1 ··· 2 eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j] #changed for kernel 3 ··· 4 b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,i] - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[i,j] 5 b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,j]- oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[j,j] 6 ···
calcEk()函数：
1 def calcEk(oS, k): 2 '''计算预测误差''' 3 fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)#changed for kernel 4 Ek = fXk - float(oS.labelMat[k]) 5 return Ek
4.2 在测试中使用核函数
前面提到径向基函数有一个用户定义的输入sigma——，我们需要确定它的大小，再利用该函数构建出一个分类器。测试代码如下：
1 def testRbf(k1=1.3): 2 ''' 3 利用核函数进行分类的径向基测试函数 4 ''' 5 dataArr,labelArr = loadDataSet('testSetRBF.txt') 6 #训练 7 b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) #C=200 important 8 datMat=mat(dataArr); labelMat = mat(labelArr).transpose() 9 svInd=nonzero(alphas.A>0)[0] 10 sVs=datMat[svInd] #获取支持向量 11 labelSV = labelMat[svInd]; 12 print("there are %d Support Vectors" % shape(sVs)[0]) 13 m,n = shape(datMat) 14 #分类 15 errorCount = 0 16 for i in range(m): 17 #数据转换 18 kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1)) 19 #将转换后的数据与前面的alpha及类标签值求积得到预测值 20 predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b 21 if sign(predict)!=sign(labelArr[i]): errorCount += 1 22 print("the training error rate is: %f" % (float(errorCount)/m)) 23 #测试 24 dataArr,labelArr = loadDataSet('testSetRBF2.txt') 25 errorCount = 0 26 datMat=mat(dataArr); labelMat = mat(labelArr).transpose() 27 m,n = shape(datMat) 28 for i in range(m): 29 kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1)) 30 predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b 31 if sign(predict)!=sign(labelArr[i]): errorCount += 1 32 print("the test error rate is: %f" % (float(errorCount)/m))
可以尝试使用不同的k1参数来观察不同的分类情况，事实表明k1会影响支持向量的数量，如果支持向量太少，就可能会得到一个很差的决策边界，如果支持向量太多，也就相当于每次都利用整个数据集进行分类，这种分类方法就变成了k近邻。
5、实例：基于SVM的数字识别
使用kNN方法进行手写体识别可以得到不错的效果，这在前面有过详细介绍，详情可以点击这里查看。使用kNN需要保留所有的训练样本，而如果用支持向量机方法则只需保留支持向量，且能获得可比的效果。本节介绍基于SVM的手写数字识别。
本质上，支持向量机是一个二类分类器，其分类结果不是+1就是-1。所以在本例中，只识别数字9，即一旦碰到9则输出类别标签-1，否则输出+1。
1 def testDigits(kTup=('rbf', 10)): 2 ''' 3 基于SVM的手写数字识别 4 ''' 5 dataArr,labelArr = loadImages('digits/trainingDigits') 6 b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, kTup) 7 datMat=mat(dataArr); labelMat = mat(labelArr).transpose() 8 svInd=nonzero(alphas.A>0)[0] 9 sVs=datMat[svInd] 10 labelSV = labelMat[svInd]; 11 print("there are %d Support Vectors" % shape(sVs)[0]) 12 m,n = shape(datMat) 13 errorCount = 0 14 for i in range(m): 15 kernelEval = kernelTrans(sVs,datMat[i,:],kTup) 16 predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b 17 if sign(predict)!=sign(labelArr[i]): errorCount += 1 18 print("the training error rate is: %f" % (float(errorCount)/m)) 19 dataArr,labelArr = loadImages('digits/testDigits') 20 errorCount = 0 21 datMat=mat(dataArr); labelMat = mat(labelArr).transpose() 22 m,n = shape(datMat) 23 for i in range(m): 24 kernelEval = kernelTrans(sVs,datMat[i,:],kTup) 25 predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b 26 if sign(predict)!=sign(labelArr[i]): errorCount += 1 27 print("the test error rate is: %f" % (float(errorCount)/m)) 28 29 def loadImages(dirName): 30 ''' 31 导入数据 32 ''' 33 from os import listdir 34 hwLabels = [] 35 trainingFileList = listdir(dirName) 36 m = len(trainingFileList) 37 trainingMat = zeros((m,1024)) 38 for i in range(m): 39 #从文件名中解析出当前图像的标签，也就是数字是几 40 fileNameStr = trainingFileList[i] 41 fileStr = fileNameStr.split('.')[0] 42 classNumStr = int(fileStr.split('_')[0]) 43 #是数字9类标签设为-1，否则为+1 44 if classNumStr == 9: hwLabels.append(-1) 45 else: hwLabels.append(1) 46 #将图片数据转换为向量 47 trainingMat[i,:] = img2vector('%s/%s' % (dirName, fileNameStr)) 48 return trainingMat, hwLabels 49 50 def img2vector(filename): 51 """ 52 将图片数据矩阵转换为向量 53 每张图片是32*32像素，也就是一共1024个字节。 54 因此转换的时候，每行表示一个样本，每个样本含1024个字节。 55 """ 56 #每个样本数据是1024=32*32个字节 57 returnVect = zeros((1,1024)) 58 fr = open(filename) 59 #循环读取32行，32列。 60 for i in range(32): 61 lineStr = fr.readline() 62 for j in range(32): 63 returnVect[0,32*i+j] = int(lineStr[j]) 64 return returnVect
testDigits()与上节testRbf()的代码几乎一样，唯一区别就是它调用了loadImages()函数来获得类别标签和数据。该函数与img2vector()函数在介绍kNN算法时均有应用到。
改变sigma的值，并尝试线性核函数，总结得到如下结果：
参考文献：
《机器学习实战》 Peter Harrington
http://www.matlabsky.com/thread-10317-1-1.html
http://blog.csdn.net/v_july_v/article/details/7624837
THE END.
