整理并分析了里面一些错误和不适之处,已经在代码中修改,原因是Python版本的更新,这次代码比较全,有学习相关章节的读者可以参考下面的代码,博主学习时,已经测试过,如有问题,请留言交流哈。
# -*-coding:utf-8 -*-
import random
from numpy import *
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):
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
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
#以上测试,没有问题。
class optStruct:
def __init__(self,dataMatIn,classLabels,C,toler,kTup):
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)))
#changed for kernel
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):
# changed for kernel
fXk=float(multiply(oS.alphas,oS.labelMat).T*(oS.K[:,k]+oS.b))
#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):
maxK=-1
maxDeltaE=0
Ej=0
oS.eCache[i]=[1,Ei]
validEcaheList=nonzero(oS.eCache[:,0].A)[0]
if (len(validEcaheList))>1:
for k in validEcaheList:
if k==i:continue
Ek=calcEk(oS,k)
deltaE=abs(Ei-Ek)
if(deltaE>maxDeltaE):
maxK=k
maxDeltaE=deltaE
Ej=Ek
return maxK,Ej
else:
j=selectJrand(i,oS.m)
Ej=calcEk(oS,j)
return j,Ej
def updateEk(oS,k):
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
#eta=2.0 * oS.X[i,:]*oS.X[j,:].T - oS.X[i,:]*oS.X[i,:].T - oS.X[j,:]*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
#changed for kernel
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]
#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)):
oS=optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler,kTup)
iter=0
entireSet=True
alphaParisChanged=0
while (iter< maxIter) and ((alphaParisChanged>0) or (entireSet)):
alphaParisChanged=0
if entireSet:
for i in range(oS.m):
alphaParisChanged+=innerL(i,oS)
print("fullset,iter: %d i: %d, pairs changed %d" %(iter,i,alphaParisChanged))
iter+=1
else:
nonBoundIs=nonzero((oS.alphas.A>0)* (oS.alphas.A <C))[0]
for i in nonBoundIs:
alphaParisChanged+=innerL(i,oS)
print("non-bound,iter: %d i:%d,pairs changed %d" %(iter,i,alphaParisChanged))
iter+=1
if entireSet:
entireSet=False
elif (alphaParisChanged==0):
entireSet=True
print("iteration number: %d" %iter)
return oS.b,oS.alphas
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
def kernelTrans(X,A,kTup):
'''
The first argument in the tuple is a string describing what type of kernel should be used.
the other arguments are optional arguments that may be needed for a kernel.
In the case of the linear kernel, a dot product is taken between the two inputs,
which are the full dataset and a row of the dataset.
In the case of the radial bias function,the Gaussian function is evaluated for every element
in the matrix in the for loop. After the for loop is finished, you apply the calculations over the entire vector.
'''
#kTup tuple contain the kernel function information:sgma.
m,n=shape(X)
K=mat(zeros((m,1)))
#linear kernel
if kTup[0]=='lin':
K=X*A.T
#radial bias function
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))
else:
raise NameError('Houstion We Have a Problem -- \
That Kernel is not recognized')
return K
def testRbf(k1=1.3):
dataArr,labelArr=loadDataSet('data/testSetRBF.txt')
b,alphas=smoP(dataArr,labelArr,200,0.001,10000,('rbf',k1))
dataMat=mat(dataArr)
labelMat=mat(labelArr).transpose()
#martix.A means transform martix into array
svInd=nonzero(alphas.A>0)[0]
sVs=dataMat[svInd]
labelSV=labelMat[svInd]
print("There are %d Support Vectors" %shape(sVs)[0])
m,n=shape(dataMat)
errorCount=0
for i in range(m):
kernelEval=kernelTrans(sVs,dataMat[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('data/testSetRBF2.txt')
errorCount=0
dataMat=mat(dataArr)
labelMat=mat(labelArr).transpose()
m,n=shape(dataMat)
for i in range(m):
kernelEval=kernelTrans(sVs,dataMat[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))
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):
from os import listdir
hwlabels=[]
trainingFileList=listdir(dirName)
m=len(trainingFileList)
trainingMat=zeros((m,1024))
for i in range(m):
fileNameStr=trainingFileList[i]
fileStr=fileNameStr.split('.')[0]
classNumStr=int(fileStr.split('_')[0])
if classNumStr==9:
hwlabels.append(-1)
else:
hwlabels.append(1)
trainingMat[i,:]=img2vector('%s/%s' %(dirName,fileNameStr))
return trainingMat,hwlabels
def testDigits(kTup=('rbf',10)):
#same as testRbf.
dataArr,labelArr=loadImages('digits/trainingDigits')
b,alphas=smoP(dataArr,labelArr,200,0.0001,10000,kTup)
dataMat=mat(dataArr)
labelMat=mat(labelArr).transpose()
#martix.A means transform martix into array
svInd=nonzero(alphas.A>0)[0]
sVs=dataMat[svInd]
labelSV=labelMat[svInd]
print("There are %d Support Vectors" %shape(sVs)[0])
m,n=shape(dataMat)
errorCount=0
for i in range(m):
kernelEval=kernelTrans(sVs,dataMat[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('digits/testDigits')
errorCount = 0
dataMat = mat(dataArr)
labelMat = mat(labelArr).transpose()
m, n = shape(dataMat)
for i in range(m):
kernelEval = kernelTrans(sVs, dataMat[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))
'''
if __name__=='__main__':
dataArr,labelArr=loadDataSet('data/testSet.txt')
b,alphas=smoP(dataArr,labelArr,0.6,0.001,40)
ws=calcWs(alphas,dataArr,labelArr)
testRbf()
testDigits()'''