第六章 支持向量机
1.算法初识
假设我们有一个简单的二维数据集,其中有两类数据点,一类是红色圆点,另一类是蓝色菱形。我们的任务是找到一个能够很好地区分这两类数据的决策边界。在这个例子中,我们可以使用支持向量机(SVM)来找到这个决策边界
首先,我们可以将以下的数据集可视化,如下所示:
红圆点: (3, 3), (4, 3), (4, 4), (5, 5)
蓝菱形: (1, 2), (2, 1), (2, 3), (3, 2)
接下来,我们需要找到最优的超平面来分隔这两个类别。通过计算,我们可以得到一个线性决策边界(一条直线),如
y
=
−
x
+
4
y = -x + 4
y=−x+4这个超平面能够最大化两个类别数据点之间的距离,即决策边界与最近的数据点(支持向量)之间的距离。在这个例子中,这些支持向量是红圆点 (3, 3) (4, 3)和蓝菱形 (3, 2)(2, 3)
在这个简单的例子中,SVM 通过找到一个线性决策边界成功将两类数据点分开。对于非线性可分的问题,可以通过核技巧将数据映射到高维空间,使得数据在高维空间中线性可分
2.基于最大间隔分离数据
基本概念
- 线性可分数据:存在一个超平面(或决策边界),可以将不同类别的数据完全分隔开来的数据集
- 分隔超平面:指在分类问题中,将不同类别的数据点分开的一个超平面。在二维空间中,分割超平面是一条直线;在三维空间中,它是一个平面;在更高维的空间中,它是一个超平面
- 间隔:某一个数据点到分隔面的距离
- 支持向量:用于确定分割超平面的训练样本点。它们是距离分割超平面最近的训练样本点
SVM简述
在一组线性可分数据中,能用一条直线将数据集分隔开,随着维数增大直线就变成了超平面,通过最大化支持向量到分隔面的距离来确定我们要寻找数据集的最佳分隔直线,SVM通过寻找最优的分隔超平面,并利用支持向量来确定分隔超平面的位置和方向,从而实现对线性可分数据的最佳分隔
3.寻找最大间隔
按照上述的例子,记分隔超平面为 l l l,它可写成 w T x + b w^Tx+b wTx+b,某一点 A A A到 l l l的法线长度为 ∣ w T A + b ∣ / ∣ ∣ w ∣ ∣ |w^TA+b|/||w|| ∣wTA+b∣/∣∣w∣∣其中 b b b类似于逻辑回归中的截距
分类器优化
SVM分类器的工作原理是将输入数据输出为一个类别标签, 比如使用类似单位阶跃函数对 w T x + b w^Tx+b wTx+b作用得到 u = f ( w T x + b ) u=f(w^Tx+b) u=f(wTx+b),当 u < 0 u<0 u<0时输出为1,反之为-1
数据点到分隔面的间隔通过 l a b e l ∗ ( w T x + b ) label*(w^Tx+b) label∗(wTx+b)来计算,因此不论是1还是-1, l a b e l ⋅ ( w T x + b ) label·(w^Tx+b) label⋅(wTx+b)都会是一个正数方便计算。当前的目标就是找到分类器定义中的 w w w和 b b b,也就是找到具有最小间隔的数据点,即支持向量,而后进行最大化,写作 arg m a x w , b { m i n n ( l a b e l ⋅ ( w T x + b ) ) ⋅ 1 ∣ ∣ w ∣ ∣ } \arg max_{w,b}\{min_n(label·(w^Tx+b))·\frac{1}{||w||}\} argmaxw,b{minn(label⋅(wTx+b))⋅∣∣w∣∣1}由于这个计算复杂因此给定一些约束条件然后求最优解,可以通过拉格朗日乘子法求解,这里的约束条件为 l a b e l ∗ ( w T x + b ) ≥ 1.0 label*(w^Tx+b)\geq1.0 label∗(wTx+b)≥1.0,因此原式可以写成 m a x a [ ∑ i = 1 m α − 1 2 ∑ i , j = 1 m l a b e l ( i ) ⋅ l a b e l ( j ) ⋅ α i ⋅ α j < x ( i ) ⋅ x ( j ) > ] max_a\left[\sum_{i=1}^m \alpha-\frac{1}{2}\sum_{i,j=1}^mlabel^{(i)}·label^{(j)}·\alpha_i·\alpha_j\left<x^{(i)}·x^{(j)}\right>\right] maxa[i=1∑mα−21i,j=1∑mlabel(i)⋅label(j)⋅αi⋅αj⟨x(i)⋅x(j)⟩]其中约束条件是 α ≥ 0 , 和 ∑ i , j = 1 m α i ⋅ l a b e l ( i ) = 0 \alpha\ge0,和\sum_{i,j=1}^m\alpha_i·label^{(i)}=0 α≥0,和i,j=1∑mαi⋅label(i)=0但是以上要满足数据都可线性分隔,实际并不满足,此时引入松弛变量,约束条件变成 C ≥ α ≥ 0 , 和 ∑ i , j = 1 m α i ⋅ l a b e l ( i ) = 0 C\ge\alpha\ge0,和\sum_{i,j=1}^m\alpha_i·label^{(i)}=0 C≥α≥0,和i,j=1∑mαi⋅label(i)=0,其中的 C C C是常数参数控制“最大化间隔”和“保证多数点的最大化间隔小于1”
4.SMO算法
本章关注SVM中的序列最小优化算法(SMO),它是将大优化问题分解为多个小优化问题,SMO算法的目标是求出一系列的 α \alpha α和 b b b,从而计算权重向量w并得到分隔超平面 l l l,它的工作原理是在每次循环中选择满足条件的两个 α \alpha α,进行优化处理,增大其中一个减小另一个。需要满足的条件有俩:1.两个 α \alpha α在间隔边界外2.两个 α \alpha α没有进行区间化处理或不在边界上
简化版SMO处理
由于改变一个 α \alpha α时可能会导致约束条件失效,因此要同时改变两个 α \alpha α,构造以下辅助函数用于在某个区间范围内随机选择一个整数,以及在 α \alpha α数值过大时进行调整
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 # 我们要选择与 i 不相等的任意 J
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
loadDataSet
:从文件中加载数据集。函数读取文件的每一行,将特征存储在dataMat
列表中,将标签存储在labelMat
列表中,并返回这两个列表
selectJrand
:随机选择一个不等于i的整数j。该函数用于选择在优化算法中的两个变量alpha的索引,以确保每次选择的alpha对不同
clipAlpha
:对
a
l
p
h
a
alpha
alpha进行修剪,确保其值在指定的上界H和下界L之间。该函数用于在优化算法中更新
a
l
p
h
a
alpha
alpha时保持其在合理范围内
可以看到结果是1或-1而不是0或1,接下来使用SMO算法
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 # 记录每轮更新的alpha对数量
for i in range(m):
fXi = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[i,:].T)) + b # 预测结果
Ei = fXi - float(labelMat[i]) # 预测结果与实际标签之间的误差
# 检查是否违反KKT条件
if ((labelMat[i] * Ei < -toler) and (alphas[i] < C)) or ((labelMat[i] * Ei > toler) and (alphas[i] > 0)):
j = selectJrand(i, m) # 随机选择另一个alpha
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])
# 更新偏置项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]):
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
以上代码是一个简化版的SMO算法,用于求解SVM的二分类问题。该算法使用拉格朗日乘子法进行优化,通过迭代更新优化变量来寻找最优的超平面
实际运行效果如下
将得到支持向量进行标记并在原数据集中进行可视化效果如下
完整Platt SMO加速优化
完整SMO与简化的SMO最主要的差别就是选择 α \alpha α的方式不同,Platt SMO算法是通过一个外循环来选择第一个 α \alpha α,要么是在所有数据集上单边扫描,或者是在非边界 α \alpha α 中实现单边扫描,非边界 α \alpha α* 是指不等于边界0或C的 α \alpha α值,这时候就要创建一个非边界 α \alpha α列表方便遍历扫描
然后是通过内循环最大化步长获得第二个 α \alpha α值,建立全局存缓保存误差值,选择最大步长的 α \alpha α值,以下是对原SMO进行处理并得到的完整Platt SMO代码
def kernelTrans(X, A, kTup): #calc the kernel or transform data to a higher dimensional space
m,n = shape(X)
K = mat(zeros((m,1)))
if kTup[0]=='lin': K = X * A.T #linear kernel
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)) #divide in NumPy is element-wise not matrix like Matlab
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): # Initialize the structure with the parameters
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))) #first column is valid flag
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):
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): #this is the second choice -heurstic, and calcs Ej
maxK = -1; maxDeltaE = 0; Ej = 0
oS.eCache[i] = [1,Ei] #set valid #choose the alpha that gives the maximum delta E
validEcacheList = nonzero(oS.eCache[:,0].A)[0]
if (len(validEcacheList)) > 1:
for k in validEcacheList: #loop through valid Ecache values and find the one that maximizes delta E
if k == i: continue #don't calc for i, waste of time
Ek = calcEk(oS, k)
deltaE = abs(Ei - Ek)
if (deltaE > maxDeltaE):
maxK = k; maxDeltaE = deltaE; Ej = Ek
return maxK, Ej
else: #in this case (first time around) we don't have any valid eCache values
j = selectJrand(i, oS.m)
Ej = calcEk(oS, j)
return j, Ej
def updateEk(oS, k):#after any alpha has changed update the new value in the cache
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
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
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)): #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)):
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
结果如上,与之前的结果相比支持向量的数量变多了,现在已经得到
α
\alpha
α值,接下来就是利用它们得到超平面
def calcWs(alphas, dataArr, classLabels):
X = mat(dataArr) # 转换为矩阵形式
labelMat = mat(classLabels).transpose() # 转换为列向量的矩阵形式
m, n = shape(X) # 获取数据矩阵的行数和列数
w = zeros((n, 1)) # 初始化权重矩阵 w,全为零向量
for i in range(m): # 遍历每个样本
w += multiply(alphas[i] * labelMat[i], X[i, :].T) # 根据支持向量的拉格朗日乘子和对应类别标签,更新权重向量 w
return w # 返回最终权重向量 w
上面的calcWs
函数用于计算在支持向量机训练过程中得到的权重向量 w,是基于支持向量的拉格朗日乘子和对应的类别标签,通过对每个支持向量进行加权累加来得到的
利用上述的代码能够对数据进行分类,比如上面的这个值<0,就把它归为标签为-1的
检验正确
5.应用核函数
一般像下面这种非可线性的数据集不可能用一条直线直接进行二分
但是通过核函数就能将某种不易分隔的数据转换成易于分类器理解的形式
核函数映射高维空间
在上述图像中,一般难以用直线进行数据分隔,但是通过某种转换就能将原数据用新的变量来表示,且这些变量能更轻易被分隔,也就是将数据从一个特征空间转换到另一个特征空间,即不同空间的映射,而实现这个映射的工具就是核函数,可以想象成一个空间间的接口比如上图通过高斯径向基核函数就能得到下图
能看到基本上两类数据被一个椭圆形给分开了,这就是核函数的魅力。在SVM中所有的运算都能写成内积形式,将内积替换成核函数的方法称为核技巧,下面就介绍一个核函数:径向基核函数
径向基核函数
径向基核函数以向量为自变量,基于距离运算输出一个标量。径向基核函数的高斯版公式为 k ( x , y ) = e x p ( − ∣ ∣ x − y ∣ ∣ 2 2 σ 2 ) k(x,y)=exp\left(\frac{-||x-y||^2}{2\sigma^2}\right) k(x,y)=exp(2σ2−∣∣x−y∣∣2)其中 σ \sigma σ用于确定到达率,即函数值降到0的速度,下述是核函数的核转换函数代码
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)) # 使用高斯径向基核函数(RBF kernel)
else:
raise NameError('未识别的核函数类型')
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) # 计算核函数值
函数kernelTrans
用于计算核函数转换后的数据。该函数接受三个参数:X是输入数据矩阵,A是某个样本点,kTup是核函数的类型和参数。函数根据给定的核函数类型,计算输入数据矩阵X与样本点A之间的核函数值,并返回结果。函数optStruct
用于保存SVM算法的优化参数和数据结构
为使用以上的核函数代码,要对之前的innerL
和calcEK
函数做修改,以下是修改后的代码
def innerL(i, oS):
Ei = calcEk(oS, i) # 计算第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)):
# 如果样本i的预测误差超出容错率范围,并且拉格朗日乘子alphas[i]的取值满足一定条件
j, Ej = selectJ(i, oS, Ei) # 选择另一个样本j以形成最大步长的对偶变量对
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] # 计算步长eta
if eta >= 0:
print("eta>=0")
return 0
oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej) / eta # 更新拉格朗日乘子alphas[j]
oS.alphas[j] = clipAlpha(oS.alphas[j], H, L) # 根据约束条件裁剪alphas[j]的取值
updateEk(oS, j) # 更新样本j的预测误差缓存
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]) # 更新拉格朗日乘子alphas[i]
updateEk(oS, i) # 更新样本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.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 # 更新决策函数的偏置项b
return 1
else:
return 0
def calcEk(oS, k):
fXk = float(multiply(oS.alphas, oS.labelMat).T * oS.K[:, k] + oS.b) # 计算样本k的预测值
Ek = fXk - float(oS.labelMat[k]) # 计算样本k的预测误差
return Ek
使用核函数
之前已经给出一个非线性可分数据的例子,现在对另一个非线性可分数据集进行核函数分类,代码如下
def testRbf(k1=1.3):
# 加载数据集
dataArr, labelArr = loadDataSet('testSetRBF.txt')
# 使用SMO算法训练SVM模型
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]
sVs = datMat[svInd] # 获取支持向量矩阵
labelSV = labelMat[svInd]
print("支持向量的个数:%d" % 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("训练集上的错误率:%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("测试集上的错误率:%f" % (float(errorCount) / m))
实现效果如上,自上而下分别是自定义参数k=1.8、1.3和0.1时的结果,0.1中每个支持向量的影响度少因此需要较多数量, 易产生过拟合情况,而比较1.8和1.3的情况可以知道1.3的参数更加合适,支持向量数量少且错误率也不差
6.SMO核函数示例
以KNN的手写识别数据集为样本进行SVM的使用,过程如下
- 收集数据:提供的文本文件。
- 准备数据:基于二值图像构造向量。
- 分析数据:对图像向量进行目测。
- 训练算法:采用两种不同的核函数,并对径向基核函数采用不同的设置来进行SMO算法。
- 测试算法:编写一个函数来测试不同的核函数并计算错误率。
- 使用算法:字面意思
相关的代码如下
def loadImages(dirName):
from os import listdir
hwLabels = []
trainingFileList = listdir(dirName) #load the training set
m = len(trainingFileList)
trainingMat = zeros((m,1024))
for i in range(m):
fileNameStr = trainingFileList[i]
fileStr = fileNameStr.split('.')[0] #take off .txt
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)):
dataArr,labelArr = loadImages('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];
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('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))
结果如下,可以通过修改其他参数达到较好的效果,但是由于跑不动就只做了一个