Machine Learning in Action 读书笔记---第6章支持向量机

最新推荐文章于 2024-09-21 18:16:51 发布

dongyuyuu

最新推荐文章于 2024-09-21 18:16:51 发布

阅读量382

点赞数

分类专栏：【机器学习】文章标签：支持向量机 svm python

本文链接：https://blog.csdn.net/id__39/article/details/119455120

版权

【机器学习】专栏收录该内容

14 篇文章 0 订阅

订阅专栏

Machine Learning in Action 读书笔记

第6章支持向量机

文章目录

Machine Learning in Action 读书笔记
一、本章专业词汇
二、支持向量机
- 1.分类器求解的优化问题
- 2.SVM的一般流程
三、SMO高效优化算法
四、在复杂数据上应用核函数
- 1.利用核函数将数据映射到高维空间
- 2.径向基核函数
五、示例：手写识别问题
六、完整代码

一、本章专业词汇

支持向量机（Support Vector Machines，SVM）：是一个现成的分类器，“现成”是指分类器不加修改即可直接使用。SVM能够对训练集之外的数据点做出很好的分类决策。
序列最小优化（Sequential Minimal Optimization，SMO）算法：实现SVM的一个算法。
核函数（kernel）：使用核函数可以将SVM扩展到更多数据集上
线性可分（linearly separable）数据：两组数据之间已经分隔得足够开，很容易就可以画出一条直线将两组数据点分开。
分隔超平面（separating hyperplane）：将线性可分的数据集分隔开来的直线称为分隔超平面。
超平面（hyperplane）：分类的决策边界。如果数据点离决策边界越远，那么其最后的预测结果也就越可信。
间隔（margin）：点到分割面的距离被称为间隔。我们希望间隔尽可能的大。
支持向量（support vector）：就是离分隔超平面最近的那些点。我们需要试着最大化支持向量到分隔面的距离。
松弛变量（slack variable）：允许有些数据点可以处于分割面的错误一侧。
数组过滤（array filtering）：数组过滤只对Numpy类型有用，不适用于python中的正则表（regular list）。例如：print(alphas[alphas > 0])
核函数（kernel）：可以将数据转换成易于分类器理解的形式
径向基函数（radial basis function）：核函数中的一种
内积（inner product）（也称点积）：指两个向量相乘，之后得到单个标量或者数值的运算。
核技巧（kernel trick）：将内积替换成核函数的方式被称为核技巧或者核“变电”（kernel substation）

二、支持向量机

优点：泛化错误率低，计算开销不大，结果易解释。
缺点：对参数调节和核函数的选择敏感，原始分类器不加修改仅适用于处理二类问题。
适用数据类型：数值型和标称型数据

1.分类器求解的优化问题

分类，使用类似于海维赛德阶跃函数（即单位阶跃函数）的函数对：
$w^Tx+b$
作用得到：
$f(u) ,u=w^Tx+b$
其中，当u < 0时，f(u)输出-1，反之输出+1。

当计算数据点到分隔面的距离并确定分隔面的放置位置时，计算间隔通过：
$label*(w^Tx+b)$
来计算。这时就能体现出-1和+1类的好处了，如果数据点处于正上方（即+1类），并且离分隔超平面很远的位置时，u会是一个很大的正数，同时 label * u 也会是一个很大的正数，而如果数据点处于负方向（-1类）并且离分隔超平面很远的位置时，此时由于类别标签为-1，则 label * u 仍然是一个很大的正数。

现在的目标是找出分类器定义中的 w 和 b 。为此我们必须找到具有最小间隔的数据点，也就是找到支持向量。一旦找到具有最小间隔的数据点，我们就需要对该间隔最大化，也可以写做：
$arg\ \mathop{max}\limits_{w,b}\ \{\mathop{\ min}\limits_{n}(label·(w^Tx+b))·\frac{1}{||w||}\}$
由于对乘积进行优化是一件很讨厌的事情，因此我们要做的是固定其中一个因子，而最大化其他因子。如果令u都为1，那么就可以通过求：
$\frac{1}{||w||}$
的最大值来得到最终解。但是并非所有数据点的u都等于1，只有那些离超平面最近的点得到的值才是1.

在上述优化问题中，给定了一些约束条件然后求最优值(固定一个因子)，因此该问题是一个带约束条件的优化问题。这里的约束条件就是：
$label*(w^Tx+b)>=1.0$
对于这类优化问题，有一个非常著名的求解方法，拉格朗日乘子法。

由于这里的约束条件都是基于数据点的，因此我们可以将超平面写成数据点的形式：
$\mathop{max}\limits_{α}\ [\ \sum_{i=1}^mα-\frac{1}{2}\sum_{i,j=1}^mlabel^{(i)}·label^{(j)}·α_i·α_j<x^{(i)},x^{(j)}>\ ]$
尖括号表示两向量的内积。其约束条件变为：
$α>=0,和\sum_{i-1}^mα_i·label^{(i)}=0$
通过松弛变量允许有些数据点可以处于分隔面的错误一侧，优化目标仍然保持不变，但是此时的约束条件变为：
$C>=α>=0,和\sum_{i-1}^mα_i·label^{(i)}=0$
这里的常数C用于控制“最大化间隔”和“保证大部分点的函数间隔小于1.0”这两个目标权重。常数C是一个参数，我们可以通过调节该参数得到不同的结果。一旦求出所有α，那么分隔超平面就可以通过这些alpha来表达。
SVM中的主要工作：就是求解这些alpha。

2.SVM的一般流程

收集数据：可以使用任意方法
准备数据：需要数值型数据
分析数据：有助于可视化分隔超平面
训练算法：SVM的大部分时间都源于训练，该过程主要实现两个参数的调优
测试算法：十分简单的计算过程就可以实现
使用算法：几乎所有分类问题都可以使用SVM，SVM本身是一个二类分类器，对多类问题应用SVM需要对代码做一些修改

三、SMO高效优化算法

1996年，John Platt发布了一个称为SMO的强大算法，用于训练SVM。Paltt的SMO算法是将大优化问题分解为多个小优化问题来求解的。

SMO算法的目标是求出一系列alpha和b，一旦求出这些alpha，就很容易计算出权重向量w并得到分隔超平面。

SMO算法的工作原理是：每次循环中选择两个alpha进行优化处理，一旦找到一对合适的alpha，那么就增大其中一个同时减小另一个。“合适”是指：满足两个条件：

这两个alpha必须要在间隔边界之外
这两个alpha还没有进行过区间化处理或者不在边界上。

1.应用简化版SMO算法处理小规模数据集

Platt SMO算法中的外循环确定要优化的最佳alpha对，而简化版的会跳过这一部分，过程如下：

首先在数据集上遍历每一个alpha
然后再剩下的alpha集合中随机选择另一个alpha，从而构建alpha对

我们要同时改变两个alpha，因为我们要满足约束条件：
$\sumα_i·label^{(i)}=0$
如果只改变一个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): # 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

2.简化版SMO算法

伪代码如下：

创建一个alpha向量并将其初始化为0向量
当迭代次数小于最大迭代次数时（外循环）：
	对数据集中的每个数据向量（内循环）：
		如果该数据向量可以被优化：
			随机选择另外一个数据向量
			同时优化这两个向量
			如果两个向量都不能被优化，退出内循环
	如果所有向量都没被优化，增加迭代数目，继续下一次循环

简化版SMO算法：

'''简化版序列最小优化算法（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)))   # 将所有alpha初始化为0
    iter = 0   # 该变量存储 在没有任何alpha改变的情况下遍历数据集的次数

    # 只有在所有数据集上遍历maxIter次，且不发生任何alpha修改之后，程序才会停止并推出while循环
    while(iter < maxIter):
        alphaPairsChanged = 0   # 每次循环都要设置为0，用于记录alpha是否已经进行优化
        for i in range(m):
            # fXi 为我们预测的类别
            fXi = float(multiply(alphas, labelMat).T * \
                        (dataMatrix * dataMatrix[i,:].T)) + b
            Ei = fXi - float(labelMat[i])  #预测结果和真实结果的误差

            # 如果误差很大，对该数据实例对应的alpha值进行优化（对正、负间隔都进行判断，同时保证alpha的值不等于0或C）
            #由于后面alpha小于0或大于C时将被调整为0或C，所以一旦在该if语句中它们等于这两个值，它们就已经在“边界”上了，不再能够减小或增大，也就不值得再对它们进行优化了
            if((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or \
                    ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):# 约束条件 C >=α >=0
                j = selectJrand(i, m)  # 利用辅助函数随机选择第二个alpha值
                fXj = float(multiply(alphas, labelMat).T * \
                            (dataMatrix * dataMatrix[j, :].T)) + b
                Ej = fXj - float(labelMat[j])  # 同样对第二个alpha值计算误差
                alphaIold = alphas[i].copy() # 复制结果，用于下面比较
                alphaJold = alphas[j].copy()
                if(labelMat[i] != labelMat[j]):  # 保证alpha在0到C之间，将alpha[j]调到0和C之间
                    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')# 输出到这里，意味着alpha对没有改变成功 1
                    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') # 输出到这里，意味着alpha对没有改变成功 2
                    continue
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta
                alphas[j] = clipAlpha(alphas[j], H, L)  # 调用了第二个辅助函数，alphas[j]数值太大时，对其进行调整，调整到L和H之间，包括边界
                if (abs(alphas[j] - alphaJold) < 0.00001):  # 检查alpha[j]是否有轻微的改变
                    print('j not moving enough')  # 输出到这里，意味着alpha对没有改变成功 3
                    continue
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])  # alpha[i]和alpha[j]一样，同样进行改变，改变的大小相同，但方向相反（一个增，另一个减）

                # 对alpha[i]和alpha[j]优化后，给这两个alpha值设置一个常数项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  # 程序执行到for循环的最后一行都没有执行continue语句，就说明已经成功的改变了一对alpha，同时增加alphaPairsChanged的值
                print('iter: %d i: %d, pairs changed：%d' % (iter, i, alphaPairsChanged))
        if alphaPairsChanged == 0: # 检查alpha的值是否做了更新
            iter += 1
        else:
            iter = 0  # 如果有更新则将iter设为0，后继续运行程序
        print('iteration number: %d' % iter)
    return b, alphas

根据获得的alpha值，可以：

得到支持向量的个数：

print(alphas[alphas > 0]) #[[0.14742212 0.17251816 0.04511926 0.36505954]]

输出支持向量的个数

print(shape(alphas[alphas > 0]))  # (1, 4)

输出哪些是支持向量

    for i in range(100):
        if alphas[i]>0.0:
            print(dataArr[i], labelArr[i])
    '''
                [4.658191, 3.507396] -1.0
                [3.457096, -0.082216] -1.0
                [5.286862, -2.358286] 1.0
                [6.080573, 0.418886] 1.03
    '''

但是使用简化版的SMO算法获得支持向量，收敛时间较长，通过构建完整的SMO算法可以加快其运行速度。

3.完整版SMO算法

完整版SMO算法在优化过程中，唯一的不同就是选择alpha的方式。完整版Platt SMO算法应用了一些能够提速的启发方法。
Platt SMO算法是通过一个外循环来选择第一个alpha值的，并且选择过程会在两种方式之间交替：

一种方式是在所有数据集上进行单遍扫描
另一种方式则是在非边界alpha中实现单遍扫描，非边界alpha指的就是那些不等于边界0或C的alpha值。

而实现非边界alpha值的扫描时，首先需要建立这些alpha值的列表，然后再对这个表进行遍历，同时，该步骤会跳过那些已知的不会改变的alpha值。

在选择了第一个alpha值后，算法会通过内循环选择第二个alpha值，在优化过程中会通过最大化步长的方式来获得第二个alpha值。
在简化版SMO算法中，我们会在选择j之后计算错误率Ej，但在完整版SMO算法中，我们会建立一个全局的缓存用于保存误差值，并从中选择使得步长或者说Ei-Ej最大的alpha值。

完整版Platt SMO算法的支持函数：

'''完整版 Platt SMO 算法'''
# 完整版SMO 的支持函数：包含一个用于清理代码的数据结构和3个用于对E进行缓存的辅助函数
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)))  # 误差缓存 m*2矩阵，eCache第一列给出的是eCache是否有效的标志位（有效意味着已经计算好了），第二列给出的是实际的E值

def calcEk(oS, k):  # 计算误差值，并返回，由于使用频繁，所以单独拎出来作为一个函数使用
    fXk = float(multiply(oS.alphas, oS.labelMat).T * \
                (oS.X * oS.X[k, :].T)) + oS.b
    Ek = fXk - float(oS.labelMat[k])
    return Ek

def selectJ(i, oS, Ei):  # 内循环中的启发式方法，用于选择第二个alpha，选择合适的alpha值以保证每次优化中采用最大步长
    maxK = -1; maxDeltaE = 0; Ej = 0
    oS.eCache[i] = [1, Ei]
    validEcacheList = 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: # 选择具有最大步长的j
                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]

# 完整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 = 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):  # 检查alpha[j]是否有轻微的改变
            print('j not moving enough')
            return 0
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (\
                        alphaJold - oS.alphas[j])  # alpha[i]和alpha[j]一样，同样进行改变，改变的大小相同，但方向相反（一个增，另一个减）

        updateEk(oS, i)  # 完整版增加内容，更新误差缓存

        # 对alpha[i]和alpha[j]优化后，给这两个alpha值设置一个常数项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
        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

innerL()函数与smoSimple()函数一摸一样，但是有一些改变：

使用了自己的数据结构，该结构在参数oS中传递
选择第二个alpha的值使用selectJ()函数而不是selectJrand()函数
在alpha值改变时，更新Ecache

接下来将上面的过程打包在一起，这就是选择第一个alpha值的外循环

# 将上述代码打包在一起，选择第一个alpha值的外循环
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
    oS = optStruct(mat(dataMatIn), 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)  # 通过innerL选择第二个alpha值
                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

上面代码中，常数C一方面要保障所有样例的间隔不小于1.0，另一方面又要使得分类间隔要尽可能大，并且要在这两方面平衡。

上面我们用了大量的时间来计算alpha值，但是如何利用alpha值进行分类呢？

首先，必须基于alpha值得到超平面，这也包括了w的计算

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

分类器分类结果测试：

 # 完整SMO算法测试
    b, alphas = smoP(dataArr, labelArr, 0.6, 0.001, 40)
    # print(b)  #[[-2.89901748]]
    # print(alphas[alphas > 0]) #[[0.06961952 0.0169055  0.0169055  0.0272699  0.04522972 0.0272699 0.0243898  0.06140181 0.06140181]]

    ws = calcWs(alphas, dataArr,labelArr)
    # print(ws)

    datMat = mat(dataArr)
    if (datMat[0] * mat(ws) + b) > 0:
        print('属于1类')
    else:
        print('属于-1类')
    # 确认分类是否正确：
    print(labelArr[0])
    '''上面的输出结果：属于-1类
                    -1.0
    '''

四、在复杂数据上应用核函数

1.利用核函数将数据映射到高维空间

从一个特征空间到另一个特征空间的映射，这种映射会将低维特征空间映射到高维空间，在高维空间中解决线性问题，也就等价于在低维空间中解决非线性问题。这种映射时通过核函数来实现的。可以将核函数想象成一个包装器或者接口，它能将数据从某个很难处理的形式转换成另一个较容易处理的形式。
SVM优化中一个特别好的地方就是，所有运算都可以写成内积（inner product）也称点积的形式。

2.径向基核函数

径向基核函数是一个常用的度量两个向量距离的核函数。径向基核函数的高斯版本，具体公式如下：
$k(x,y)=exp(\frac{-||x-y||^2}{2σ^2})$
σ是用户定义的用于确定到达率（reach）或者说函数值跌落到0的速度参数。
核抓换函数：

'''核转换函数'''
def kernelTrans(X, A, kTup):  # 函数调用：kernelEval = kernelTrans(sVs, datMat[i, :], ('rbf', k1)) sVs为支持向量
    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 # 返回k向量

'''使用核函数需要对innerL()及calcEk()函数进行修改，修改后的函数如下：'''
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]   # 此处进行了修改 1
        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):  # 检查alpha[j]是否有轻微的改变
            print('j not moving enough')
            return 0
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (\
                        alphaJold - oS.alphas[j])  # alpha[i]和alpha[j]一样，同样进行改变，改变的大小相同，但方向相反（一个增，另一个减）

        updateEk(oS, i)

        # 对alpha[i]和alpha[j]优化后，给这两个alpha值设置一个常数项b  # 修改处 2
        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 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 testRbf(k1=1.3):  # 输入参数是高斯径向基函数中的一个用户定义变量
    dataArr, labelArr = loadDataSet('testSet.txt')
    b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1))
    datMat = mat(dataArr)
    labelMat = mat(labelArr).transpose()
    svInd = nonzero(alphas.A>0)[0] # alphas.A>0返回True or False; nonzero:返回一个二维数组，第一个数组从行维度来描述非零元素的索引值，第二个数组从列的维度来描述
    # print('svInd:', svInd)  #svInd: [17 39 55 57 64 70 74 94] alpha为100*1维矩阵，列只有0列，故返回非零alpha的行数
    sVs = datMat[svInd]
    # print('sVs:', sVs)  #sVs: [[ 4.658191  3.507396].....  # 返回非零alpha所在行的元素特征值
    labelSV = labelMat[svInd]
    # print('labelSV:', labelSV)   # labelSV: [[-1.]...... # 返回支持向量的类别
    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))  # sVs：支持向量的特征
        # print('kernelEval:', kernelEval)
        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))

上述代码中：首先，程序从文件中读入数据集，然后再该数据集上运行Platt SMO算法，其中核函数的类型为 ’ rbf '。
最后输出测试错误率。训练错误率、支持向量个数，观察它们随k1的变化而变化的情况。
通过修改发现，如果降低σ的值，训练错误率就会降低，但是测试错误率却会上升。
SVM的优点在于它能对数据进行高效分类，如果支持向量太少，就可能会得到一个很差的决策边界；如果支持向量太多，也就相当于每次都利用整个数据集进行分类，这种类方法称为k近邻。

五、示例：手写识别问题

代码见完整代码。

六、完整代码

import random
from numpy import *

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): # 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

'''简化版序列最小优化算法（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   # 该变量存储 在没有任何alpha改变的情况下遍历数据集的次数

    # 只有在所有数据集上遍历maxIter次，且不发生任何alpha修改之后，程序才会停止并推出while循环
    while(iter < maxIter):
        alphaPairsChanged = 0   # 每次循环都要设置为0，用于记录alpha是否已经进行优化
        for i in range(m):
            # fXi 为我们预测的类别
            fXi = float(multiply(alphas, labelMat).T * \
                        (dataMatrix * dataMatrix[i,:].T)) + b
            Ei = fXi - float(labelMat[i])  #预测结果和真实结果的误差

            # 如果误差很大，对该数据实例对应的alpha值进行优化（对正、负间隔都进行判断，同时保证alpha的值不等于0或C）
            #由于后面alpha小于0或大于C时将被调整为0或C，所以一旦在该if语句中它们等于这两个值，它们就已经在“边界”上了，不再能够减小或增大，也就不值得再对它们进行优化了
            if((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or \
                    ((labelMat[i]*Ei > toler) and (alphas[i] > 0)): # 约束条件 C >=α >=0
                j = selectJrand(i, m)  # 利用辅助函数随机选择第二个alpha值
                fXj = float(multiply(alphas, labelMat).T * \
                            (dataMatrix * dataMatrix[j, :].T)) + b
                Ej = fXj - float(labelMat[j])  # 同样对第二个alpha值计算误差
                alphaIold = alphas[i].copy() # 复制结果，用于下面比较
                alphaJold = alphas[j].copy()
                if(labelMat[i] != labelMat[j]):  # 保证alpha在0到C之间，将alpha[j]调到0和C之间
                    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') # 输出到这里，意味着alpha对没有改变成功
                    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] -= labelMat[j]*(Ei - Ej)/eta
                alphas[j] = clipAlpha(alphas[j], H, L)
                if (abs(alphas[j] - alphaJold) < 0.00001):  # 检查alpha[j]是否有轻微的改变
                    print('j not moving enough')
                    continue
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])  # alpha[i]和alpha[j]一样，同样进行改变，改变的大小相同，但方向相反（一个增，另一个减）

                # 对alpha[i]和alpha[j]优化后，给这两个alpha值设置一个常数项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  # 程序执行到for循环的最后一行都没有执行continue语句，就说明已经成功的改变了一对alpha，同时增加alphaPairsChanged的值
                print('iter: %d i: %d, pairs changed：%d' % (iter, i, alphaPairsChanged))
        if alphaPairsChanged == 0: # 检查alpha的值是否做了更新
            iter += 1
        else:
            iter = 0  # 如果有更新则将iter设为0，后继续运行程序
        print('iteration number: %d' % iter)
    return b, alphas

'''完整版 Platt SMO 算法'''
# 完整版SMO 的支持函数：包含一个用于清理代码的数据结构和3个用于对E进行缓存的辅助函数
class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler, kTup): #使用核函数新添加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):
            # X为描述所用核函数类型的一个字符串，其他两个参数都为核函数可能需要的可选参数
            self.K[:, i] = kernelTrans(self.X, self.X[i, :], kTup)  # 全局的矩阵k值只需要计算一次，当想要使用核函数时，就可以对它进行调用，省去了很多冗余的计算开销


# def calcEk(oS, k):
#     fXk = float(multiply(oS.alphas, oS.labelMat).T * \
#                 (oS.X * oS.X[k, :].T)) + oS.b
#     Ek = fXk - float(oS.labelMat[k])
#     return Ek

def selectJ(i, oS, Ei):
    maxK = -1; maxDeltaE = 0; Ej = 0
    oS.eCache[i] = [1, Ei]
    validEcacheList = 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

def updateEk(oS, k):
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1, Ek]

# 完整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)
#         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):  # 检查alpha[j]是否有轻微的改变
#             print('j not moving enough')
#             return 0
#         oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (\
#                         alphaJold - oS.alphas[j])  # alpha[i]和alpha[j]一样，同样进行改变，改变的大小相同，但方向相反（一个增，另一个减）
#
#         updateEk(oS, i)
#
#         # 对alpha[i]和alpha[j]优化后，给这两个alpha值设置一个常数项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
#         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

# 将上述代码打包在一起，选择第一个alpha值的外循环
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
    # oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler)
    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

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 kernelTrans(X, A, kTup):  # 函数调用：kernelEval = kernelTrans(sVs, datMat[i, :], ('rbf', k1)) sVs为支持向量
    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 # 返回k向量
'''使用核函数需要对innerL()及calcEk()函数进行修改，修改后的函数如下：'''
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]   # 此处进行了修改 1
        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):  # 检查alpha[j]是否有轻微的改变
            print('j not moving enough')
            return 0
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (\
                        alphaJold - oS.alphas[j])  # alpha[i]和alpha[j]一样，同样进行改变，改变的大小相同，但方向相反（一个增，另一个减）

        updateEk(oS, i)

        # 对alpha[i]和alpha[j]优化后，给这两个alpha值设置一个常数项b  # 修改处 2
        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 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 testRbf(k1=1.3):
    dataArr, labelArr = loadDataSet('testSet.txt')
    b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1))
    datMat = mat(dataArr)
    labelMat = mat(labelArr).transpose()
    svInd = nonzero(alphas.A>0)[0] # alphas.A>0返回True or False; nonzero:返回一个二维数组，第一个数组从行维度来描述非零元素的索引值，第二个数组从列的维度来描述
    # print('svInd:', svInd)  #svInd: [17 39 55 57 64 70 74 94] alpha为100*1维矩阵，列只有0列，故返回非零alpha的行数
    sVs = datMat[svInd]
    # print('sVs:', sVs)  #sVs: [[ 4.658191  3.507396].....  # 返回非零alpha所在行的元素特征值
    labelSV = labelMat[svInd]
    # print('labelSV:', labelSV)   # labelSV: [[-1.]...... # 返回支持向量的类别
    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))  # sVs：支持向量的特征
        # print('kernelEval:', kernelEval)
        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的手写数字识别'''
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)           #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))

if __name__ == '__main__':
    dataArr, labelArr = loadDataSet('testSet.txt')
    # print(dataArr)  # [[3.542485, 1.977398], [3.018896, 2.556416],
    # print(labelArr)

    #简化版SMO算法 测试
    # b, alphas = smoSimple(dataArr, labelArr, 0.6, 0.001, 40)

    # print(b)  #[[-3.78586587]]
    # print(alphas[alphas > 0]) #[[0.14742212 0.17251816 0.04511926 0.36505954]]

    # 输出 支持向量的个数
    # print(shape(alphas[alphas > 0]))  # (1, 4)

    # 输出是 支持向量的点
    # for i in range(100):
    #     if alphas[i]>0.0:
    #         print(dataArr[i], labelArr[i])
    '''
                [4.658191, 3.507396] -1.0
                [3.457096, -0.082216] -1.0
                [5.286862, -2.358286] 1.0
                [6.080573, 0.418886] 1.03
    '''
    # 完整SMO算法测试
    # b, alphas = smoP(dataArr, labelArr, 0.6, 0.001, 40)
    # print(b)  #[[-2.89901748]]
    # print(alphas[alphas > 0]) #[[0.06961952 0.0169055  0.0169055  0.0272699  0.04522972 0.0272699 0.0243898  0.06140181 0.06140181]]

    # ws = calcWs(alphas, dataArr,labelArr)
    # print(ws)

    # datMat = mat(dataArr)
    # if (datMat[0] * mat(ws) + b) > 0:
    #     print('属于1类')
    # else:
    #     print('属于-1类')
    # # 确认分类是否正确：
    # print(labelArr[0])
    '''上面的输出结果：属于-1类
                    -1.0
    '''
    # testRbf(2.0)

    # 手写数字测试
    testDigits(('rbf', 20))