基于最大间隔分隔数据
import matplotlib import matplotlib.pyplot as plt from numpy import * xcord0 = [] ycord0 = [] xcord1 = [] ycord1 = [] markers =[] colors =[] fr = open('F:\machinelearninginaction\Ch06\testSet.txt')#this file was generated by 2normalGen.py for line in fr.readlines(): lineSplit = line.strip().split(' ') xPt = float(lineSplit[0]) yPt = float(lineSplit[1]) label = int(lineSplit[2]) if (label == 0): xcord0.append(xPt) ycord0.append(yPt) else: xcord1.append(xPt) ycord1.append(yPt) fr.close() fig = plt.figure() ax = fig.add_subplot(221) xcord0 = []; ycord0 = []; xcord1 = []; ycord1 = [] for i in range(300): [x,y] = random.uniform(0,1,2) if ((x > 0.5) and (y < 0.5)) or ((x < 0.5) and (y > 0.5)): xcord0.append(x); ycord0.append(y) else: xcord1.append(x); ycord1.append(y) ax.scatter(xcord0,ycord0, marker='s', s=90) ax.scatter(xcord1,ycord1, marker='o', s=50, c='red') plt.title('A') ax = fig.add_subplot(222) xcord0 = random.standard_normal(150); ycord0 = random.standard_normal(150) xcord1 = random.standard_normal(150)+2.0; ycord1 = random.standard_normal(150)+2.0 ax.scatter(xcord0,ycord0, marker='s', s=90) ax.scatter(xcord1,ycord1, marker='o', s=50, c='red') plt.title('B') ax = fig.add_subplot(223) xcord0 = [] ycord0 = [] xcord1 = [] ycord1 = [] for i in range(300): [x,y] = random.uniform(0,1,2) if (x > 0.5): xcord0.append(x*cos(2.0*pi*y)); ycord0.append(x*sin(2.0*pi*y)) else: xcord1.append(x*cos(2.0*pi*y)); ycord1.append(x*sin(2.0*pi*y)) ax.scatter(xcord0,ycord0, marker='s', s=90) ax.scatter(xcord1,ycord1, marker='o', s=50, c='red') plt.title('C') ax = fig.add_subplot(224) xcord1 = zeros(150); ycord1 = zeros(150) xcord0 = random.uniform(-3,3,350); ycord0 = random.uniform(-3,3,350); xcord1[0:50] = 0.3*random.standard_normal(50)+2.0; ycord1[0:50] = 0.3*random.standard_normal(50)+2.0 xcord1[50:100] = 0.3*random.standard_normal(50)-2.0; ycord1[50:100] = 0.3*random.standard_normal(50)-3.0 xcord1[100:150] = 0.3*random.standard_normal(50)+1.0; ycord1[100:150] = 0.3*random.standard_normal(50) ax.scatter(xcord0,ycord0, marker='s', s=90) ax.scatter(xcord1,ycord1, marker='o', s=50, c='red') plt.title('D') plt.show()
寻找最大间隔
分类器求解的优化问题
这里的类别标签为什么采用-1和+1，而不是0和 1呢？这是由于-1和+1仅仅相差一个符号，
方便数学上的处理。我们可以通过一个统一公式来表示间隔或者数据点到分隔超平面的距离，同
时不必担心数据到底是属于-1还是+1类。
S V M 应用的一般框架
from numpy import * from time import sleep def loadDataSet(fileName): dataMat = []; labelMat = [] fr = open(fileName) for line in fr.readlines(): lineArr = line.strip().split(' ') dataMat.append([float(lineArr[0]), float(lineArr[1])]) labelMat.append(float(lineArr[2])) return dataMat,labelMat 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): if aj > H: aj = H if L > aj: aj = L return aj
dataMat,labelMat = loadDataSet('F:\machinelearninginaction\Ch06\testSet.txt') print(labelMat)
可以看得出来，这里采用的类别标签是-1和1
def smoSimple(dataMatIn, classLabels, C, toler, maxIter): dataMatrix = mat(dataMatIn); labelMat = mat(classLabels).transpose() b = 0; m,n = shape(dataMatrix) alphas = mat(zeros((m,1))) iter = 0 while (iter < maxIter): alphaPairsChanged = 0 for i in range(m): fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b Ei = fXi - float(labelMat[i])#if checks if an example violates KKT conditions 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: print("L==H") continue eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T if eta >= 0: print("eta>=0") continue alphas[j] -= labelMat[j]*(Ei - Ej)/eta 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])#update i by the same amount as j #the update is in the oppostie direction 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 print("iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)) if (alphaPairsChanged == 0): iter += 1 else: iter = 0 print("iteration number: %d" % iter) return b,alphas
该函数有5个输人参数，分别
是：数据集、类别标签、常数C 、容错率和取消前最大的循环次数。
b,alphas = smoSimple(dataMat,labelMat, 0.6, 0.001, 40) print(b) print(alphas[alphas>0])
在原始数据集上对这些支持向量画圈之后的结果
import matplotlib import matplotlib.pyplot as plt from numpy import * from matplotlib.patches import Circle xcord0 = [] ycord0 = [] xcord1 = [] ycord1 = [] markers =[] colors =[] fr = open('F:\machinelearninginaction\Ch06\testSet.txt')#this file was generated by 2normalGen.py for line in fr.readlines(): lineSplit = line.strip().split(' ') xPt = float(lineSplit[0]) yPt = float(lineSplit[1]) label = int(lineSplit[2]) if (label == -1): xcord0.append(xPt) ycord0.append(yPt) else: xcord1.append(xPt) ycord1.append(yPt) fr.close() fig = plt.figure() ax = fig.add_subplot(111) ax.scatter(xcord0,ycord0, marker='s', s=90) ax.scatter(xcord1,ycord1, marker='o', s=50, c='red') plt.title('Support Vectors Circled') circle = Circle((4.6581910000000004, 3.507396), 0.5, facecolor='none', edgecolor=(0,0.8,0.8), linewidth=3, alpha=0.5) ax.add_patch(circle) circle = Circle((3.4570959999999999, -0.082215999999999997), 0.5, facecolor='none', edgecolor=(0,0.8,0.8), linewidth=3, alpha=0.5) ax.add_patch(circle) circle = Circle((6.0805730000000002, 0.41888599999999998), 0.5, facecolor='none', edgecolor=(0,0.8,0.8), linewidth=3, alpha=0.5) ax.add_patch(circle) #plt.plot([2.3,8.5], [-6,6]) #seperating hyperplane b = -3.75567 w0=0.8065 w1=-0.2761 x = arange(-2.0, 12.0, 0.1) y = (-w0*x - b)/w1 ax.plot(x,y) ax.axis([-2,12,-8,6]) plt.show()
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 #linear kernel elif kTup[0]=='rbf': for j in range(m): deltaRow = X[j,:] - A K[j] = deltaRow*deltaRow.T 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 optStruct: 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 calcEk(oS, k): fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b) Ek = fXk - float(oS.labelMat[k]) return Ek def selectJ(i, oS, Ei): #this is the second choice -heurstic, and calcs Ej maxK = -1; maxDeltaE = 0; Ej = 0 oS.eCache[i] = [1,Ei] #set valid #choose the alpha that gives the maximum delta E validEcacheList = nonzero(oS.eCache[:,0].A)[0] if (len(validEcacheList)) > 1: for k in validEcacheList: #loop through valid Ecache values and find the one that maximizes delta E if k == i: continue #don't calc for i, waste of time Ek = calcEk(oS, k) deltaE = abs(Ei - Ek) if (deltaE > maxDeltaE): maxK = k; maxDeltaE = deltaE; Ej = Ek return maxK, Ej else: #in this case (first time around) we don't have any valid eCache values j = selectJrand(i, oS.m) Ej = calcEk(oS, j) return j, Ej 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): 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 = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,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) 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.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 smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)): #full Platt SMO oS = optStruct(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 += innerL(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 += innerL(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
b,alphas = smoP(dataMat,labelMat, 0.6, 0.001, 40) print(b) print(alphas[alphas>0])
下面列出的一个小函数可
以用于实现上述任务：
def calcWs(alphas,dataArr,classLabels): X = mat(dataArr); labelMat = mat(classLabels).transpose() m,n = shape(X) w = zeros((n,1)) for i in range(m): w += multiply(alphas[i]*labelMat[i],X[i,:].T) return w
w = calcWs(alphas,dataMat,labelMat) print(w)
r = mat(dataMat[0])*w+b print(r)
print(labelMat[0])
利用核函数将数据映射到高维空间
数据点处于一个圆中，人类的大脑能够意识到这一点。然而，对于分类器而言，
它只能识别分类器的结果是大于0还是小于0。如果只在1和^轴构成的坐标系中插人直线进行分类
的话，我们并不会得到理想的结果。我们或许可以对圆中的数据进行某种形式的转换，从而得到
某些新的变量来表示数据。在这种表示情况下，我们就更容易得到大于0或者小于0的测试结果。
在这个例子中，我们将数据从一个特征空间转换到另一个特征空间。在新空间下，我们可以很容
易利用巳有的工具对数据进行处理。数学家们喜欢将这个过程称之为从一个特征空间到另一个特
择空间的映射。在通常情况下，这种映射会将低维特征空间映射到高维空间。
这种从某个特征空间到另一个特征空间的映射是通过核函数来实现的。
核函数并不仅仅应用于支持向量机，很多其他的机器学习算法也都用到核函数。
径向基核函数
上述高斯核函数将数据从其特征空间映射到更高维的空间，具体来说这里是映射到一个无穷
维的空间。
高斯核函数只是一个常用的核函数，使用
者并不需要确切地理解数据到底是如何表现的，而且使用高斯核函数还会得到一个理想的结果。
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 #linear kernel elif kTup[0]=='rbf': for j in range(m): deltaRow = X[j,:] - A K[j] = deltaRow*deltaRow.T 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
在测试中使用核函数
利用核函数进行分类的径向基测试函数
def testRbf(k1=1.3): dataArr,labelArr = loadDataSet('F:\machinelearninginaction\Ch06\testSetRBF.txt') b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) #C=200 important datMat=mat(dataArr); labelMat = mat(labelArr).transpose() svInd=nonzero(alphas.A>0)[0] sVs=datMat[svInd] #get matrix of only support vectors 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', k1)) 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('F:\machinelearninginaction\Ch06\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', k1)) 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))
