一、基于最大间隔分割数据
支持向量机 优点:泛化错误率低,计算开销不大,结果易解释 缺点:对参数调节和核函数的选择敏感,原始分类器不加修改仅适用于处理二类问题。 适用数据类型:数值型和标称型数据。
数据点在二维平面上,分隔超平面就只是一条直线,但数据集是三维时,那么分隔超平面就是一个平面。依此类推,如果数据集是
N
(
N
≥
2
N(N\geq2
N(N≥2)维时,那么就需要一个
N
?
1
N-1
N?1维的对象来分隔数据。该对象被称为超平面,也就是分类的决策边界。
理想状态是分布在超平面一侧的所有数据都属于某个类别,而分布在另一侧的所有数据则属于另一个类别。
间隔应该尽可能地大,这是因为如果我们犯错或者在有限数据上训练分类器的话,大的间隔可以增加分类器的鲁棒性。
支持向量到分割面的距离应该最大化。
二、寻找最大间隔
如何求解数据集的最佳分隔直线?以下图为例: 分隔超平面的形式为:
w
τ
x
+
b
{\bf w^\tau x} + b
wτx+b 点A到分隔超平面的法线长度:
∣
w
τ
A
+
b
∣
/
∣
∣
w
∣
∣
|{\bf w^\tau A} + b| / | | {\bf w} | |
∣wτA+b∣/∣∣w∣∣
最大化间隔的目标就是找出分类器定义中的w和b。为此,我们必须找到具有最小间隔的数据点,而这些数据点也就是前面提到的支持向量。一旦找到具有最小间隔的数据点,我们就需要对该间隔最大化。这就可以写作:
arg?max
?
w
,
b
?
{
min
?
n
?
(
l
a
b
e
l
?
(
w
τ
x
+
b
)
)
?
1
∣
∣
w
∣
∣
}
{\underset {w, b}{\operatorname {arg\,max} }}\, \left\{ {\underset {n}{\operatorname {min} }}\, (label · ({\bf w^\tau x} + b)) · {\frac{1}{| | {\bf w} | |}} \right\}
w,bargmax?{nmin?(label?(wτx+b))?∣∣w∣∣1?}
直接求解上述问题相当困难,所以我们将它转换成为另一种更容易求解的形式。
m
a
x
α
[
∑
i
=
1
m
α
?
1
2
∑
i
,
j
=
1
m
l
a
b
e
l
(
i
)
?
l
a
b
e
l
(
j
)
?
a
i
?
a
j
?
x
(
i
)
,
x
(
j
)
?
]
{\underset {\alpha}{max}}\left[ \sum^{m}_{i=1}\alpha - \frac{1}{2}\sum^{m}_{i,j=1}label^{(i)} · label^{(j)} · a_i · a_j \left\langle x^{(i)}, x^{(j)} \right\rangle \right]
αmax?[i=1∑m?α?21?i,j=1∑m?label(i)?label(j)?ai??aj??x(i),x(j)?]
约束条件为:
C
≥
α
≥
0
和
∑
i
=
1
m
α
i
?
l
a
b
e
l
(
i
)
=
0
C \geq \alpha \geq 0 和 \sum^m_{i=1}\alpha_i · label^{(i)} = 0
C≥α≥0和i=1∑m?αi??label(i)=0
其中常数C 用于控制 “最大化间隔” 和 “保证大部分点的函数间隔小于1.0” 这两个目标的权重。在优化算法的实现代码中,常数C 是一个参数,因此可以通过调节该参数得到不同的结果。一旦求出了所有的
α
\alpha
α,那么分隔超平面就可以通过这些
α
\alpha
α 来表达。
SVM的一般流程 (1) 收集数据:可以使用任意方法。 (2) 准备数据:需要数值型数据。 (3) 分析数据:有助于可视化分隔超平面。 (4) 训练算法:SVM的大部分时间都源自训练,该过程主要实现两个参数的调优。 (5) 测试算法:十分简单的计算过程就可以实现。 (6) 使用算法:几乎所有分类问题都可以使用SVM,值得一提的是,SVM本身是一个二类分类器,对多类问题应用SVM需要对代码做一些修改。
三、简化版SMO算法
简化版SMO算法,省略了确定要优化的最佳
α
\alpha
α 对的步骤,而是首先在数据集上进行遍历每一个
α
\alpha
α,再在剩余的数据集中找到另外一个
α
\alpha
α,构成要优化的
α
\alpha
α 对,同时对其进行优化,这里的同时是要确保公式:
∑
α
i
?
l
a
b
e
l
(
i
)
=
0
\sum\alpha_i * label^{(i)} = 0
∑αi??label(i)=0
所以改变一个
α
\alpha
α 显然会导致等式失效,所以这里需要同时改变两个
α
\alpha
α。代码如下:
from time import sleep
import matplotlib.pyplot as plt
import numpy as np
import random
import types
"""
函数说明:读取数据
Parameters:
fileName - 文件名
Returns:
dataMat - 数据矩阵
labelMat - 数据标签
"""
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
"""
函数说明:随机选择alpha
Parameters:
i - alpha
m - alpha参数个数
Returns:
j -
"""
def selectJrand(i, m):
j = i
while (j == i):
j = int(random.uniform(0, m))
return j
"""
函数说明:修剪alpha
Parameters:
aj - alpha值
H - alpha上限
L - alpha下限
Returns:
aj - alpah值
"""
def clipAlpha(aj,H,L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj
"""
函数说明:简化版SMO算法
Parameters:
dataMatIn - 数据矩阵
classLabels - 数据标签
C - 松弛变量
toler - 容错率
maxIter - 最大迭代次数
Returns:
无
"""
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
dataMatrix = np.mat(dataMatIn); labelMat = np.mat(classLabels).transpose()
b = 0; m,n = np.shape(dataMatrix)
alphas = np.mat(np.zeros((m,1)))
iter_num = 0
while (iter_num < maxIter):
alphaPairsChanged = 0
for i in range(m):
fXi = float(np.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)
fXj = float(np.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("alpha_j变化太小"); continue
alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
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("第%d次迭代 样本:%d, alpha优化次数:%d" % (iter_num,i,alphaPairsChanged))
if (alphaPairsChanged == 0): iter_num += 1
else: iter_num = 0
print("迭代次数: %d" % iter_num)
return b,alphas
"""
函数说明:分类结果可视化
Parameters:
dataMat - 数据矩阵
w - 直线法向量
b - 直线解决
Returns:
无
"""
def showClassifer(dataMat, w, b):
data_plus = []
data_minus = []
for i in range(len(dataMat)):
if labelMat[i] > 0:
data_plus.append(dataMat[i])
else:
data_minus.append(dataMat[i])
data_plus_np = np.array(data_plus)
data_minus_np = np.array(data_minus)
plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1], s=30, alpha=0.7)
plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1], s=30, alpha=0.7)
x1 = max(dataMat)[0]
x2 = min(dataMat)[0]
a1, a2 = w
b = float(b)
a1 = float(a1[0])
a2 = float(a2[0])
y1, y2 = (-b- a1*x1)/a2, (-b - a1*x2)/a2
plt.plot([x1, x2], [y1, y2])
for i, alpha in enumerate(alphas):
if abs(alpha) > 0:
x, y = dataMat[i]
plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
plt.show()
"""
函数说明:计算w
Parameters:
dataMat - 数据矩阵
labelMat - 数据标签
alphas - alphas值
Returns:
无
"""
def get_w(dataMat, labelMat, alphas):
alphas, dataMat, labelMat = np.array(alphas), np.array(dataMat), np.array(labelMat)
w = np.dot((np.tile(labelMat.reshape(1, -1).T, (1, 2)) * dataMat).T, alphas)
return w.tolist()
if __name__ == '__main__':
dataMat, labelMat = loadDataSet('testSet.txt')
b,alphas = smoSimple(dataMat, labelMat, 0.6, 0.001, 40)
w = get_w(dataMat, labelMat, alphas)
showClassifer(dataMat, w, b)
通过多次运行测试后,发现该简化SMO算法运行时间不算太短,而这只是一个仅有100个点的小规模数据集,这就意味着在更大的数据集上,收敛时间会变得更长,所以要将通过完整的SMO算法进行加速。
四、加速优化的完整版 Platt SMO 算法
在这两个版本(简化版和完整版)中,实现
α
\alpha
α 的更改和代数运算的优化环节一模一样。在优化过程中,唯一的不同就是 选择
α
\alpha
α的方式。完整版的Platt SMO算法应用了一些能够提速的启发方法。
Platt SMO算法是通过一个 外循环 来选择第一个
α
\alpha
α值的,并且其选择过程会在两种方式之间进行交替:一种方式是在所有数据集上进行单遍扫描,另一种方式则是在非边界
α
\alpha
α中实现单遍扫描。而所谓非边界
α
\alpha
α指的就是那些不等于边界0或C的
α
\alpha
α值。对整个数据集的扫描相当容易,而实现非边界
α
\alpha
α值的扫描时,首先需要建立这些
α
\alpha
α值的列表,然后再对这个表进行遍历。同时,该步骤会跳过那些已知的不会改变的
α
\alpha
α值。
在选择第一个
α
\alpha
α值后,算法会通过一个 内循环 来选择第二个
α
\alpha
α值。在优化过程中,会通过 最大化步长 的方式来获得第二个
α
\alpha
α值。在简化版SMO算法中,我们会在选择j 之后计算错误率
E
j
E_j
Ej?。但在这里,我们会建立一个全局的缓存用于保存误差值,并从中选择使得步长或者说
E
i
?
E
j
E_i - E_j
Ei??Ej? 最大的
α
\alpha
α值。
import matplotlib.pyplot as plt
import numpy as np
import random
"""
Parameters:
dataMatIn - 数据矩阵
classLabels - 数据标签
C - 松弛变量
toler - 容错率
"""
class optStruct:
def __init__(self, dataMatIn, classLabels, C, toler):
self.X = dataMatIn
self.labelMat = classLabels
self.C = C
self.tol = toler
self.m = np.shape(dataMatIn)[0]
self.alphas = np.mat(np.zeros((self.m,1)))
self.b = 0
self.eCache = np.mat(np.zeros((self.m,2)))
"""
Parameters:
fileName - 文件名
Returns:
dataMat - 数据矩阵
labelMat - 数据标签
"""
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
"""
Parameters:
oS - 数据结构
k - 标号为k的数据
Returns:
Ek - 标号为k的数据误差
"""
def calcEk(oS, k):
fXk = float(np.multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T) + oS.b)
Ek = fXk - float(oS.labelMat[k])
return Ek
"""
Parameters:
i - alpha_i的索引值
m - alpha参数个数
Returns:
j - alpha_j的索引值
"""
def selectJrand(i, m):
j = i
while (j == i):
j = int(random.uniform(0, m))
return j
"""
Parameters:
i - 标号为i的数据的索引值
oS - 数据结构
Ei - 标号为i的数据误差
Returns:
j, maxK - 标号为j或maxK的数据的索引值
Ej - 标号为j的数据误差
"""
def selectJ(i, oS, Ei):
maxK = -1; maxDeltaE = 0; Ej = 0
oS.eCache[i] = [1,Ei]
validEcacheList = np.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
"""
Parameters:
oS - 数据结构
k - 标号为k的数据的索引值
Returns:
无
"""
def updateEk(oS, k):
Ek = calcEk(oS, k)
oS.eCache[k] = [1,Ek]
"""
Parameters:
aj - alpha_j的值
H - alpha上限
L - alpha下限
Returns:
aj - 修剪后的alpah_j的值
"""
def clipAlpha(aj,H,L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj
"""
Parameters:
i - 标号为i的数据的索引值
oS - 数据结构
Returns:
1 - 有任意一对alpha值发生变化
0 - 没有任意一对alpha值发生变化或变化太小
"""
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)
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.X[i,:] * oS.X[j,:].T - oS.X[i,:] * oS.X[i,:].T - oS.X[j,:] * oS.X[j,:].T
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("alpha_j变化太小")
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
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
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
"""
Parameters:
dataMatIn - 数据矩阵
classLabels - 数据标签
C - 松弛变量
toler - 容错率
maxIter - 最大迭代次数
Returns:
oS.b - SMO算法计算的b
oS.alphas - SMO算法计算的alphas
"""
def smoP(dataMatIn, classLabels, C, toler, maxIter):
oS = optStruct(np.mat(dataMatIn), np.mat(classLabels).transpose(), C, toler)
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("全样本遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter,i,alphaPairsChanged))
iter += 1
else:
nonBoundIs = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
for i in nonBoundIs:
alphaPairsChanged += innerL(i,oS)
print("非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter,i,alphaPairsChanged))
iter += 1
if entireSet:
entireSet = False
elif (alphaPairsChanged == 0):
entireSet = True
print("迭代次数: %d" % iter)
return oS.b,oS.alphas
"""
Parameters:
dataMat - 数据矩阵
w - 直线法向量
b - 直线解决
Returns:
无
"""
def showClassifer(dataMat, classLabels, w, b):
data_plus = []
data_minus = []
for i in range(len(dataMat)):
if classLabels[i] > 0:
data_plus.append(dataMat[i])
else:
data_minus.append(dataMat[i])
data_plus_np = np.array(data_plus)
data_minus_np = np.array(data_minus)
plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1], s=30, alpha=0.7)
plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1], s=30, alpha=0.7)
x1 = max(dataMat)[0]
x2 = min(dataMat)[0]
a1, a2 = w
b = float(b)
a1 = float(a1[0])
a2 = float(a2[0])
y1, y2 = (-b- a1*x1)/a2, (-b - a1*x2)/a2
plt.plot([x1, x2], [y1, y2])
for i, alpha in enumerate(alphas):
if abs(alpha) > 0:
x, y = dataMat[i]
plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
plt.show()
"""
Parameters:
dataArr - 数据矩阵
classLabels - 数据标签
alphas - alphas值
Returns:
w - 计算得到的w
"""
def calcWs(alphas,dataArr,classLabels):
X = np.mat(dataArr); labelMat = np.mat(classLabels).transpose()
m,n = np.shape(X)
w = np.zeros((n,1))
for i in range(m):
w += np.multiply(alphas[i]*labelMat[i],X[i,:].T)
return w
if __name__ == '__main__':
dataArr, classLabels = loadDataSet('testSet.txt')
b, alphas = smoP(dataArr, classLabels, 0.6, 0.001, 40)
w = calcWs(alphas,dataArr, classLabels)
showClassifer(dataArr, classLabels, w, b)
优化后: 并且通过代码计时,发现优化后的程序运行时间相较于之前快了很多。
五、在复杂数据上应用核函数
核函数的目的主要是为了解决非线性分类问题,通过核技巧将低维的非线性特征转化为高维的线性特征,从而可以通过线性模型来解决非线性的分类问题。 在图中,数据点处于一个圆中,人类的大脑能够意识到这一点。然而,对于分类器而言,它只能识别分类器的结果是大于0还是小于0。如果只在x和y轴构成的坐标系中插入直线进行分类的话,我们并不会得到理想的结果。但是或许可以对圆中的数据进行某种形式的转换,从而得到某些新的变量来表示数据。在这种表示情况下,我们就更容易得到大于0或者小于0的测试结果。在通常情况下,这种映射是通过 核函数 来实现的,会将低维特征空间映射到高维空间
这种从某个特征空间到另一个特征空间的映射是通过核函数来实现的。它能把数据从某个很难处理的形式转换成为另一个较容易处理的形式。如果上述特征空间映射的说法听起来很让人迷糊的话,那么可以将它想象成为另外一种距离计算的方法。距离计算的方法有很多种,核函数一样具有多种类型。经过空间转换之后,我们可以在高维空间中解决线性问题,这也就等价于在低维空间中解决非线性问题。
将数据集一分为二,作为训练集和测试集进行测试。训练集如下: 可见,数据明显是线性不可分的。下面我们根据公式,编写核函数,并增加初始化参数kTup用于存储核函数有关的信息,同时我们只要将之前的内积运算变成核函数的运算即可。最后编写testRbf()函数,用于测试。
import matplotlib.pyplot as plt
import numpy as np
import random
"""
Parameters:
dataMatIn - 数据矩阵
classLabels - 数据标签
C - 松弛变量
toler - 容错率
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 = np.shape(dataMatIn)[0]
self.alphas = np.mat(np.zeros((self.m,1)))
self.b = 0
self.eCache = np.mat(np.zeros((self.m,2)))
self.K = np.mat(np.zeros((self.m,self.m)))
for i in range(self.m):
self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)
"""
Parameters:
X - 数据矩阵
A - 单个数据的向量
kTup - 包含核函数信息的元组
Returns:
K - 计算的核K
"""
def kernelTrans(X, A, kTup):
m,n = np.shape(X)
K = np.mat(np.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 = np.exp(K/(-1*kTup[1]**2))
else: raise NameError('核函数无法识别')
return K
"""
Parameters:
fileName - 文件名
Returns:
dataMat - 数据矩阵
labelMat - 数据标签
"""
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
"""
Parameters:
oS - 数据结构
k - 标号为k的数据
Returns:
Ek - 标号为k的数据误差
"""
def calcEk(oS, k):
fXk = float(np.multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
Ek = fXk - float(oS.labelMat[k])
return Ek
"""
Parameters:
i - alpha_i的索引值
m - alpha参数个数
Returns:
j - alpha_j的索引值
"""
def selectJrand(i, m):
j = i
while (j == i):
j = int(random.uniform(0, m))
return j
"""
Parameters:
i - 标号为i的数据的索引值
oS - 数据结构
Ei - 标号为i的数据误差
Returns:
j, maxK - 标号为j或maxK的数据的索引值
Ej - 标号为j的数据误差
"""
def selectJ(i, oS, Ei):
maxK = -1; maxDeltaE = 0; Ej = 0
oS.eCache[i] = [1,Ei]
validEcacheList = np.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
"""
Parameters:
oS - 数据结构
k - 标号为k的数据的索引值
Returns:
无
"""
def updateEk(oS, k):
Ek = calcEk(oS, k)
oS.eCache[k] = [1,Ek]
"""
Parameters:
aj - alpha_j的值
H - alpha上限
L - alpha下限
Returns:
aj - 修剪后的alpah_j的值
"""
def clipAlpha(aj,H,L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj
"""
Parameters:
i - 标号为i的数据的索引值
oS - 数据结构
Returns:
1 - 有任意一对alpha值发生变化
0 - 没有任意一对alpha值发生变化或变化太小
"""
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)
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]
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("alpha_j变化太小")
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.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
"""
Parameters:
dataMatIn - 数据矩阵
classLabels - 数据标签
C - 松弛变量
toler - 容错率
maxIter - 最大迭代次数
kTup - 包含核函数信息的元组
Returns:
oS.b - SMO算法计算的b
oS.alphas - SMO算法计算的alphas
"""
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup = ('lin',0)):
oS = optStruct(np.mat(dataMatIn), np.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("全样本遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter,i,alphaPairsChanged))
iter += 1
else:
nonBoundIs = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
for i in nonBoundIs:
alphaPairsChanged += innerL(i,oS)
print("非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter,i,alphaPairsChanged))
iter += 1
if entireSet:
entireSet = False
elif (alphaPairsChanged == 0):
entireSet = True
print("迭代次数: %d" % iter)
return oS.b,oS.alphas
"""
Parameters:
k1 - 使用高斯核函数的时候表示到达率
Returns:
无
"""
def testRbf(k1 = 0.3):
dataArr,labelArr = loadDataSet('testSetRBF.txt')
b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 100, ('rbf', k1))
datMat = np.mat(dataArr); labelMat = np.mat(labelArr).transpose()
svInd = np.nonzero(alphas.A > 0)[0]
sVs = datMat[svInd]
labelSV = labelMat[svInd];
print("支持向量个数:%d" % np.shape(sVs)[0])
m,n = np.shape(datMat)
errorCount = 0
for i in range(m):
kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
predict = kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b
if np.sign(predict) != np.sign(labelArr[i]): errorCount += 1
print("训练集错误率: %.2f%%" % ((float(errorCount)/m)*100))
dataArr,labelArr = loadDataSet('testSetRBF2.txt')
errorCount = 0
datMat = np.mat(dataArr); labelMat = np.mat(labelArr).transpose()
m,n = np.shape(datMat)
for i in range(m):
kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
predict=kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b
if np.sign(predict) != np.sign(labelArr[i]): errorCount += 1
print("测试集错误率: %.2f%%" % ((float(errorCount)/m)*100))
if __name__ == '__main__':
testRbf()
>>>
······
全样本遍历:第2次迭代 样本:99, alpha优化次数:0
迭代次数: 3
支持向量个数:29
训练集错误率: 0.00%
测试集错误率: 4.00%
可以尝试不同的K1值,经过实验会发现K1越大,过拟合越严重。支持向量的数目存在一个最优值。SVM的优点在于它能对数据进行高效分类。如果支持向量太少,就可能会得到一个很差的决策边界;如果支持向量太多,也就相当于每次都利用整个数据集进行分类,这种分类方法称为k近邻。
六、根据建立的气候模型来预测是否会崩溃
6.1 数据集
本数据集是在UCI Machine Learning Repository网站下载 一本有540条数据,每条数据有17个特征。现将数据集中80%数据作为训练集,剩下的20%数据作为测试集。
6.2 预测
使用SVM代码如下:
import matplotlib.pyplot as plt
import numpy as np
import random
"""
数据结构,维护所有需要操作的值
Parameters:
dataMatIn - 数据矩阵
classLabels - 数据标签
C - 松弛变量
toler - 容错率
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 = np.shape(dataMatIn)[0]
self.alphas = np.mat(np.zeros((self.m,1)))
self.b = 0
self.eCache = np.mat(np.zeros((self.m,2)))
self.K = np.mat(np.zeros((self.m,self.m)))
for i in range(self.m):
self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)
"""
通过核函数将数据转换更高维的空间
Parameters:
X - 数据矩阵
A - 单个数据的向量
kTup - 包含核函数信息的元组
Returns:
K - 计算的核K
"""
def kernelTrans(X, A, kTup):
m,n = np.shape(X)
K = np.mat(np.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 = np.exp(K/(-1*kTup[1]**2))
else: raise NameError('核函数无法识别')
return K
"""
读取数据
Parameters:
fileName - 文件名
Returns:
dataMat - 数据矩阵
labelMat - 数据标签
"""
def loadDataSet(fileName):
dataMat = []; labelMat = []
fr = open(fileName)
for line in fr.readlines():
lineArr = line.strip().split('\t')
temp = []
for i in range(18):
temp.append(float(lineArr[i]))
dataMat.append(temp)
labelMat.append(float(lineArr[-1]))
return dataMat,labelMat
"""
计算误差
Parameters:
oS - 数据结构
k - 标号为k的数据
Returns:
Ek - 标号为k的数据误差
"""
def calcEk(oS, k):
fXk = float(np.multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
Ek = fXk - float(oS.labelMat[k])
return Ek
"""
函数说明:随机选择alpha_j的索引值
Parameters:
i - alpha_i的索引值
m - alpha参数个数
Returns:
j - alpha_j的索引值
"""
def selectJrand(i, m):
j = i
while (j == i):
j = int(random.uniform(0, m))
return j
"""
内循环启发方式2
Parameters:
i - 标号为i的数据的索引值
oS - 数据结构
Ei - 标号为i的数据误差
Returns:
j, maxK - 标号为j或maxK的数据的索引值
Ej - 标号为j的数据误差
"""
def selectJ(i, oS, Ei):
maxK = -1; maxDeltaE = 0; Ej = 0
oS.eCache[i] = [1,Ei]
validEcacheList = np.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
"""
计算Ek,并更新误差缓存
Parameters:
oS - 数据结构
k - 标号为k的数据的索引值
Returns:
无
"""
def updateEk(oS, k):
Ek = calcEk(oS, k)
oS.eCache[k] = [1,Ek]
"""
修剪alpha_j
Parameters:
aj - alpha_j的值
H - alpha上限
L - alpha下限
Returns:
aj - 修剪后的alpah_j的值
"""
def clipAlpha(aj,H,L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj
"""
优化的SMO算法
Parameters:
i - 标号为i的数据的索引值
oS - 数据结构
Returns:
1 - 有任意一对alpha值发生变化
0 - 没有任意一对alpha值发生变化或变化太小
"""
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)
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]
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("alpha_j变化太小")
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.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
"""
完整的线性SMO算法
Parameters:
dataMatIn - 数据矩阵
classLabels - 数据标签
C - 松弛变量
toler - 容错率
maxIter - 最大迭代次数
kTup - 包含核函数信息的元组
Returns:
oS.b - SMO算法计算的b
oS.alphas - SMO算法计算的alphas
"""
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup = ('lin',0)):
oS = optStruct(np.mat(dataMatIn), np.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("全样本遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter,i,alphaPairsChanged))
iter += 1
else:
nonBoundIs = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
for i in nonBoundIs:
alphaPairsChanged += innerL(i,oS)
print("非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter,i,alphaPairsChanged))
iter += 1
if entireSet:
entireSet = False
elif (alphaPairsChanged == 0):
entireSet = True
print("迭代次数: %d" % iter)
return oS.b,oS.alphas
"""
Parameters:
kTup - 包含核函数信息的元组
Returns:
无
"""
def testDigits(kTup=('rbf', 10)):
dataArr, labelArr = loadDataSet('pop_failures_training.txt')
b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10, kTup)
datMat = np.mat(dataArr); labelMat = np.mat(labelArr).transpose()
svInd = np.nonzero(alphas.A>0)[0]
sVs=datMat[svInd]
labelSV = labelMat[svInd];
print("支持向量个数:%d" % np.shape(sVs)[0])
m,n = np.shape(datMat)
errorCount = 0
for i in range(m):
kernelEval = kernelTrans(sVs,datMat[i,:],kTup)
predict=kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b
if np.sign(predict) != np.sign(labelArr[i]): errorCount += 1
print("训练集错误率: %.2f%%" % (float(errorCount)/m))
dataArr, labelArr = loadDataSet('pop_failures_test.txt')
errorCount = 0
datMat = np.mat(dataArr); labelMat = np.mat(labelArr).transpose()
m,n = np.shape(datMat)
for i in range(m):
kernelEval = kernelTrans(sVs,datMat[i,:],kTup)
predict=kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b
if np.sign(predict) != np.sign(labelArr[i]): errorCount += 1
print("测试集错误率: %.2f%%" % (float(errorCount)/m))
if __name__ == '__main__':
testDigits()
>>>
L==H
全样本遍历:第0次迭代 样本:396, alpha优化次数:0
L==H
全样本遍历:第0次迭代 样本:397, alpha优化次数:0
L==H
全样本遍历:第0次迭代 样本:398, alpha优化次数:0
L==H
全样本遍历:第0次迭代 样本:399, alpha优化次数:0
迭代次数: 1
支持向量个数:0
训练集错误率: 0.91%
测试集错误率: 0.92%
|