支持向量(support vector)就是离分隔超平面最近的点
SVM特点:几乎所有分类问题都可以解决,SVM本身是一个二类分类器,对多类问题应用SVM需要对代码做一些修改
序列最小优化算法(Sequential Minimal Optimization,SMO)算法:
支持向量机的学习问题可以形式化为求解凸二次规划问题。
SMO算法将原来的大问题不断分解为子问题并对子问题求解,进而达到求解原问题的目的。
整个SMO算法包含两个部分:求解两个变量二次规划的解析方法和选择变量的启发方法。
先看代码:(简要的SMO实现代码)
# encoding: utf-8
import random
import numpy as np
from numpy import *
#6-1 SMO算法中的辅助函数
def loadDataSet(fileName):#读取txt中数据
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):#i是第一个alpha的下标,m是所有alpha的数目
j = i
while (j == i ):
j = int(random.uniform(0,m))
return j
def clipAlpha(aj,H,L):#调整大于H或者小于L的alpha值
if aj > H:
aj = H
if L > aj:
aj = L
return aj
#6-2简化版SMO算法
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
#输入的五个参数分别为数据集、类别标签、常数C、容错率和退出前最大的循环次数
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 #这个参数用于记录alphas参数是否改变,改变的话即为1
for i in range(m):
fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b#我们预测的类别
Ei = fXi - float(labelMat[i])#预测与实际之间的误差
#是否可以继续优化,判断依据是误差率与实际预测的积与容错率相比较
if ((labelMat[i]*Ei < -toler) and(alphas[i] < C)) or \
((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
j = selectJrand(i,m)#j是随机取得的index,是我们取得一个随机样本
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是alpha[j]的最优修改量
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]
alphas[j] -= labelMat[j]*(Ei - Ej)/eta
#通过门限函数来选取范围
alphas[j] = clipAlpha(alphas[j],H,L)
#检查alphas[j]是否有轻微的改变
if (abs(alphas[j] - alphaJold) < 0.00001):
print('j not moving enough');continue
#alphas[i]进行和alphas[j]同样的改变,但两者改变的方向相反,若一个增加,另一个就减小
alphas[i] += labelMat[j] * labelMat[i] * (alphaJold -alphas[j])
#为两个alphas位置常数b
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]):#在0到C的范围内
b = b1
elif (0 < alphas[j]) and (C > alphas[j]):#同样alphas[j]也这样判断
b = b2
else:b = (b1 + b2)/2.0
#如果到这都没有跳出这个循环,说明两个alphas已经修改过了
alphaPairsChanged += 1
print('iter: %d i:%d, pairs changed %d' %(iter,i,alphaPairsChanged))
#检查alpha是否做了更新,如果有更新则将iter归0
if (alphaPairsChanged == 0):iter += 1
else: iter = 0
print('iteration number: %d'% iter)
return b,alphas
一、两个变量二次规划的求解方法
我们首先选取一个alpha值 i,并计算它的预测值和误差。
g(x)就是我们的预测值,但不是实际的预测值,这里面肯定有误差
误差就是Ei = g(xi) - yi
在简化版里我们随机选择另一个alpha值 j,选之前判断一下第一个alpha值能否需要继续优化,不需要就不用往下走了。
同样的方法,我们再计算j对应的预测值和误差。
然后就是比较这两个alpha的预测值即y1和y2,确定最大范围H和最小范围L:
1.当y1 = y2时,
2.当y1 != y2时,
算出最优化分量η(eta),
是输入空间到特征空间的映射。
未经剪辑时的
剪辑后在L、H和 中选取。也同样要改变,数值上和alpha2一致,但方向相反。
同时更新b1和b2的值:
如果和同时满足那么,否则 b = (b1 + b2)/2.0(中点)
二、变量的选择方法
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)))
self.b = 0
self.eCache = mat(zeros((self.m,2)))
def calcEk(oS, k):#计算E值并返回(误差)
#fXk = float(multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T)) + oS.b
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):#用于选择第二个alpha或者说内循环的alpha值,选择其中使得改变最大的那个值
#第2个变量的选择标准是:希望能使alpha2有足够大的变化
maxK = -1;maxDeltaE = 0; Ej = 0
oS.eCache[i] = [1,Ei] #根据Ei更新误差缓存
#nonzero()返回一个array,但这里[0],使得返回的是非0元素的行号
validEcacheList = nonzero(oS.eCache[:,0].A)[0]#储存着非0元素索引
if (len(validEcacheList)) > 1:
for k in validEcacheList:
#通过for循环遍历选择其中使得改变最大的那个值
if k == i:continue #如果坐标与第一个alpha索引一样,退出
Ek = calcEk(oS,k) #计算出这个alpha值对应的误差
deltaE = abs(Ei - Ek) #两个误差的差值
if (deltaE > maxDeltaE):#判断语句作用就是寻找最大差值对应的alpha索引
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]
#6-4完整Platt SMO算法中的优化历程(寻找决策边界的优化)
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)#选择第二个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])
if L==H:print('L==H');return 0
# eta是alpha[j]的最优修改量
#eta = 2.0 *oS.X[i,:]*oS.X[j,:].T - oS.X[i,:]*oS.X[i,:].T - oS.X[j,:]*oS.X[j,:].T
eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j]
if eta >= 0:print('eta >= 0');return 0
# 计算出一个新的alphas[j]
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
# alphas[i]进行和alphas[j]同样的改变,但两者改变的方向相反,若一个增加,另一个就减小
oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])
updateEk(oS,i)
# 为两个alphas位置常数b
#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
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.X[i, :] * oS.X[j, :].T - oS.labelMat[j] * \
#(oS.alphas[j] - alphaJold) * oS.X[j, :] * oS.X[j, :].T
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
#6-5完整版Platt SMO的外循环代码
def smop(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin',0)):
#输入的五个参数分别为数据集、类别标签、常数C、容错率和退出前最大的循环次数
#构建一个数据结构来容纳所有数据
oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler,kTup)
iter = 0 #进行各项的初始化
entireSet = True; alphaPairsChanged = 0
#循环推出条件:1.迭代次数大于最大迭代次数 2.遍历整个集合都未对任意alpha进行修改
while (iter < maxIter) and ((alphaPairsChanged > 0)or (entireSet)):
alphaPairsChanged = 0
if entireSet:
for i in range (oS.m):#遍历alpha
alphaPairsChanged += innerL(i,oS)
#调用innerL()来选择第二个alpha,并在可能时对其进行优化处理。如果alpha改变返回1,没有改变,返回0
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
1.第1个变量的选择:SMO称选择第1个变量的过程为外层循环。外层循环在训练样本中选取违反KKT条件最严重的样本点
2.第2个变量的选择:SMO称选择第1个变量的过程为内层循环。第2个变量选择的标准是希望是alpha2有足够大的变化。
简化版:采取的仅仅是随机选取alpha2。而完整版里采用的是selectJ()函数来选取。
三、核函数来解决非线性分类
#6-6核转换函数
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 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)))
for i in range(self.m):
self.K[:, i] = kernelTrans(self.X, self.X[i, :], kTup)
#6-8利用核函数进行分类的径向基测试函数
def testRbf(k1=0.31):#该参数就是核函数中的σ的大小
dataArr,labelArr = loadDataSet('testSetRBF.txt')#读取数据集
b,alphas = smop(dataArr,labelArr,200,0.0001,10000,('rbf',k1))#使用SMO算法,其中核函数类型为rbf
datMat = mat(dataArr); labelMat = mat(labelArr).transpose()
svInd = nonzero(alphas.A > 0)[0]#找出非0的alpha值的索引
sVs = datMat[svInd] #找到支持向量
labelSV = labelMat[svInd] #alpha的类别标签值
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))
四、手写识别系统SVM来分类
#6-9基于SVM的手写数字识别
#为了使用前面两个分类器,因为前面的例子都是只有一行的,我们必须将图像格式化处理为一个向量。
# 我们将把一个32*32的二进制图像矩阵转换为1*1024的向量,这样前两节使用的分类器就可以处理数字图像信息了。
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)
#os.listdir() 方法用于返回指定的文件夹包含的文件或文件夹的名字的列表。
# 这个列表以字母顺序。 它不包括 '.' 和'..' 即使它在文件夹中。
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('digits/trainingDigits/%s' % fileNameStr)
return trainingMat,hwLabels
def testDigits(KTup = ('rbf',10)):
dataArr,labelArr = loadImages('digits/trainingDigits')#读取训练集数据
b,alphas = smop(dataArr,labelArr,200,0.0001,10000,KTup)
datMat = mat(dataArr);labelMat = mat(labelArr).transpose()
svInd = nonzero(alphas.A > 0)[0]
sVs = datMat[svInd] # 找到支持向量
labelSV = labelMat[svInd] # alpha的类别标签值
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,:], 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
datMat = mat(dataArr)
labelMat = mat(labelArr).transpose()
m, n = shape(datMat)
for i in range(m):
kernelEval = kernelTrans(sVs, datMat[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)
开始测试:
1.(rbf,100)时:
2.(rbf,10)时:
3.('rbf',0.1)时:
四、使用sklearn封装的库中的SVM
参考文章:https://blog.csdn.net/yuanlulu/article/details/81011849
https://www.cnblogs.com/pinard/p/6117515.html
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('digits/%s' %(dirName))
#os.listdir() 方法用于返回指定的文件夹包含的文件或文件夹的名字的列表。
# 这个列表以字母顺序。 它不包括 '.' 和'..' 即使它在文件夹中。
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:#类别标签只有-1和1,当等于9时,为-1,,否则则为1。
hwLabels.append(-1)
else:
hwLabels.append(1)
trainingMat[i,:] = img2vector('digits/%s/%s' % (dirName,fileNameStr))
return trainingMat,hwLabels
TraindataArr,TrainlabelArr = loadImages('trainingDigits')#读取训练集数据
dataArr, labelArr = loadImages('testDigits') # 读取测试集
#print (TrainlabelArr)
#clf = svm.SVC(C=1.0, kernel='rbf', degree=3, gamma='auto', coef0=0.0, shrinking=True, probability=False, tol=0.001, cache_size=200, class_weight=None, verbose=False, max_iter=-1, decision_function_shape='ovr', random_state=None)
thresholds = np.linspace(0,0.001,100)
param_grid = {'gamma': thresholds}
clf = GridSearchCV(svm.SVC(kernel='rbf'), param_grid, cv=5)
clf.fit(TraindataArr,TrainlabelArr)
print("best param: {0}\nbest score: {1}".format(clf.best_params_,clf.best_score_))
y_pred = clf.predict(dataArr)
print("y_pred:{}".format(y_pred))
print("y_test:{}".format(labelArr))
print("查准率:",metrics.precision_score(y_pred, labelArr))
print("召回率:",metrics.recall_score(y_pred, labelArr))
print("F1:",metrics.f1_score(y_pred, labelArr))
通过GridSearchCV函数(sklearn里面的)来寻找rbf模型中的参数最合适的情况。运行结果如下:
以及召回率、准确率等等评价参数:https://blog.csdn.net/quiet_girl/article/details/70830796