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

野性酷女 2022-09-03 08:27 210阅读 0赞

# Machine Learning in Action 读书笔记 #

第6章 支持向量机

--------------------

### 文章目录 ###

*  Machine Learning in Action 读书笔记
 *  一、本章专业词汇
 *  二、支持向量机
 *   *  1.分类器求解的优化问题
     *  2.SVM的一般流程
 *  三、SMO高效优化算法
 *   *  1.应用简化版SMO算法处理小规模数据集
     *  2.简化版SMO算法
     *  3.完整版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 T x + b w^Tx+b wTx\+b  
作用得到：  
 f ( u ) , u = w T x + b f(u) ,u=w^Tx+b f(u),u=wTx\+b  
其中，当u < 0时，f(u)输出-1，反之输出+1。  
  
  
当计算数据点到分隔面的距离并确定分隔面的放置位置时，**计算间隔**通过：  
 l a b e l ∗ ( w T x + b ) label\*(w^Tx+b) label∗(wTx\+b)  
来计算。这时就能体现出-1和+1类的好处了，如果数据点处于正上方（即+1类），并且离分隔超平面很远的位置时，u会是一个很大的正数，同时 label \* u 也会是一个很大的正数，而如果数据点处于负方向（-1类）并且离分隔超平面很远的位置时，此时由于类别标签为-1，则 label \* u 仍然是一个很大的正数。  
  
  
现在的目标是找出分类器定义中的 **w** 和 **b** 。为此我们必须找到具有最小间隔的数据点，也就是**找到支持向量**。一旦找到具有最小间隔的数据点，我们就需要对该**间隔最大化**，也可以写做：  
 a r g m a x w , b \{ m i n n ( l a b e l ⋅ ( w T x + b ) ) ⋅ 1 ∣ ∣ w ∣ ∣ \} arg\\ \\mathop\{max\}\\limits\_\{w,b\}\\ \\\{\\mathop\{\\ min\}\\limits\_\{n\}(label·(w^Tx+b))·\\frac\{1\}\{||w||\}\\\} argw,bmax\{ nmin(label⋅(wTx\+b))⋅∣∣w∣∣1\}  
由于对乘积进行优化是一件很讨厌的事情，因此我们要做的是固定其中一个因子，而最大化其他因子。如果令u都为1，那么就可以通过求：  
 1 ∣ ∣ w ∣ ∣ \\frac\{1\}\{||w||\} ∣∣w∣∣1  
的最大值来得到最终解。但是并非所有数据点的u都等于1，只有那些离超平面最近的点得到的值才是1.  
  
  
在上述优化问题中，给定了一些约束条件然后求最优值(固定一个因子)，因此该问题是一个带约束条件的优化问题。这里的约束条件就是：  
 l a b e l ∗ ( w T x + b ) > = 1.0 label\*(w^Tx+b)>=1.0 label∗(wTx\+b)>=1.0  
对于这类优化问题，有一个非常著名的求解方法，拉格朗日乘子法。  
  
  
由于这里的约束条件都是基于数据点的，因此我们可以将超平面写成数据点的形式：  
 m a x α \[ ∑ 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 ) > \] \\mathop\{max\}\\limits\_\{α\}\\ \[\\ \\sum\_\{i=1\}^mα-\\frac\{1\}\{2\}\\sum\_\{i,j=1\}^mlabel^\{(i)\}·label^\{(j)\}·α\_i·α\_j<x^\{(i)\},x^\{(j)\}>\\ \] αmax\[i=1∑mα−21i,j=1∑mlabel(i)⋅label(j)⋅αi⋅αj<x(i),x(j)>\]  
尖括号表示两向量的内积。其约束条件变为：  
 α > = 0 , 和 ∑ i − 1 m α i ⋅ l a b e l ( i ) = 0 α>=0,和\\sum\_\{i-1\}^mα\_i·label^\{(i)\}=0 α>=0,和i−1∑mαi⋅label(i)=0  
通过松弛变量允许有些数据点可以处于分隔面的错误一侧，优化目标仍然保持不变，但是此时的约束条件变为：  
 C > = α > = 0 , 和 ∑ i − 1 m α i ⋅ l a b e l ( i ) = 0 C>=α>=0,和\\sum\_\{i-1\}^mα\_i·label^\{(i)\}=0 C>=α>=0,和i−1∑mαi⋅label(i)=0  
这里的**常数C**用于控制“最大化间隔”和“保证大部分点的函数间隔小于1.0”这两个目标权重。常数C是一个参数，我们可以通过调节该参数得到不同的结果。一旦求出所有α，那么分隔超平面就可以通过这些alpha来表达。  
SVM中的主要工作：就是求解这些alpha。

## 2.SVM的一般流程 ##

1.  收集数据：可以使用任意方法
2.  准备数据：需要**数值型数据**
3.  分析数据：有助于可视化分隔超平面
4.  训练算法：SVM的大部分时间都源于训练，**该过程主要实现两个参数的调优**
5.  测试算法：十分简单的计算过程就可以实现
6.  使用算法：几乎所有分类问题都可以使用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对，而简化版的会跳过这一部分，过程如下：

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

我们要同时改变两个alpha，因为我们要满足约束条件：  
 ∑ α i ⋅ l a b e l ( i ) = 0 \\sumα\_i·label^\{(i)\}=0 ∑α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值，可以：

1.  得到支持向量的个数：
    
        print(alphas[alphas > 0]) #[[0.14742212 0.17251816 0.04511926 0.36505954]]
2.  输出支持向量的个数
    
        print(shape(alphas[alphas > 0]))  # (1, 4)
3.  输出哪些是支持向量
    
        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 ) = e x p ( − ∣ ∣ x − y ∣ ∣ 2 2 σ 2 ) k(x,y)=exp(\\frac\{-||x-y||^2\}\{2σ^2\}) k(x,y)=exp(2σ2−∣∣x−y∣∣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))

--------------------

本章内容并没有完全吃懂，还需要返回来继续打磨…