通过前面两篇文章的学习,这里我们简单的用python代码实现SVM,通过前面两篇文章,我们知道在求解SVM模型之前我们需要知道alpha(拉格朗日乘子),然后计算w和b,这样才能得到模型,其相应的公式如下:
计算alpha:
计算w:
计算b:
(1)定义存储所有的数据的数据结构
'''
功能:计算核函数
输入:X-整个数据集
A-第i行数据
kTup-核信息元组
输出:k-核矩阵
'''
def kernelTrans(X,A,kTup):
m,n=shape(X)
k=mat(zeros((m,1)))
if kTup[0]=='lin':
k=x*A.T
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('Houston we Have a problem--that kernel is not recognized')
return k;
class optStructK:
def __init__(self,dataMatIn,classLabels,C,toler,kTup):
self.X=dataMatIn#样本
self.labelMat=classLabels#标签
self.C=C#alpha的上边界
self.tol=toler#容错率
self.m=shape(dataMatIn)[0]#样本的数量
self.alphas=mat(zeros((self.m,1)))#alphas-拉格朗日乘子
self.b=0#模型的b值
self.eCache=mat(zeros((self.m,2)))#误差缓存
self.K=mat(zeros((self.m,self.m)))#核函数矩阵
for i in range(self.m):
self.K[:,i]=kernelTrans(self.X,self.X[i,:],kTup)
(2)SMO优化
'''
功能:计算E值
输入:oS-输入数据
k-第k个样本
输出: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
'''
功能:用于调整alphas的大小
'''
def clipAlpha(aj,H,L):
if aj>H:
aj=H
if aj<L:
aj=L
return aj
'''
功能:更新Ek
'''
def updateEk(oS,k):
Ek=calcEk(oS,k)
oS.eCache[k]=[1,Ek]
'''
功能:内循环中启发式方法选择第二个alphas的下标
输入:i-第一个alphas的下标
oS-输入数据
Ei-第i个样本的误差
输出:maxK-具有最大变化的第二个alphas的下标
Ej-第二个alphas的预测值与标签的差值
'''
def selectJ(i,oS,Ei):
#初始化
maxK=-1;maxDeltaE=0;Ej=0
oS.eCache[i]=[1,Ei]
#将变化量转化为非零表
validEcacheList=nonzero(oS.eCache[:,0].A)[0]
#如果存在第二个alphas
if(len(validEcacheList))>1:
#遍历所有alphas使得alphas变化最大的下标
for k in validEcacheList:
if k==i:continue
Ek=calcEk(oS,k)
deltaE=abs(Ei-Ek)
#选择具有最大步长的j
if(deltaE>maxDeltaE):
maxK=k;maxDeltaE=deltaE;Ej=Ek
return maxK,Ej
#如果不存在alphas,则采用随机选取法
else:
j=selectJrand(i,oS.m)
Ej=calcEk(oS,j)
return j,Ej
'''
功能:随机选择第二个alpha j的下标
输入:i-第一个alphas的下标
m-样本数量
输出:j-第二个alphas的下标
'''
def selectJrand(i,m):
j=i
while(j==i):
j=int(random.uniform(0,m))
return j
'''
功能:计算E值
输入:oS-输入数据
k-第k个样本
输出: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
'''
功能:计算两个变量的二次规划
输入:i-第一个变量alphas
oS-输入数据
'''
def innerLK(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)#使用启发式方法选取第二个alphas
alphaIold=oS.alphas[i].copy();alphaJold=oS.alphas[j].copy()#初始化两个alphas
#计算第二个alphas的剪切边界
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
#计算2K12-K11-K22
eta=2.0*oS.K[i,j]-oS.K[i,i]-oS.K[j,j]
if eta>=0:
print 'eta>=0';return 0
#计算第二个alphas,并对其进行约束
oS.alphas[j]-=oS.labelMat[j]*(Ei-Ej)/eta
oS.alphas[j]=clipAlpha(oS.alphas[j],H,L)
updateEk(oS,j)
if(abs(oS.alphas[j]-alphaJold)<0.00001):
print 'j not moving enought';return 0
#计算第一个alphas
oS.alphas[i]+=oS.labelMat[j]*oS.labelMat[i]*(alphaJold-oS.alphas[j])
updateEk(oS,i)
#计算b
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]
#如果0<a1<C,则b=b1
if(0<oS.alphas[i])and(oS.C>oS.alphas[i]):
oS.b=b1
#如果0<a2<C,则b=b2
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
'''
功能:SMO算法
输入:dataMatIn-输入样本
classLabels-输入标签
C-alphas的上边界
toler-误差
maxIter-最大迭代次数
输出:b-b
alphas-alphas
'''
def smoPK(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):
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
(3)使用SVM
'''
功能:计算w
'''
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 testRbf(k=1.3):
dataArr,labelArr=loadDataSet('./testSetRBF.txt')
b,alphas=smo.smoPK(dataArr,labelArr,200,0.001,10000,('rbf',k))
dataMat=mat(dataArr);labelMat=mat(labelArr).transpose()
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 = smo.kernelTrans(sVs,dataMat[i,:],('rbf', k))
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)
b,alphas=smo.smoPK(dataArr,labelArr,200,0.001,10000,('rbf',k))
dataMat=mat(dataArr);labelMat=mat(labelArr).transpose()
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 = smo.kernelTrans(sVs,dataMat[i,:],('rbf', k))
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)
dataArr1,labelArr1 = loadDataSet('testSetRBF2.txt')
errorCount = 0
datMat=mat(dataArr1); labelMat = mat(labelArr1).transpose()
m,n = shape(datMat)
for i in range(m):
kernelEval = smo.kernelTrans(sVs,datMat[i,:],('rbf', k))
predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
if sign(predict)!=sign(labelArr1[i]): errorCount += 1
print "the test error rate is: %f" % (float(errorCount)/m)
#显示数据集
dataMat0 = array(dataArr)
n = shape(dataMat0)[0]
xcord1 = []; ycord1 = []
xcord2 = []; ycord2 = []
for i in range(n):
if int(labelArr[i])== 1:
xcord1.append(dataMat0[i,0]); ycord1.append(dataMat0[i,1])
else:
xcord2.append(dataMat0[i,0]); ycord2.append(dataMat0[i,1])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
ax.scatter(xcord2, ycord2, s=30, c='green')
for i in range(n):
if alphas[i]>0.0:
circle = Circle((float(dataArr[i][0]), float(dataArr[i][1])), 0.1, facecolor='blue', \
edgecolor=(0,0.8,0.8), linewidth=3, alpha=0.5)
print dataArr[i],labelArr[i]
ax.add_patch(circle)
plt.show()