支持向量机

#!/usr/bin/python
# -*- coding: utf-8 -*-
#coding=utf-8

from numpy import *
from time import sleep

#加载数据
#打开文件并逐行解析，从而得到每行的类标签和整个数据矩阵
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的j值
def selectJrand(i, m):
j = i
while(j == i):
j = int(random.uniform(0, m))
return j

#调整aj
def clipAlpha(aj, H, L):
if aj > H:
aj = H
if aj < L:
aj = L
return aj

#核转换函数，将低维空间的数据映射到高维空间
#采用径向基函数的高斯版本
#输入参数为2个数值型变量和1个元祖。元祖中第一个参数为核函数类型，其它2个为核函数可能需要的可选参数
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

#kTup是包含核函数信息的元祖
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)))  #误差缓存，第一列是是否合法的标志位
self.K = mat(zeros((self.m, self.m)))  #矩阵K被构建
for i in range(self.m):              #调用函数kernelTrans()进行填充
self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)

#计算E值并返回
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

#选择第二个alpha或者内循环的alpha值
def selectJ(i, oS, Ei):
maxK = -1
maxDeltaE = 0
Ej = 0
oS.eCache[i] = [1, Ei]
validEcacheList = nonzero(oS.eCache[:,0].A)[0]  #合法E值的下标
#寻找具有最大步长的j
if (len(validEcacheList)) > 1:
for k in validEcacheList:
if k == i:
continue
Ek = calcEk(oS, k)
deltaE = abs(Ei - Ek)
if(deltaE > maxDeltaE):
maxDeltaE = deltaE
maxK = k
Ej = Ek
return maxK, Ej
else:  #如果是第一次运行，eCache中无合法值
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)
#如果误差较大且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)):
#采用启发式方法选择第二个alpha
j, Ej = selectJ(i, oS, Ei)
alphaIold = oS.alphas[i].copy()
alphaJold = oS.alphas[j].copy()
#保证alpha在0和C之间
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是alphas[j]的最优修改量
eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j]
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)  #更新误差缓存
if(abs(oS.alphas[j] - alphaJold) < 0.00001):
print "j not moving enough"
return 0
#i的修改量与j相同，但方向相反
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.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.alphas[i] < oS.C):
oS.b = b1
elif (0 < oS.alphas[j]) and(oS.alphas[j] < oS.C):
oS.b = b2
else:
oS.b = (b1 + b2) / 2.0
return 1
else:
return 0

#SMO外循环代码
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
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:  #遍历所有值
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:
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
elif (alphaPairsChanged == 0):
entireSet = True
print "iteration number: %d" % iter
return oS.b, oS.alphas

#计算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

#!/usr/bin/python
# -*- coding: utf-8 -*-
#coding=utf-8

from numpy import *
from time import sleep

#加载数据
#打开文件并逐行解析，从而得到每行的类标签和整个数据矩阵
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的j值
def selectJrand(i, m):
j = i
while(j == i):
j = int(random.uniform(0, m))
return j

#调整aj
def clipAlpha(aj, H, L):
if aj > H:
aj = H
if aj < L:
aj = L
return aj

#核转换函数，将低维空间的数据映射到高维空间
#采用径向基函数的高斯版本
#输入参数为2个数值型变量和1个元祖。元祖中第一个参数为核函数类型，其它2个为核函数可能需要的可选参数
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

#kTup是包含核函数信息的元祖
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)))  #误差缓存，第一列是是否合法的标志位
self.K = mat(zeros((self.m, self.m)))  #矩阵K被构建
for i in range(self.m):              #调用函数kernelTrans()进行填充
self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)

#计算E值并返回
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

#选择第二个alpha或者内循环的alpha值
def selectJ(i, oS, Ei):
maxK = -1
maxDeltaE = 0
Ej = 0
oS.eCache[i] = [1, Ei]
validEcacheList = nonzero(oS.eCache[:,0].A)[0]  #合法E值的下标
#寻找具有最大步长的j
if (len(validEcacheList)) > 1:
for k in validEcacheList:
if k == i:
continue
Ek = calcEk(oS, k)
deltaE = abs(Ei - Ek)
if(deltaE > maxDeltaE):
maxDeltaE = deltaE
maxK = k
Ej = Ek
return maxK, Ej
else:  #如果是第一次运行，eCache中无合法值
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)
#如果误差较大且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)):
#采用启发式方法选择第二个alpha
j, Ej = selectJ(i, oS, Ei)
alphaIold = oS.alphas[i].copy()
alphaJold = oS.alphas[j].copy()
#保证alpha在0和C之间
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是alphas[j]的最优修改量
eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j]
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)  #更新误差缓存
if(abs(oS.alphas[j] - alphaJold) < 0.00001):
print "j not moving enough"
return 0
#i的修改量与j相同，但方向相反
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.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.alphas[i] < oS.C):
oS.b = b1
elif (0 < oS.alphas[j]) and(oS.alphas[j] < oS.C):
oS.b = b2
else:
oS.b = (b1 + b2) / 2.0
return 1
else:
return 0

#SMO外循环代码
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
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:  #遍历所有值
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:
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
elif (alphaPairsChanged == 0):
entireSet = True
print "iteration number: %d" % iter
return oS.b, oS.alphas

#计算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(k1=1.3):
dataArr,labelArr = loadDataSet('testSetRBF.txt')
b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) #C=200 重要
datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
svInd=nonzero(alphas.A>0)[0]  #找出非零的alpha值，从而得到所需要的支持向量
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', 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('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)

>>> import svmMLiA
>>> dataArr, labelArr = loadDataSet('testSet.txt')
>>> b, alphas = smoP(dataArr, labelArr, 0.6, 0.001, 40)
>>> ws = calcWs(alphas, dataArr, labelArr)
>>> dataMat = mat(dataArr)
>>> labelMat = mat(labelArr)
>>> dataMat[0] * mat(ws) + b
matrix([[-1.11362722]])   #小于0的应该是-1；大于0为1
>>> labelArr[0]
-1.0
#继续检查其它分类结果
>>> dataMat[2] * mat(ws) + b
matrix([[ 2.37754259]])
>>> labelArr[2]
1.0
>>> dataMat[1] * mat(ws) + b
matrix([[-1.58551706]])
>>> labelArr[1]
-1.0
>>>

>>> import svmMLiA
>>> testRbf()

... ...
L==H
L==H
L==H
L==H
fullSet, iter : 6 i : 99, pairs changed 0
iteration number: 7
there are 27 Support Vectors
the training error rate is: 0.030000
the test error rate is: 0.040000