# coding=utf-8
from numpy import *
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
#加载数据集
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
#i为alpha_1在样本集合中的下标 ,需要在找一个不同的alpha_2 其下表不能为i
def selectJrand(i,m):
j=i
while(j==i):
j=int(random.uniform(0,m))
return j
#alpha需要满足一定的约束条件 0<=alpha<=C sum(alpha*y)=0
def clipAlpha(aj,H,L):
if aj>H:
aj=H
if L>aj:
aj=L
return aj
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
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)))#生成一个全是0的alphas向量
self.b=0
self.eCache=mat(zeros((self.m,2)))#误差缓存 eCache第一列给出eCache是否有效的标志位,第二列给出的是实际值
def calcEk(oS,k):#计算给定alpha的时候的误差 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):
#选择第二个/内循环的alpha值 选择合适的第二个alphas值 是的每次有优化的时候的步长都最大
maxK=-1
maxDeltaE=0
Ej=0
oS.eCache[i]=[1,Ei]#先设置为有效的
validEcacheList=nonzero(oS.eCache[:,0].A)[0]
if (len(validEcacheList))>1:
for k in validEcacheList:
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]
#SMO内部函数
def innerL(i,oS):
Ei=calcEk(oS,i)#计算第一个alpha的误差
# 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)):
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)):
#SMO必须要满足KKT条件
j,Ej=selectJ(i,oS,Ei)#选择具有最大步长的第二个alpha
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])
eta=2.0*oS.X[i,:]*oS.X[j,:].T-oS.X[i,:]*oS.X[i,:].T-oS.X[j,:]*oS.X[j,:].T
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):
return 0
oS.alphas[i]+=oS.labelMat[j]*oS.labelMat[i]*(alphaJold-oS.alphas[j])
updateEk(oS,i)
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)#oS.K[i,j]
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)#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
#选择第一个alpha的策略,在所有的数据集上进行单边扫描,在非边界上进行扫描
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)):# entireSet控制在边界上和整个数据集上进行选择alpha
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,oS.C
dataArr,labelArr=loadDataSet("testSet.txt")
b,alphas,C=smoP(dataArr,labelArr,0.6,0.01,40)
ws=calcWs(alphas,dataArr,labelArr)
print ws
from matplotlib.patches import Circle
xcord0=[]
ycord0=[]
xcord1=[]
ycord1=[]
markers=[]
colors=[]
fr=open("testSet.txt")
for line in fr.readlines():
lineSplit = line.strip().split('\t')
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()
xd=[]
yd=[]
m=shape(dataArr)[0]
for i in range(100):
if alphas[i]>0.0 and alphas[i]!=C:
xd.append(dataArr[i])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord0,ycord0, marker='s', s=50)
ax.scatter(xcord1,ycord1, marker='o', s=50, c='red')
plt.title('SVM')
totalm=shape(xd)[0];
for i in range(totalm):
circle = Circle((xd[i][0], xd[i][1]), 0.3, facecolor='none', edgecolor=(0,0.8,0.8), linewidth=2, alpha=0.5)
ax.add_patch(circle)
w0=float(ws[0])
w1=float(ws[1])
x=arange(-2.0,12.0,0.1)
x=mat(x)
z=(-w0*x-b)/w1
z=mat(z)
ax.plot(x.tolist()[0],z.tolist()[0])
ax.axis([-2,12,-8,6])
plt.show()
速度特别快