感谢：
本次代码修改了原书中的个别bug，增加了原书中没有展示的数据可视化代码
那么，正文就从这里开始啦！（我的代码都是可直接运行的，只要环境正确）
1、svm01_SMO_base.py
'''
简化版的Platt SMO算法
100个测试点，训练需要数秒
数据可视化参考自：
javascript:void(0)
'''
from numpy import *
import matplotlib.pyplot as plt
# 处理数据
def loadData(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的情况下再随机选择一个不一样的alpha
# m 是可选下标长度，i，j皆为下标
def selectJrand(i, m):
j = i
while (j == i):
j = int(random.uniform(0, m))
return j
# 调整大于H，或小于L的alpha值
def clipAlpha(aj, H, L):
if aj > H:
aj = H
if aj < L:
aj = L
return aj
'''
参数：数据集，类别标签，松弛变量C，容错率， 退出前最大的循环次数
C表示不同优化问题的权重，如果C很大，分类器会力图通过分离超平面将样例都正确区分，如果C很小又很难很好的分类，C值需要平衡两者
'''
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
for i in range(m): # 内循环
fXi = float(multiply(alphas, labelMat).T *
(dataMatrix * dataMatrix[i, :].T)) + b # 预测类别，multiply：矩阵对应元素相乘
Ei = fXi - float(labelMat[i]) # 预测误差 = 预测类别 - 实际类别，如果误差很大就要进行优化，如果在0-C之外，则已在边界，没有优化空间了，无需优化
# 如果该数据向量可以被优化，则随机选择另一个数据向量，同时进行优化
if ((labelMat[i] * Ei < -toler) and (alphas[i] < C)) or \
((labelMat[i] * Ei > toler) and (alphas[i] > 0)):
# 随机选择另一个数据向量
j = selectJrand(i, m)
fXj = float(multiply(alphas, labelMat).T *
(dataMatrix * dataMatrix[j, :].T)) + b # 预测类别
Ej = fXj - float(labelMat[j]) # 预测误差 = 预测类别 - 实际类别
# 进行深拷贝，便于之后的新旧值比较
alphaIold = alphas[i].copy()
alphaJold = alphas[j].copy()
# 保证alpha在0-C之间
if labelMat[i] != labelMat[j]:
L = max(0, alphas[j] - alphas[i])
H = min(C, C + alphas[j] - alphas[i])
else:
'''
sum = alphas[i] + alphas[j]
如果 sum > C , 则 H = C，L = sum - C
如果 sum < C , 则 H = sum, L = 0
'''
L = max(0, alphas[i] + alphas[j] - C)
H = min(C, alphas[i] + alphas[j])
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
# eta >= 0的情况比较少，并且优化过程计算复杂，所以此处做了简化处理，直接跳过了
if eta >= 0:
print("eta >= 0")
continue
# 获取alphas[j]的优化值
alphas[j] -= labelMat[j] * (Ei - Ej) / eta
alphas[j] = clipAlpha(alphas[j], H, L)
# 比对原值，看变化是否明显，如果优化并不明显则退出
if abs(alphas[j] - alphaJold) < 1e-5:
print("j not moving enough")
continue
# 对alphas[i]进行和alphas[j]量相同、方向相反的优化
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("iter:", iter, "i:", i, ", pairs changed", alphaPairsChanged)
'''
如果所有向量都没有优化，则增加迭代数目，继续下一次循环
当所有的向量都不再优化，并且达到最大迭代次数才退出，一旦有向量优化，iter都要归零
所以iter存储的是在没有任何alpha改变的情况下遍历数据集的次数
'''
if alphaPairsChanged == 0:
iter += 1
else:
iter = 0
print("iteration number:", iter)
return b, alphas
# 分析数据，数据可视化
def plot(dataMat, labelMat):
dataArr = array(dataMat)
num = shape(dataArr)[0]
xcord1 = []
ycord1 = [] # 标签为1的数据点坐标
xcord0 = []
ycord0 = [] # 标签为0的数据点坐标
for i in range(num):
if int(labelMat[i]) == 1:
xcord1.append(dataArr[i, 0])
ycord1.append(dataArr[i, 1])
else:
xcord0.append(dataArr[i, 0])
ycord0.append(dataArr[i, 1])
fig = plt.figure()
fig.set_size_inches(18.5, 18.5)
ax = fig.add_subplot(111)
ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
ax.scatter(xcord0, ycord0, s=30, c='green', marker='o')
plt.show()
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):
# 大部分的alpha为0，少数几个是值不为0的支持向量
w += multiply(alphas[i] * labelMat[i], X[i, :].T)
return w
# 画图并圈出支持向量的数据点
def plot_support(dataMat, labelMat, support_points, ws, b):
dataArr = array(dataMat)
num = shape(dataArr)[0]
xcord1 = []
ycord1 = [] # 标签为1的数据点坐标
xcord0 = []
ycord0 = [] # 标签为0的数据点坐标
for i in range(num):
if int(labelMat[i]) == 1:
xcord1.append(dataArr[i, 0])
ycord1.append(dataArr[i, 1])
else:
xcord0.append(dataArr[i, 0])
ycord0.append(dataArr[i, 1])
fig = plt.figure()
fig.set_size_inches(18.5, 18.5)
ax = fig.add_subplot(111)
ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
ax.scatter(xcord0, ycord0, s=30, c='green', marker='o')
# 圈出支持向量
for i in range(len(support_points)):
x, y = support_points[i][0]
plt.scatter([x], [y], s=150, c='none', alpha=.5, linewidth=1.5, edgecolor='red')
# 绘制分离超平面
x1 = max(dataMat)[0]
x2 = min(dataMat)[0]
b = float(b)
a1 = float(ws[0][0])
a2 = float(ws[1][0])
y1, y2 = (-b - a1 * x1) / a2, (-b - a1 * x2) / a2
plt.plot([x1, x2], [y1, y2])
ax.grid(True)
plt.show()
if __name__ == "__main__":
dataArr, labelArr = loadData("data/testSet.txt")
# 分析数据，画图
# plot(dataArr, labelArr)
b, alphas = smoSimple(dataArr, labelArr, 0.6, 0.001, 40)
print(b) # [[-3.84798638]]、[[-3.8362795]]
# numpy 特有的 "数组过滤"语法，输出大于零的元素组成的矩阵
ws = calcWs(alphas, dataArr, labelArr)
print(alphas[alphas > 0])
# 支持向量的个数
print(shape(alphas[alphas > 0]))
# 输出为支持向量的数据点
support_points = []
for i in range(100):
if alphas[i] > 0.0:
support_points.append([dataArr[i], labelArr[i]])
print(support_points)
# 画图并圈出支持向量，画出分离超平面
plot_support(dataArr, labelArr, support_points, ws, b)
运行结果：
2、svm02_SMO_full.py
'''
Platt SMO 算法完整版
相比简化版，不同在于alpha的选择方式上
应用了一些能够提速的启发式方法
'''
from numpy import *
import matplotlib.pyplot as plt
from svm01_SMO_base import *
# 采用仅包含一个init方法的类来存储各个重要值，方便输入和调用（调用起来比字典还要方便一点）
class Entity:
def __init__(self, dataMatIn, classLabels, C, toler):
self.X = dataMatIn
self.labelMat = classLabels
self.C = C
self.toler = toler
self.m = shape(dataMatIn)[0]
self.alphas = mat(zeros((self.m, 1)))
self.b = 0
# eCache用于缓存误差，第一列是是否有效的标志位，第二列是实际的E值
self.eCache = mat(zeros((self.m, 2)))
# 计算E值并返回，由于该过程出现频繁，因此单独立为函数
def calcEk(ent, k):
fXk = float(multiply(ent.alphas, ent.labelMat).T *
(ent.X * ent.X[k, :].T)) + ent.b
Ek = fXk - float(ent.labelMat[k])
return Ek
# 用于选择第二个alpha值，返回索引和误差值
def selectJ(i, ent, Ei):
maxK = -1
maxDeltaE = 0
Ej = 0
ent.eCache[i] = [1, Ei]
validEcacheList = nonzero(ent.eCache[:, 0].A)[0] # 获取有效值的行索引
# 从现有的有效误差中寻找第二个alpha
if len(validEcacheList) > 1:
# 如果现有的有效误差至少两个（除了第一个alpha以外至少还有一个有效值），遍历所有有效误差，寻找步长最大化的那个alpha
for k in validEcacheList:
if k == i:
continue
Ek = calcEk(ent, k)
deltaE = abs(Ei - Ek)
if deltaE > maxDeltaE:
maxK = k
maxDeltaE = deltaE
Ej = Ek
return maxK, Ej
else:
# 如果现有有效误差只有1个（即第一个alpha），则说明这是第一次循环，直接在所有数据里面进行随机选择
j = selectJrand(i, ent.m)
Ej = calcEk(ent, j)
return j, Ej
# 计算误差值，并存入缓存
def updateEk(ent, k):
Ek = calcEk(ent, k)
ent.eCache[k] = [1, Ek]
# 内循环
def innerL(i, ent):
Ei = calcEk(ent, i)
if ((ent.labelMat[i] * Ei < -ent.toler) and (ent.alphas[i] < ent.C)) or \
((ent.labelMat[i] * Ei > ent.toler) and (ent.alphas[i] > 0)):
j, Ej = selectJ(i, ent, Ei)
alphaIold = ent.alphas[i].copy()
alphaJold = ent.alphas[j].copy()
if ent.labelMat[i] != ent.labelMat[j]:
L = max(0, ent.alphas[j] - ent.alphas[i])
H = min(ent.C, ent.C + ent.alphas[j] - ent.alphas[i])
else:
L = max(0, ent.alphas[i] + ent.alphas[j] - ent.C)
H = min(ent.C, ent.alphas[i] + ent.alphas[j])
if L == H:
print("L == H")
return 0
# eta是alpha[j]的最优修改量
eta = 2.0 * ent.X[i, :] * ent.X[j, :].T - \
ent.X[i, :] * ent.X[i, :].T - \
ent.X[j, :] * ent.X[j, :].T
# eta >= 0的情况比较少，并且优化过程计算复杂，所以此处做了简化处理，直接跳过了
if eta >= 0:
print("eta >= 0")
return 0
ent.alphas[j] -= ent.labelMat[j] * (Ei - Ej) / eta
ent.alphas[j] = clipAlpha(ent.alphas[j], H, L)
updateEk(ent, j)
# 比对原值，看变化是否明显，如果优化并不明显则退出
if abs(ent.alphas[j] - alphaJold) < 1e-5:
print("j not moving enough")
return 0
ent.alphas[i] += ent.labelMat[j] * ent.labelMat[i] * (alphaJold - ent.alphas[j])
updateEk(ent, i)
b1 = ent.b - Ei - ent.labelMat[i] * (ent.alphas[i] - alphaIold) * ent.X[i, :] * ent.X[i, :].T \
- ent.labelMat[j] * (ent.alphas[j] - alphaJold) * ent.X[i, :] * ent.X[j, :].T
b2 = ent.b - Ej - ent.labelMat[i] * (ent.alphas[i] - alphaIold) * ent.X[i, :] * ent.X[j, :].T \
- ent.labelMat[j] * (ent.alphas[j] - alphaJold) * ent.X[j, :] * ent.X[j, :].T
if (0 < ent.alphas[i]) and (ent.C > ent.alphas[i]):
ent.b = b1
elif (0 < ent.alphas[j]) and (ent.C > ent.alphas[j]):
ent.b = b2
else:
ent.b = (b1 + b2) / 2.0
return 1
else:
return 0
'''
外循环函数
参数：数据集，类别标签，松弛变量C，容错率， 总共最大的循环次数，核函数
C表示不同优化问题的权重，如果C很大，分类器会力图通过分离超平面将样例都正确区分，如果C很小又很难很好的分类，C值需要平衡两者
'''
def smop(dataMatIn, classLabels, C, toler, maxIter):
ent = Entity(mat(dataMatIn), mat(classLabels).transpose(), C, toler)
iter = 0
entireSet = True
alphaPairsChanged = 0
# 循环条件1：iter未超过最大迭代次数 && 循环条件2：上次循环有收获或者遍历方式为遍历所有值
# 因为最终的遍历方式是趋向于遍历非边界值，如果仍在遍历所有值说明还需要更多的训练
while iter < maxIter and ((alphaPairsChanged > 0) or entireSet):
alphaPairsChanged = 0
# 遍历所有的值
if entireSet:
for i in range(ent.m):
alphaPairsChanged += innerL(i, ent)
print("fullSet, iter:", iter, ", i:", i, ", pairs changed:", alphaPairsChanged)
iter += 1
# 遍历非边界值
else:
# 取出非边界值的行索引
nonBoundIs = nonzero((ent.alphas.A > 0) * (ent.alphas.A < ent.C))[0]
for i in nonBoundIs:
alphaPairsChanged += innerL(i, ent)
print("non-bound, iter:", iter, ", i:", i, ", pairs changed: ", alphaPairsChanged)
iter += 1
# 如果上次循环是遍历所有值，那么下次循环改为遍历非边界值
if entireSet:
entireSet = False
# 如果上次循环是遍历非边界值但没有收获，则改为遍历所有值
elif alphaPairsChanged == 0:
entireSet = True
print("iteration number: ", iter)
return ent.b, ent.alphas
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):
# 大部分的alpha为0，少数几个是值不为0的支持向量
w += multiply(alphas[i] * labelMat[i], X[i, :].T)
return w
if __name__ == "__main__":
dataArr, labelArr = loadData("data/testSet.txt")
# 此处取70%的数据作为训练数据
trainArr = dataArr[0:70]
trainLabel = labelArr[0: 70]
# C表示不同优化问题的权重，如果C很大，分类器会力图通过分离超平面将样例都正确区分，如果C很小又很难很好的分类，C值需要平衡两者
b, alphas = smop(trainArr, trainLabel, 0.6, 0.001, 30)
# print(b, alphas) # [[-2.89901748]]
# 计算w值
ws = calcWs(alphas, trainArr, trainLabel)
# print(ws) # [[ 0.7929368 ],[-0.16014993]]
dataMat = mat(dataArr)
# 此处将后面的30%数据作为测试数据（因为数据点本身分布得就很均匀，所以准确率很高）
errorCount = 0.0
testMat = dataMat[70: 100]
testLabel = labelArr[70: 100]
for i in range(len(testMat)):
predict_label = (testMat[i] * mat(ws) + b)[0][0]
# print(predict_label)
if sign(predict_label) != sign(testLabel[i]):
print("预测错误")
errorCount += 1.0
print("错误率为：", (errorCount / len(dataMat)))
运行结果：
3、svm03_radial_basis_func.py
'''
引入径向基函数（一种常见的核函数）来将低维空间的数据映射到高维空间，从而解决低维空间中的非线性问题
note1：smop函数第一行代码需要改变，是书中没有提及的
note2：如果支持向量太少会得到很差的决策边界，如果支持向量太多也相当于每次都利用整个数据集进行分类，这种分类方法接近于k近邻
'''
from svm02_SMO_full import *
'''
核转换函数
kTup是一个元组，第一个参数是核函数类型，第二个参数是核函数可能需要的可选参数
'''
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 Entity:
def __init__(self, dataMatIn, classLabels, C, toler, kTup):
self.X = dataMatIn
self.labelMat = classLabels
self.C = C
self.toler = toler
self.m = shape(dataMatIn)[0]
self.alphas = mat(zeros((self.m, 1)))
self.b = 0
# eCache用于缓存误差，第一列是是否有效的标志位，第二列是实际的E值
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)
# 计算E值并返回，由于该过程出现频繁，因此单独立为函数
def calcEk(ent, k):
fXk = float(multiply(ent.alphas, ent.labelMat).T * ent.K[:, k] + ent.b)
Ek = fXk - float(ent.labelMat[k])
return Ek
# 内循环
# 原来的ent.X[i, :] * ent.X[j, :].T改成了 ent.K[j, j]，只有类似格式进行了变动，其他没变
def innerL(i, ent):
Ei = calcEk(ent, i)
if ((ent.labelMat[i] * Ei < -ent.toler) and (ent.alphas[i] < ent.C)) or \
((ent.labelMat[i] * Ei > ent.toler) and (ent.alphas[i] > 0)):
j, Ej = selectJ(i, ent, Ei)
alphaIold = ent.alphas[i].copy()
alphaJold = ent.alphas[j].copy()
if ent.labelMat[i] != ent.labelMat[j]:
L = max(0, ent.alphas[j] - ent.alphas[i])
H = min(ent.C, ent.C + ent.alphas[j] - ent.alphas[i])
else:
L = max(0, ent.alphas[i] + ent.alphas[j] - ent.C)
H = min(ent.C, ent.alphas[i] + ent.alphas[j])
if L == H:
print("L == H")
return 0
# eta是alpha[j]的最优修改量
eta = 2.0 * ent.K[i, j] - ent.K[i, i] - ent.K[j, j]
# eta >= 0的情况比较少，并且优化过程计算复杂，所以此处做了简化处理，直接跳过了
if eta >= 0:
print("eta >= 0")
return 0
ent.alphas[j] -= ent.labelMat[j] * (Ei - Ej) / eta
ent.alphas[j] = clipAlpha(ent.alphas[j], H, L)
updateEk(ent, j)
# 比对原值，看变化是否明显，如果优化并不明显则退出
if abs(ent.alphas[j] - alphaJold) < 1e-5:
print("j not moving enough")
return 0
ent.alphas[i] += ent.labelMat[j] * ent.labelMat[i] * (alphaJold - ent.alphas[j])
updateEk(ent, i)
b1 = ent.b - Ei - ent.labelMat[i] * (ent.alphas[i] - alphaIold) * ent.K[i, i] \
- ent.labelMat[j] * (ent.alphas[j] - alphaJold) * ent.K[i, j]
b2 = ent.b - Ej - ent.labelMat[i] * (ent.alphas[i] - alphaIold) * ent.K[i, j] \
- ent.labelMat[j] * (ent.alphas[j] - alphaJold) * ent.K[j, j]
if (0 < ent.alphas[i]) and (ent.C > ent.alphas[i]):
ent.b = b1
elif (0 < ent.alphas[j]) and (ent.C > ent.alphas[j]):
ent.b = b2
else:
ent.b = (b1 + b2) / 2.0
return 1
else:
return 0
'''
外循环函数
参数：数据集，类别标签，松弛变量C，容错率， 总共最大的循环次数，核函数
C表示不同优化问题的权重，如果C很大，分类器会力图通过分离超平面将样例都正确区分，如果C很小又很难很好的分类，C值需要平衡两者
'''
def smop(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
ent = Entity(mat(dataMatIn), mat(classLabels).transpose(), C, toler, kTup)
iter = 0
entireSet = True
alphaPairsChanged = 0
# 循环条件1：iter未超过最大迭代次数 && 循环条件2：上次循环有收获或者遍历方式为遍历所有值
# 因为最终的遍历方式是趋向于遍历非边界值，如果仍在遍历所有值说明还需要更多的训练
while iter < maxIter and ((alphaPairsChanged > 0) or entireSet):
alphaPairsChanged = 0
# 遍历所有的值
if entireSet:
for i in range(ent.m):
alphaPairsChanged += innerL(i, ent)
print("fullSet, iter:", iter, ", i:", i, ", pairs changed:", alphaPairsChanged)
iter += 1
# 遍历非边界值
else:
# 取出非边界值的行索引
nonBoundIs = nonzero((ent.alphas.A > 0) * (ent.alphas.A < ent.C))[0]
for i in nonBoundIs:
alphaPairsChanged += innerL(i, ent)
print("non-bound, iter:", iter, ", i:", i, ", pairs changed: ", alphaPairsChanged)
iter += 1
# 如果上次循环是遍历所有值，那么下次循环改为遍历非边界值
if entireSet:
entireSet = False
# 如果上次循环是遍历非边界值但没有收获，则改为遍历所有值
elif alphaPairsChanged == 0:
entireSet = True
print("iteration number: ", iter)
return ent.b, ent.alphas
def testRbf(k1=1.3):
dataArr, labelArr = loadData('data/testSetRBF.txt')
# 参数：数据集，类别标签，松弛变量C，容错率， 总共最大的循环次数，核函数
b, alphas = smop(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1))
dataMat = mat(dataArr)
labelMat = mat(labelArr).transpose()
svInd = nonzero(alphas.A > 0)[0] # 支持向量的行索引
sVs = dataMat[svInd] # 支持向量的数据点
labelSV = labelMat[svInd] # 支持向量的类别标签
print("There are", shape(sVs)[0], "support vectors")
m, n = shape(dataMat)
errorCount = 0
for i in range(m):
kernelEval = kernelTrans(sVs, dataMat[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", (float(errorCount) / m))
dataArr, labelArr = loadData('data/testSetRBF2.txt')
dataMat = mat(dataArr)
m, n = shape(dataMat)
errorCount = 0
for i in range(m):
kernelEval = kernelTrans(sVs, dataMat[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", (float(errorCount) / m))
if __name__ == "__main__":
testRbf()
运行结果：
4、svm04_handwriting_recognition.py
'''
基于SVM的数字识别
1、收集数据
2、准备数据：基于二值图像构造向量
3、分析数据：对图像向量进行目测
4、训练算法：采用两种不同的核函数，并对径向基核函数采用不同的设置来运行SMO算法
5、测试算法：测试，计算错误率
6、使用算法：要求一些图像处理的知识
note: 书中最大迭代次数设置为10000，实际运行要很久，此处改小为10。
'''
from svm03_radial_basis_func import *
# 将数字图像信息转换为矩阵
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])
# 将所有数据排在一行存储到 returnVect 矩阵
return returnVect
'''
将手写识别问题转换为二分类问题，9类别为-1，非9类别为1
'''
def loadImage(dirName):
from os import listdir
hwLabels = []
# 获取所有文件
trainingFileList = listdir(dirName)
m = len(trainingFileList)
trainingMat = zeros((m, 1024))
for i in range(m):
fileFullName = trainingFileList[i]
fileName = fileFullName.split('.')[0]
# 文件名的前缀表示类别标签
classNumStr = int(fileName.split('_')[0])
if classNumStr == 9:
hwLabels.append(-1)
else:
hwLabels.append(1)
trainingMat[i, :] = img2vector(dirName + '/' + fileFullName)
return trainingMat, hwLabels
def testDigits(kTup=('rbf', 10)):
dataArr, labelArr = loadImage('handwriting_image_data2/trainingDigits')
b, alphas = smop(dataArr, labelArr, 200, 0.0001, 10, kTup)
dataMat = mat(dataArr)
labelMat = mat(labelArr).transpose()
svInd = nonzero(alphas.A > 0)[0]
sVs = dataMat[svInd]
labelSV = labelMat[svInd]
print("There are", shape(sVs)[0], "support vectors")
m, n = shape(dataMat)
errorCount = 0
for i in range(m):
kernelEval = kernelTrans(sVs, dataMat[i, :], kTup)
predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
if sign(predict) != sign(labelArr[i]):
errorCount += 1
print("the training error rate is", (float(errorCount) / m))
dataArr, labelArr = loadImage('handwriting_image_data2/testDigits')
dataMat = mat(dataArr)
m, n = shape(dataMat)
errorCount = 0
for i in range(m):
kernelEval = kernelTrans(sVs, dataMat[i, :], kTup)
predict = kernelEval.T * multiply(labelSV, alphas[svInd]) + b
if sign(predict) != sign(labelArr[i]):
errorCount += 1
print("the test error rate is", (float(errorCount) / m))
if __name__ == "__main__":
testDigits() # training error rate：0.08706467661691543，test error rate：0.05913978494623656
运行结果：（运行稍久）
