【机器学习】支持向量机的底层实现-Python

xiaoxiao2021-02-28 49

支持向量机数学原理

1. 支持向量机原理（一）线性支持向量机 2. 支持向量机原理（二）线性支持向量机的软间隔最大化模型 3. 支持向量机原理（三）线性不可分支持向量机与核函数 4. 支持向量机原理（四）SMO算法原理

应用简化版SMO算法处理小规模数据集　　

　　我们首先从一个已经准备好的文件”testSet.txt”中读入数据，然后运用简化版的SMO算法进行实验，简化版SMO的具体流程请看代码注释。文件”testSet.txt”的数据集如下图所示：　　

''' def loadDataSet(fileName):从数据集文件fileName中读入数据，返回数据列表dataMat和标签列表labelMat ''' def loadDataSet(fileName): dataMat = []; labelMat = [] fr = open(fileName) for line in fr.readlines(): lineArr = line.strip().split('\t') dataMat.append([float(lineArr[0]), float(lineArr[1])]) labelMat.append(float(lineArr[2])) return dataMat,labelMat ''' def selectJrand(i,m):从[1,m)中，随机返回一个不等于i的数 ''' def selectJrand(i,m): j=i #we want to select any J not equal to i while (j==i): j = int(random.uniform(0,m)) return j ''' def clipAlpha(aj,H,L): 这个是剪辑函数，如果aj大于边界H，则返回H；如果aj小于边界L，则返回L；如果L<=aj<=H，则直接返回aj ''' def clipAlpha(aj,H,L): if aj > H: aj = H if L > aj: aj = L return aj ''' def smoSimple(dataMatIn, classLabels, C, toler, maxIter):简化版的SMO算法。 dataMatIn是数据列表，classLabels是标签列表，toler是松弛因子，C是toler的惩罚系数，maxIter是最大迭代次数。该函数首先遍历每个ai，然后再寻找当前ai对应合适的aj。如果遍历所有的ai都没有找到合适的aj，则重新遍历ai，但是该过程的最大次数为maxIter。该函数返回列表alphas、最佳的b ''' def smoSimple(dataMatIn, classLabels, C, toler, maxIter): dataMatrix = mat(dataMatIn); labelMat = mat(classLabels).transpose() #将classLabels从行向量转化为列向量 b = 0; m,n = shape(dataMatrix) alphas = mat(zeros((m,1))) iter = 0 #iter表示寻找ai-aj对的连续失败次数 while (iter < maxIter): alphaPairsChanged = 0 for i in range(m): #这是外循环，用于遍历所有的参数ai fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b # multiply(alphas,labelMat).T表示列向量alphas与列向量labelMat对应元素相乘，然后再转置为行向量 Ei = fXi - float(labelMat[i]) if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)): # ‘if’检查是否违反KKT条件,如果没有违反则继续选择第二个参数aj j = selectJrand(i,m) #随机从区间[i,m)中选出一个aj 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: print ("L==H"); continue #重新选择ai eta = dataMatrix[i,:]*dataMatrix[i,:].T + dataMatrix[j,:]*dataMatrix[j,:].T - 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T if eta <= 0: print ("eta<=0"); continue #eta<=0不能符号条件，进行下一轮,重新选择ai alphas[j] += labelMat[j]*(Ei - Ej)/eta #计算第二个参数aj alphas[j] = clipAlpha(alphas[j],H,L) if (abs(alphas[j] - alphaJold) < 0.00001): print ("j not moving enough"); continue alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j]) #通过aj和ai的关系，由aj求出ai 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 #ai与aj已经求出，alphaPairsChanged += 1 print ("iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)) if (alphaPairsChanged == 0): iter += 1 #如果遍历所有的ai都没有寻找到合适的aj，则表示寻找ai-aj对的连续失败次数iter+=1 else: iter = 0 #找到了合适的ai-aj对，iter重置为0 print ("iteration number: %d" % iter) return b,alphas

　　使用如下代码运行上述简化版的SMO算法，算出SVM的b，和拉格朗日参数 (α1,α2,...,αm)组成的向量α，其中下标m=样本的数量.

if __name__=='__main__': dataMatIn,classLabels=loadDataSet('testSet.txt') b,alphas=smoSimple(dataMatIn,classLabels,0.6,0.001,40) print(b,alphas1[alphas1>0])

　　运行结果如下图所示：　　

完整版的SMO算法处理小规模数据集　　

　　完整版的SMO算法的相关细节请查看支持向量机原理（四）SMO算法原理。完整版的SMO算法实现如下代码：

''' class optStruct:不使用核函数的版本。结构体，里面存有数据集及其标签、还存有松弛因子、惩罚系数、所有的拉格朗日系数alphas、所有的参数Ei ''' class optStruct: def __init__(self,dataMatIn, classLabels, C, toler): # 使用参数初始化结构体 self.X = dataMatIn self.labelMat = classLabels self.C = C self.tol = toler self.m = shape(dataMatIn)[0] self.alphas = mat(zeros((self.m,1))) self.b = 0 self.eCache = mat(zeros((self.m,2))) #第一列是有效标志位，例如第k行的值为[1,Ek]，则表明Ek已经计算过了。当[0,Ek]，表示Ek还没有计算 ''' def calcE(oS, k):不使用核函数的版本。使用结构体oS中的数据，计算Ek ''' def calcEk(oS, k): fXk = float(multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T)) + oS.b Ek = fXk - float(oS.labelMat[k]) return Ek ''' def selectJ(i, oS, Ei):不使用核函数的版本。给定ai和Ei，从结构体oS中的数据集中，找出能够使得|Ei-Ej|最大的aj，并返回下标j和Ej ''' def selectJ(i, oS, Ei): maxK = -1; maxDeltaE = 0; Ej = 0 oS.eCache[i] = [1,Ei] #设置标志位为‘1’，表示该Ei已经计算好了 # #set valid #choose the alpha that gives the maximum delta E validEcacheList = nonzero(oS.eCache[:,0].A)[0] #返回array类型数组中非零值的下标 #oS.eCache[:,0].A其中‘.A’的作用是将列向量oS.eCache[:,0]转化为array类型 # 例如，b1 = array([True, False, True, False]) # print(nonzero(b1)[0]),返回值是非零元素的下标，本句的结果是[0 2] if (len(validEcacheList)) > 1: for j in validEcacheList: if j == i: continue #j==i，继续寻找 Ek = calcEk(oS, j) deltaE = abs(Ei - Ek) if (deltaE > maxDeltaE): maxJ = j; maxDeltaE = deltaE; Ej = Ek return maxK, Ej else: #第一轮oS.eCache[i]的所有标志位都为'0'，因为还没有开始计算 j = selectJrand(i, oS.m) Ej = calcEk(oS, j) return j, Ej ''' def updateEk(oS, k):不使用核函数的版本。计算Ek，并添加入结构体oS中的oS.eCache[k] ''' def updateEk(oS, k):#after any alpha has changed update the new value in the cache Ek = calcEk(oS, k) oS.eCache[k] = [1,Ek] ''' def innerL(i, oS):不使用核函数的版本。该函数给定参数ai的下标i和包含数据集的结构体oS，寻找出合适的参数aj并算出ai存入结构体oS中，并计算相应的参数b，最后返回参数0或1 （1表示成功找到合适的参数对ai-aj）。函数innerL(i,oS)完全和函数smoSimple(dataMatIn, classLabels, C, toler, maxIter)一样，唯一的区别就是函数innerL(i,oS)使用了自己的数据结构oS、选择j时使用函数selectJ(i, oS, Ei) ''' def innerL(i, oS): Ei = calcEk(oS, i) 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)): j,Ej = selectJ(i, oS, Ei) #this has been changed from selectJrand alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy(); if (oS.labelMat[i] != oS.labelMat[j]): L = max(0, oS.alphas[j] - oS.alphas[i]) H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i]) else: L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C) H = min(oS.C, oS.alphas[j] + oS.alphas[i]) if L==H: print( "L==H"); return 0 eta =oS.X[i,:]*oS.X[i,:].T + oS.X[j,:]*oS.X[j,:].T-2.0 * oS.X[i,:]*oS.X[j,:].T if eta <= 0: print ("eta<=0"); return 0 oS.alphas[j] += oS.labelMat[j]*(Ei - Ej)/eta oS.alphas[j] = clipAlpha(oS.alphas[j],H,L) updateEk(oS, j) #added this for the Ecache if (abs(oS.alphas[j] - alphaJold) < 0.00001): print ("j not moving enough"); return 0 oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j updateEk(oS, i) #added this for the Ecache #the update is in the oppostie direction 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,:].T 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,:].T if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1 elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2 else: oS.b = (b1 + b2)/2.0 return 1 else: return 0 ''' def smoP(dataMatIn, classLabels, C, toler, maxIter): 该函数是不使用核函数的完整版SMO。参数ai首先在所有的alphas中遍历来寻找aj，然后参数ai再在非边界值的alphas中遍历来寻找aj，这两个过程交替进行。该函数最后返回结构体oS中的b和向量alphas ''' def smoP(dataMatIn, classLabels, C, toler, maxIter): oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler) iter = 0 #变量iter统计两种情况，如果遍历所有的ai就iter+=1，或者遍历所有的非边界值ai就iter+=1 entireSet = True; alphaPairsChanged = 0 while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)): alphaPairsChanged = 0 if entireSet: #当entireSet==True，就遍历所有的alphas作为ai(i是从0到m-1) for i in range(oS.m): alphaPairsChanged += innerL(i,oS) print ("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)) iter += 1 else:#遍历所有的非边界值的alphas作为ai nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0] #返回oS.alphas列向量中，非边界值的alphas下标 #nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]返回oS.alphas列向量中，非零元素值且在区间(0,C)的元素下标。 #例如a=array([-1,2,3]) # print(a>0) # print(a<2.5) # print((a>0)*(a<2.5)) # print(nonzero((a>0)*(a<2.5))) # print(nonzero((a>0)*(a<2.5))[0]) #结果是 #[False True True] #[ True True False] #[False True False] #(array([1], dtype=int64),) #[1] for i in nonBoundIs: #nonBoundIs是非边界值的alphas下标，alphas下标与x的下标相等 alphaPairsChanged += innerL(i,oS) print ("non-bound, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)) iter += 1 if entireSet: entireSet = False #只要本轮是entireSet循环，无论是否找到参数对ai-aj，都关闭entireSet的循环 elif (alphaPairsChanged == 0): entireSet = True #无论本轮是entireSet循环还是非边界值循环，只要本轮找不到合适的参数对ai-aj，就打开entireSet循环 print ("iteration number: %d" % iter) return oS.b,oS.alphas

　　还是使用简化版SMO算法使用的数据集，我们运行以下代码，看看完整版的SMO算法效果如何：

if __name__=='__main__': dataMatIn,classLabels=loadDataSet('testSet.txt') b,alphas=smoP(dataMatIn,classLabels,0.6,0.001,40) print(b,alphas[alphas>0])

我们可以明显看出来，简化版SMO算法需要迭代40次，完整版SMO算法只需要迭代5次就能出结果。

在复杂数据上应用核函数　

　　这里我们使用另一个数据集“testSetRBF.txt”来进行加上核函数的完整版SMO算法。本例子使用径向基核函数(RBF)作为核函数，其表达式为：

k(x,y)=exp(−||x−y||22σ2)(1) (1) k ( x , y ) = e x p ( − | | x − y | | 2 2 σ 2 ) 　　代码实现的每一部分都已经标上注释，且在关键语句也给出了测试代码的注释来帮助理解，具体代码如下：

''' def kernelTrans(X, A, kTup):核函数。 X是数据集矩阵，A是某个样本Xi的行向量，kTup('type',参数delta)是控制核函数类型的元组。该函数根据不同的参数kTup，以不同的方式计算数据集矩阵X与样本A的径向基函数结果。例如X有m行（每行是一个样本的行向量），A是某个样本Xi的行向量，则该函数返回一个列向量[K1i，K2i，...,Kmi] 如果kTup[0]=='lin'，则执行线性核。如果kTup[0]=='rbf'，则执行径向基函数核其中kTup[1]是参数delta ''' def kernelTrans(X, A, kTup): #calc the kernel or transform data to a higher dimensional space m,n = shape(X) K = mat(zeros((m,1))) if kTup[0]=='lin': K = X * A.T #元组kTup是lin的话，就执行线性核 elif kTup[0]=='rbf': #元组是径向基函数rbf的话，就执行rbf核 for j in range(m): deltaRow = X[j,:] - A #行向量X[j,:]和行向量A对应元素相减 K[j] = deltaRow*deltaRow.T #计算|xj-y|^2，xj和y分别表示两个样本 K = exp(K/(-1*kTup[1]**2)) #divide in NumPy is element-wise not matrix like Matlab else: raise NameError('Houston We Have a Problem -- That Kernel is not recognized') return K ''' class optStructK:使用核函数的版本。这是结构体，里面存有数据集及其标签、还存有松弛因子、惩罚系数、所有的拉格朗日系数alphas、所有的参数Ei。相比没有使用核函数的版本，里面新增一个核函数矩阵。 ''' class optStructK: def __init__(self,dataMatIn, classLabels, C, toler, kTup): # Initialize the structure with the parameters self.X = dataMatIn self.labelMat = classLabels self.C = C self.tol = toler self.m = shape(dataMatIn)[0] self.alphas = mat(zeros((self.m,1))) self.b = 0 self.eCache = mat(zeros((self.m,2))) #first column is valid flag self.K = mat(zeros((self.m,self.m))) for i in range(self.m): #构建核函数矩阵 self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup) ''' def calcEkK(oS, k):使用核函数的版本。从结构体oS中的数据集中，计算Ek ''' def calcEkK(oS, k): fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b) Ek = fXk - float(oS.labelMat[k]) return Ek ''' def selectJK(i, oS, Ei):使用核函数的版本。给定ai和Ei，从结构体oS中的数据集中，找出能够使得|Ei-Ej|最大的aj，并返回下标j和Ej ''' def selectJK(i, oS, Ei): maxK = -1; maxDeltaE = 0; Ej = 0 oS.eCache[i] = [1,Ei] #设置标志位为‘1’，表示该Ei已经计算好了 # #set valid #choose the alpha that gives the maximum delta E validEcacheList = nonzero(oS.eCache[:,0].A)[0] #返回array类型数组中非零值的下标 #oS.eCache[:,0].A其中‘.A’的作用是将列向量oS.eCache[:,0]转化为array类型 # 例如，b1 = array([True, False, True, False]) # print(nonzero(b1)[0]),返回值是非零元素的下标，本句的结果是[0 2] if (len(validEcacheList)) > 1: for j in validEcacheList: if j == i: continue #j==i，继续寻找 Ek = calcEkK(oS, j) deltaE = abs(Ei - Ek) if (deltaE > maxDeltaE): maxJ = j; maxDeltaE = deltaE; Ej = Ek return maxK, Ej else: #第一轮oS.eCache[i]的所有标志位都为'0'，因为还没有开始计算 j = selectJrand(i, oS.m) Ej = calcEkK(oS, j) return j, Ej ''' def updateEk(oS, k):使用核函数的版本。计算Ek，并将Ek相应的存入 oS.eCache[k] ''' def updateEkK(oS, k):#after any alpha has changed update the new value in the cache Ek = calcEkK(oS, k) oS.eCache[k] = [1,Ek] ''' def innerLK(i, oS):使用核函数的版本。该函数给定参数ai的下标i和包含数据集的结构体oS，寻找出合适的参数aj并算出ai存入结构体oS中，并计算相应的参数b，最后返回参数0或1 （1表示成功找到合适的参数对ai-aj）。函数innerL(i,oS)完全和函数smoSimple(dataMatIn, classLabels, C, toler, maxIter)一样，唯一的区别就是函数innerL(i,oS)使用了自己的数据结构oS、选择j时使用函数selectJ(i, oS, Ei) ''' def innerLK(i, oS): Ei = calcEkK(oS, i) 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)): j,Ej = selectJK(i, oS, Ei) #this has been changed from selectJrand alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy(); if (oS.labelMat[i] != oS.labelMat[j]): L = max(0, oS.alphas[j] - oS.alphas[i]) H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i]) else: L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C) H = min(oS.C, oS.alphas[j] + oS.alphas[i]) if L==H: print ("L==H"); return 0 eta = oS.K[i,i] + oS.K[j,j]-2.0 * oS.K[i,j] #changed for kernel if eta <= 0: print ("eta<=0"); return 0 oS.alphas[j] += oS.labelMat[j]*(Ei - Ej)/eta oS.alphas[j] = clipAlpha(oS.alphas[j],H,L) updateEkK(oS, j) #added this for the Ecache if (abs(oS.alphas[j] - alphaJold) < 0.00001): print ("j not moving enough"); return 0 oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j updateEkK(oS, i) #added this for the Ecache #the update is in the oppostie direction 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] 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] if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1 elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2 else: oS.b = (b1 + b2)/2.0 return 1 else: return 0 ''' def smoPK(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):该函数是使用核函数的完整版SMO。参数ai首先在所有的alphas中遍历来寻找aj，然后参数ai再在非边界值的alphas中遍历来寻找aj，这两个过程交替进行。该函数最后返回结构体oS中的b和向量alphas ''' def smoPK(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)): #full Platt SMO oS = optStructK(mat(dataMatIn),mat(classLabels).transpose(),C,toler, kTup) iter = 0 entireSet = True; alphaPairsChanged = 0 while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)): alphaPairsChanged = 0 if entireSet: #go over all for i in range(oS.m): alphaPairsChanged += innerLK(i,oS) print ("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)) iter += 1 else:#go over non-bound (railed) alphas nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0] for i in nonBoundIs: alphaPairsChanged += innerLK(i,oS) print ("non-bound, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)) iter += 1 if entireSet: entireSet = False #toggle entire set loop elif (alphaPairsChanged == 0): entireSet = True print ("iteration number: %d" % iter) return oS.b,oS.alphas ''' def calcWs(alphas,dataArr,classLabels):由数组alphas、数据集dataArr、标签列表classLabels，计算出参数w。 ''' def calcWs(alphas,dataArr,classLabels): X = mat(dataArr); labelMat = mat(classLabels).transpose() m,n = shape(X) w = zeros((n,1)) # w是n行一列的向量，n是一个样本的属性数字 for i in range(m): w += multiply(alphas[i]*labelMat[i],X[i,:].T) #对于数组而言，multiply()和*一样，而对于矩阵而言，multiply是对应元素相乘，而*是矩阵乘法 # 这里的乘法*是对于数组而言，第i号元素alphas[i]*labelMat[i]的结果是一个数字 #X[i,:]是第i个样本Xi的所有属性组成的向量 #例如， Xi=array([[1,2,3]]) # print(Xi) # print(Xi.T) # print(multiply(3,Xi.T)) #结果是： # [[1 2 3]] # [[1] # [2] # [3]] # [[3] # [6] # [9]] return w ''' def testRbf(delta=1.3):delta是径向基核函数的参数。该函数将'testSetRBF.txt'作为训练集，然后得出b,alphas。然后给出在训练集'testSetRBF.txt'的分类错误率。然后使用训练集得到的分割超平面，在测试集'testSetRBF2.txt'中计算分类错误率。 ''' def testRbf(delta=1.3): dataArr,labelArr = loadDataSet('testSetRBF.txt') b,alphas = smoPK(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', delta)) #C=200 important datMat=mat(dataArr); labelMat = mat(labelArr).transpose() svInd=nonzero(alphas.A>0)[0] #获取支持向量的下标 sVs=datMat[svInd] #创建只有支持向量的矩阵 labelSV = labelMat[svInd]; print ("there are %d Support Vectors" % shape(sVs)[0]) m,n = shape(datMat) errorCount = 0 for i in range(m): kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', delta)) predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b if sign(predict)!=sign(labelArr[i]): errorCount += 1 print ("the training error rate is: %f" % (float(errorCount)/m)) dataArr,labelArr = loadDataSet('testSetRBF2.txt') errorCount = 0 datMat=mat(dataArr); labelMat = mat(labelArr).transpose() m,n = shape(datMat) for i in range(m): kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', delta))#计算支持向量与数据集的核函数结果 predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b if sign(predict)!=sign(labelArr[i]): errorCount += 1 print ("the test error rate is: %f" % (float(errorCount)/m))

　　运行如下代码并查看结果：

if __name__=='__main__': testRbf()

在手写数据集应用核函数　　

　　下面是读取‘手写数据集’的代码，并调用SVM算法进行测试：

''' def img2vector(filename):filename为32*32的图片，将32*32位数据按行为单位，摆成一个1024位的列向量 ''' def img2vector(filename): returnVect = zeros((1,1024)) fr = open(filename) for i in range(32): lineStr = fr.readline() for j in range(32): returnVect[0,32*i+j] = int(lineStr[j]) return returnVect ''' def loadImages(dirName):获取文件夹dirName下的所有手写数据集文件数据，最后返回数据集矩阵trainingMat（每行有1024位，代表一个32*32的图片），类别列表hwLabels（数字'9'的标签是+1，其余数字的标签是-1） ''' def loadImages(dirName): from os import listdir hwLabels = [] trainingFileList = listdir(dirName) #将文件夹dirName下的所有文件路径记录在列表trainingFileList m = len(trainingFileList) trainingMat = zeros((m,1024)) for i in range(m): fileNameStr = trainingFileList[i] fileStr = fileNameStr.split('.')[0] #take off .txt classNumStr = int(fileStr.split('_')[0]) #获取该文件的标签 if classNumStr == 9: hwLabels.append(-1) #如果是数字'9',则归类为-1 else: hwLabels.append(1) #如果数字不等于'9'，则归类为 1 trainingMat[i,:] = img2vector('%s/%s' % (dirName, fileNameStr)) return trainingMat, hwLabels ''' def testDigits(kTup=('rbf', 10)):使用带有核函数版本的SVM对手写数据集进行分类 ''' def testDigits(kTup=('rbf', 100)): dataArr,labelArr = loadImages('trainingDigits') b,alphas = smoPK(dataArr, labelArr, 200, 0.0001, 10000, kTup) datMat=mat(dataArr); labelMat = mat(labelArr).transpose() svInd=nonzero(alphas.A>0)[0] sVs=datMat[svInd] labelSV = labelMat[svInd]; print ("there are %d Support Vectors" % shape(sVs)[0]) m,n = shape(datMat) errorCount = 0 for i in range(m): kernelEval = kernelTrans(sVs,datMat[i,:],kTup) predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b if sign(predict)!=sign(labelArr[i]): errorCount += 1 print ("the training error rate is: %f" % (float(errorCount)/m)) dataArr,labelArr = loadImages('testDigits') errorCount = 0 datMat=mat(dataArr); labelMat = mat(labelArr).transpose() m,n = shape(datMat) for i in range(m): kernelEval = kernelTrans(sVs,datMat[i,:],kTup) predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b if sign(predict)!=sign(labelArr[i]): errorCount += 1 print ("the test error rate is: %f" % (float(errorCount)/m) )

　　下图是使用线性核函数对手写数据集验证的代码和结果：　　

转载请注明原文地址: https://www.6miu.com/read-2614270.html

技术

最新回复(0)