机器学习实战第六章 支持向量机

支持向量机

什么是支持向量机

支持向量机(Support Vector Machine,SVM)是一种监督学习算法,用于进行二分类和回归分析。它的目标是找到一个最优的超平面(在二维空间中即为一条直线),将不同类别的样本点分隔开来,并且使得离超平面最近的样本点到该超平面的距离最大化。

  • 主要思想:支持向量机的主要思想是将原始数据映射到高维特征空间,并在该空间中寻找最优超平面,使得不同类别的样本点能够被最大间隔地分开。这种映射可以通过核函数(Kernel Function)来实现,它能够计算两个样本在高维空间中的内积,避免了直接进行昂贵的高维计算。
  • 补充说明:
    • 离超平面最近的样本点被称为支持向量,因为它们对于确定超平面起到了重要的支持作用。
    • 支持向量机可以通过引入软间隔和惩罚系数来处理非完全线性可分的问题。
    • 支持向量是指,空间中到间隔平面最近的点

问题数学化

由上述问题可知,我们的目标是要将间隔最大化,那么我们可以有以下公式来表示超平面:
w T x + b = 0 w^Tx+b = 0 wTx+b=0
其中w为法向量,b为位移项。
由于w和b等比例缩放并不会影响什么,我们可以由此假设支持向量到间隔平面的距离为1,那么两个异类支持向量到超平面的距离之和就是(此处略去一些证明,感兴趣的话可以看一下西瓜书,讲的非常好;但是没有说凸优化、对偶问题,这一部分尤其不理解) r = 2 ∣ ∣ w ∣ ∣ r = \frac{2}{||w||} r=∣∣w∣∣2
那么我们的目标就是
m a x w , b 2 w , s . t . y i ( w T x i + b ) ≥ 1 , i = 1 , 2 , . . . , m \underset{w,b}{max}\frac{2}{w}, \quad \quad \quad s.t.y_i(w^Tx_i+b)\geq1,\quad i = 1,2,...,m w,bmaxw2,s.t.yi(wTxi+b)1,i=1,2,...,m
也即 m i n w , b 1 2 ∣ ∣ w ∣ ∣ 2 , s . t . y i ( w T x i + b ) ≥ 1 , i = 1 , 2 , . . . , m \underset{w,b}{min}\frac{1}{2}||w||^2, \quad \quad \quad s.t.y_i(w^Tx_i+b)\geq1,\quad i = 1,2,...,m w,bmin21∣∣w2,s.t.yi(wTxi+b)1,i=1,2,...,m
为了求解该问题,可以使用拉格朗日乘子法的对偶问题来解决该问题(由于本人数学基础薄弱学了几天仍未完全理解,因此此处仅仅是记录下来,如若后续理解,立刻回来补充,感兴趣的话可以去参考西瓜书)。总之,我们最后得到的结果就是,我们要求解 m a x α ∑ i = 1 m α i − 1 2 ∑ i = 1 m ∑ j = 1 m α i α j y i y j x i T x j \underset{\alpha}{max}\sum_{i=1}^m\alpha_i-\frac{1}{2}\sum_{i=1}^m\sum_{j=1}^m\alpha_i\alpha_jy_iy_j{x_i}^Tx_j αmaxi=1mαi21i=1mj=1mαiαjyiyjxiTxj
其中 α \alpha α为拉格朗日乘子,然后约束条件有以下几条(包含了kkt条件,主要是这里搞不懂):
∑ i = 1 m a i y i = 0 α i ≥ 0 y i f ( x i ) − 1 ≥ 0 α i ( y i f ( x i ) − 1 ) = 0 \sum_{i=1}^ma_iy_i=0 \\ \alpha_i\geq0 \\ y_if(x_i)-1 \geq 0 \\ \alpha_i(y_if(x_i)-1)=0 i=1maiyi=0αi0yif(xi)10αi(yif(xi)1)=0
由于这个问题依旧很难解决,因此就有了SMO算法

SMO

根据上面给出的结论,只有 α \alpha α是变量,那么我们可以采取每次只优化一对 α \alpha α来逐渐逼近最终结果。
由于书上是先讲了一个简单的实现,再讲正式实现,因此我们也采取同样的方法以便于理解。

简化版SMO

前面我们知道了,每次只优化一对 α \alpha α来逐渐逼近最终结果,但是如何取 α \alpha α还是有讲究的(虽然不遵循这个规则应该也能得到正确的结果),我们主要采用以下两条规则:

  1. 选取的两个 α \alpha α必须要在边界之外
  2. 这两个 α \alpha α还没有进行过区间化处理或者不在边界上

此外,由于总有一些点无法正确分割,因此引入松弛变量的概念,这一概念能够允许一部分点出去分割面错误的一侧,体现在约束条件上就是约束条件变为:
C ≥ α ≥ 0 , 以及 ∑ i − 1 m α y i = 0 C\geq\alpha\geq0,以及\sum_{i-1}^m\alpha y_i = 0 Cα0,以及i1mαyi=0

此外还有一些需要注意的细节,下述代码中的eta是通过求解二次优化问题求解的来(数学基础薄弱,暂时无法理解)。同时,绝大部分代码都是为了满足约束而写的,下面将在代码中直接将需要注意的点标注出来:

# dataMatIn: 数据集,classLabels: 类别标签,C: 松弛变量,toler: 容错率,maxIter: 最大迭代次数
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()
    b = 0
    m, n = shape(dataMatrix)
    # 初始化alpha参数为0
    alphas = mat(zeros((m, 1)))
    iter = 0
    while iter < maxIter:
        alphasPairsChanged = 0  # 用于记录alpha是否已经进行优化
        for i in range(m):
            # ❶ 计算预测值与真实值的误差
            # fXi是预测的类别,Ei是误差。最后计算出来的Ei要拿去和toler比较,如果误差很大,就要对alpha进行优化
            fXi = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[i, :].T)) + b
            Ei = fXi - float(labelMat[i])
            # 优化alpha,更设定一定的容错率。若误差很大,且alpha值不在边界0或C上,则对alpha进行优化
            if ((labelMat[i] * Ei < -toler) and (alphas[i] < C)) or ((labelMat[i] * Ei > toler) and (alphas[i] > 0)):
                # ❷ 随机选择第⼆个alpha
                # 计算出误差,用于比较,看是否需要对alpha进行优化
                j = selectJrand(i, m)
                fXj = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[j, :].T)) + b
                Ej = fXj - float(labelMat[j])

                # 保存更新前的alpha值,后续需要用到来计算新的alpha值
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy()

                # ❸ 计算上下界L和H,用于将alpha[j]调整到0到C之间
                if labelMat[i] != labelMat[j]:
                    L = max(0, alphas[j] - alphas[i])
                    H = min(C, C + alphas[j] - alphas[i])
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, alphas[j] + alphas[i])

                # 如果L==H,则说明alpha[j]已经在边界上,不需要再优化了,直接执行下一次循环
                if L == H:
                    print("L==H")
                    continue

                # ❹ 计算eta,用于计算新的alpha[j],这里的eta是alpha[j]的最优修改量,可以用数学上的推导来证明
                # 前文中提到的分解为一系列步骤的优化问题,每一步都可以通过求导来求解,这里就是求解alpha[j]的最优修改量
                eta = 2.0 * dataMatrix[i, :] * dataMatrix[j, :].T - \
                      dataMatrix[i, :] * dataMatrix[i, :].T - \
                      dataMatrix[j, :] * dataMatrix[j, :].T 
                
                # 如果eta>=0,则不做任何改变,直接执行下一次循环,因为eta是alpha[j]的最优修改量,如果eta>=0,说明alpha[j]的改变量太⼩,不值得优化
                if eta >= 0:
                    print("eta>=0")
                    continue
                    
                # ❺ 更新alpha[j],由于更新后的alpha[j]必须满⾜0<=alpha[j]<=C,所以紧接着要对alpha[j]进行调整
                alphas[j] -= labelMat[j] * (Ei - Ej) / eta
                alphas[j] = clipAlpha(alphas[j], H, L)
                # 如果alpha[j]的改变量太⼩,则不做任何改变,直接执行下一次循环
                if abs(alphas[j] - alphaJold) < 0.00001:
                    print("j not moving enough")
                    continue
                    
                # ❻ 更新alpha[i],并同时更新b,这里的更新alpha[i]的公式是由前文提到的alpha[i]和alpha[j]的约束条件推导出来的,即他们进行一些运算后,和为常数
                alphas[i] += labelMat[j] * labelMat[i] * (alphaJold - alphas[j])
                
                # 更新b,但是要先计算出两个常数项b1和b2,这两个常数项的意义是:如果alpha[i]和alpha[j]都满⾜0<alpha<C,那么b1和b2是相等的
                # 如果alpha[i]或alpha[j]有⼀个等于0或C,那么b1和b2是不相等的,但是他们的平均值是相等的
                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
                
                # 如果0<alpha[i]<C,则b=b1;如果0<alpha[j]<C,则b=b2;如果alpha[i]和alpha[j]都等于0或者都等于C,则b=(b1+b2)/2
                if (0 < alphas[i]) and (C > alphas[i]):  # 如果0<alpha[i]<C,则b=b1
                    b = b1
                elif (0 < alphas[j]) and (C > alphas[j]):  # 如果0<alpha[j]<C,则b=b2
                    b = b2
                else:  # 如果alpha[i]和alpha[j]都等于0或者都等于C,则b=(b1+b2)/2
                    b = (b1 + b2) / 2.0
                alphasPairsChanged += 1  # 记录alpha是否已经进行优化
                print("iter: %d i:%d, pairs changed %d" % (iter, i, alphasPairsChanged))
                
        # ❼ 如果alpha没有更新,即alphasPairsChanged==0,则进行下一次循环,如果alpha有更新,则将iter置0后继续循环
        if alphasPairsChanged == 0:  # 如果alpha没有更新,则迭代次数+1,继续下一次循环
            iter += 1
        else:  # 如果alpha有更新,则迭代次数置0,继续循环
            iter = 0
        print("iteration number: %d" % iter)
    return b, alphas  # 返回常数项b和alpha向量

最终跑出来的结果如下图所示,由于数据集刚好可分,因此分割线画的很漂亮。此外还有一些其他的辅助函数此处不再给出,在本文末尾会给出完整代码。
在这里插入图片描述

完整版SMO

值得注意的是书上的代码是错误的,因此后面的代码将会有些许不同。

此外,完整版的SMO在原理上与简易版并没有什么不同,主要是在实现上有一些不同。具体来说有以下几点:

  1. 在选取第一个 α \alpha α方面:简易版采取遍历所有 α \alpha α的方式来选取第一个 α \alpha α,而完整版则是从非边界的 α \alpha α值中选取。就实现层面而言,是建立了一个 α \alpha α的列表,并跳过那些已经处于边界(即值为0或C)的 α \alpha α
  2. 在选取第二个 α \alpha α方面:原本是随机选取第二个 α \alpha α,现在则是通过最大化步长来实现。具体来讲,就是用eCache存储误差,然后选取误差最大的来进行下一步计算

由于鄙人不才,因为我感觉计算出来的结果不太对,主要疑惑的点在于这个eCache究竟有什么用。按照作者的说法,eCache是用来存储误差的,但是奇怪的点就在于每次修改 α \alpha α后误差应该会有变化才对,原本的误差应该就失效了,但是没有看到类似的语句。然后updateEk这个函数同样每次只更新一个值的误差,无法理解这个东西有什么用。

举个简单的例子:在这里插入图片描述
这段代码是为了测试更新i前后误差j的变化量的,根据下图的运行结果可以看到误差很大:在这里插入图片描述

但是我已经在这一章耗了太久了,留着后面完全搞懂了再回来修改,下面我会将完整代码以及结果贴出来,欢迎诸位批评指正。

import matplotlib.pyplot as plt
import numpy as np
import random


class optStruct:
    '''数据结构保存数据'''

    def __init__(self, dataMatIn, classLabels, C, toler):  # 初始化参数
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = np.shape(dataMatIn)[0]
        self.alphas = np.mat(np.zeros((self.m, 1)))
        self.b = 0
        self.eCache = np.mat(np.zeros((self.m, 2)))  # 第一列为有效标志flag


def loadDataSet(fileName):
    """
    读取数据
    Parameters:
        fileName - 文件名
    Returns:
        dataMat - 数据矩阵
        labelMat - 数据标签
    """
    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 calcEk(oS, k):
    '''计算预测误差'''
    fXk = float(np.multiply(oS.alphas, oS.labelMat).T * (oS.X * oS.X[k, :].T) + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek


def selectJrand(i, m):
    """
    函数说明:随机选择alpha_j的索引值

    Parameters:
        i - alpha_i的索引值
        m - alpha参数个数
    Returns:
        j - alpha_j的索引值
    """
    j = i  # 选择一个不等于i的j
    while (j == i):
        j = int(random.uniform(0, m))
    return j


def selectJ(i, oS, Ei):
    '''
     选择第二个alpha,即内循环的alpha值
     依据最大步长(max(abs(Ei - Ej)))选择
     '''
    maxK = -1;
    maxDeltaE = 0;
    Ej = 0
    # 存储Ei,第一位为有效标记
    oS.eCache[i] = [1, Ei]
    # # 构建一个非零表,返回非零E值所对应的index
    validEcacheList = np.nonzero(oS.eCache[:, 0].A)[0]

    if (len(validEcacheList)) > 1:
        for k in validEcacheList:
            if k == i:
                continue
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k
                maxDeltaE = deltaE
                Ej = Ek
        return maxK, Ej
    else:  # 第一次循环随机选一个alpha值
        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]


def clipAlpha(aj, H, L):
    """
    修剪alpha_j
    Parameters:
        aj - alpha_j的值
        H - alpha上限
        L - alpha下限
    Returns:
        aj - 修剪后的alpah_j的值
    """
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj



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
        #        eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j] #changed for kernel
        if eta >= 0: print("eta>=0"); return 0
        oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej) / eta
        oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
        updateEk(oS, j)  # 更新误差缓存
        if (abs(oS.alphas[j] - alphaJold) < 0.00001): print("j not moving enough"); return 0
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j])
        updateEk(oS, i)  # 更新误差缓存
        b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.X[i, :] * oS.X[i, :].T - oS.labelMat[j] * (
                    oS.alphas[j] - alphaJold) * oS.X[i, :] * oS.X[j, :].T
        b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.X[i, :] * oS.X[j, :].T - oS.labelMat[j] * (
                    oS.alphas[j] - alphaJold) * oS.X[j, :] * oS.X[j, :].T
        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]):
            oS.b = b1
        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]):
            oS.b = b2
        else:
            oS.b = (b1 + b2) / 2.0
        # alpha值发生变动返回1,否则返回0
        return 1
    else:
        return 0


def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):  # full Platt SMO
    '''
    主函数--外循环
    '''
    # 用构建的数据结构容纳所有数据
    oS = optStruct(np.mat(dataMatIn), np.mat(classLabels).transpose(), C, toler)
    iter = 0
    entireSet = True
    alphaPairsChanged = 0
    # 第一个alpha值的选择会在两种方式之间进行交替
    # 一种方式是在所有数据集上进行单遍扫描,另一种方式则是在非边界alpha中实现单遍扫描
    # 这里非边界指的是那些不等于边界0或C的alpha值
    # 同时,这里会跳过那些已知的不会改变的alpha值
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:
            # 在数据集上遍历任意可能的alpha
            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:  # 遍历非边界alpha值
            nonBoundIs = np.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  # 控制循环的flag
        elif (alphaPairsChanged == 0):
            entireSet = True
        print("iteration number: %d" % iter)
    return oS.b, oS.alphas


def showClassifer(dataMat, classLabels, w, b):
    """
    分类结果可视化
    Parameters:
        dataMat - 数据矩阵
        w - 直线法向量
        b - 直线解决
    Returns:
        无
    """
    # 绘制样本点
    data_plus = []  # 正样本
    data_minus = []  # 负样本
    for i in range(len(dataMat)):
        if classLabels[i] > 0:
            data_plus.append(dataMat[i])
        else:
            data_minus.append(dataMat[i])
    data_plus_np = np.array(data_plus)  # 转换为numpy矩阵
    data_minus_np = np.array(data_minus)  # 转换为numpy矩阵
    plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1], s=30, alpha=0.7)  # 正样本散点图
    plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1], s=30, alpha=0.7)  # 负样本散点图
    # 绘制直线
    x1 = max(dataMat)[0]
    x2 = min(dataMat)[0]
    a1, a2 = w
    b = float(b)
    a1 = float(a1[0])
    a2 = float(a2[0])
    y1, y2 = (-b - a1 * x1) / a2, (-b - a1 * x2) / a2
    plt.plot([x1, x2], [y1, y2])
    # 找出支持向量点
    for i, alpha in enumerate(alphas):
        if abs(alpha) > 0:
            print(alpha)
            x, y = dataMat[i]
            plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
    plt.show()


def calcWs(alphas, dataArr, classLabels):
    '''
    实际计算出的alpha值为0,而非零alpha所对应的也就是支持向量
    本计算函数遍历所有alpha,但最终起作用的只有支持向量
    '''
    X = np.mat(dataArr)
    labelMat = np.mat(classLabels).transpose()
    m, n = np.shape(X)
    w = np.zeros((n, 1))
    for i in range(m):
        w += np.multiply(alphas[i] * labelMat[i], X[i, :].T)
    return w


if __name__ == '__main__':
    random.seed(5)
    dataArr, classLabels = loadDataSet('testSet.txt')
    b, alphas = smoP(dataArr, classLabels, 0.6, 0.001, 100)
    w = calcWs(alphas, dataArr, classLabels)
    showClassifer(dataArr, classLabels, w, b)

运行结果:
在这里插入图片描述
从这里就能够理解我为什么说感觉计算结果不太对了,黄色区域中间的那个支持向量,不太可能成为支持向量才对。但是网上似乎没人讨论这个问题,而且书上的图也有几个并不是最近的点成为了支持向量,我有点怀疑是不是我自己搞错了。不过值得一提的是,他的分割线看起来是没什么问题的,这个倒是可以接受。

核函数

假定我们前面的代码是没有问题的,那么在二维平面上,它也无法将这样子的点用一条直线分开来:
在这里插入图片描述
看过三体的同志都知道,低纬度不可能做到的事情,高维未必就做不到,上面这些点,如果用核函数升维之后,可以看作类似以下样子:请添加图片描述

可以看到,它应该可以用一个平面来分割(平面升维为空间,那么直线就要升维为平面)

此处我们采用径向核函数来进行下一步的计算,其能够根据向量输出一个标量,其高斯版本公式为:
k ( x , y ) = e x p ( − ∣ ∣ x − y ∣ ∣ 2 2 σ 2 ) k(x,y) = exp(\frac{-||x-y||^2}{2\sigma^2}) k(x,y)=exp(2σ2∣∣xy2)

核函数处理后的数据与原始数据的运算没有什么区别,所以此处只贴出关键代码,完整代码参见本文结尾:

def kernelTrans(X, A, kTup):
    m, n = np.shape(X)
    K = np.mat(np.zeros((m, 1)))
    if kTup[0] == 'lin':
        K = X * A.T
    elif kTup[0] == 'rbf':
        for j in range(m):
            deltaRow = X[j, :] - A
            K[j] = deltaRow * deltaRow.T
        # ❶ 元素间的除法
        K = np.exp(K / (-1 * kTup[1] ** 2))
    else:
        raise NameError('Houston We Have a Problem -- That Kernel is not recognized')
    return K


class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler, kTup):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = np.shape(dataMatIn)[0]
        self.alphas = np.mat(np.zeros((self.m, 1)))
        self.b = 0
        self.eCache = np.mat(np.zeros((self.m, 2)))
        self.K = np.mat(np.zeros((self.m, self.m)))
        for i in range(self.m):
            self.K[:, i] = kernelTrans(self.X, self.X[i, :], kTup)


def innerL(i, oS):
    '''
    内循环
    :param i: 选择第i个alpha进行优化
    :param oS: 数据结构
    :return:
    '''
    Ei = calcEk(oS, i)
    # 违背KKT条件
    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)):
        # 启发式选择第二个alpha
        j, Ej = selectJ(i, oS, Ei)
        # 保存更新前的alpha值,使用深拷贝
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()
        # 保证alpha在0和C之间
        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])
        # 如果上下界相等,不做任何改变,直接返回0
        if L == H:
            print("L==H")
            return 0
        # 计算eta
        eta = 2.0 * oS.K[i, j] - oS.K[i, i] - oS.K[j, j]  # changed for kernel
        if eta >= 0:
            print("eta>=0")
            return 0
        # 更新alpha[j]
        oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej) / eta
        # 修剪alpha[j]
        oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
        # 更新Ej至误差缓存
        updateEk(oS, j)  # added this for the Ecache
        # alpha[j]变化太小
        if (abs(oS.alphas[j] - alphaJold) < 0.00001):
            print("j not moving enough")
            return 0
        # 更新alpha[i]
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j])
        # 更新Ei至误差缓存
        updateEk(oS, i)  # added this for the Ecache
        # 更新b1和b2
        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]
        # 根据b1和b2更新b
        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
        # alpha值发生变动返回1,否则返回0
        return 1
    else:
        return 0


运行结果如下图所示:
在这里插入图片描述
可以看到支持向量多数都分布在数据边缘,不规则是因为映射到高维空间后的样子未必与现在相同,没有分割线的原因是二维平面没有办法展示出高维平面。
在这里插入图片描述
同时可以看到,错误率处于一个比较低的水准。

实验

下面我们做一个手写识别的实验来加深理解。

主要需要修改的部分如下:


def img2vector(filename):
    '''
    将图片转换为向量
    '''
    returnVect = np.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)  # 加载目录下的所有文件
    m = len(trainingFileList)  # 文件个数
    trainingMat = np.zeros((m, 1024))  # 初始化数据集
    for i in range(m):
        fileNameStr = trainingFileList[i]  # 获得文件名
        fileStr = fileNameStr.split('.')[0]  # 获得不带后缀的文件名
        classNumStr = int(fileStr.split('_')[0])  # 获得分类数字
        if classNumStr == 9:  # 如果是9,则设为-1
            hwLabels.append(-1)
        else:  # 否则设为1
            hwLabels.append(1)
        trainingMat[i, :] = img2vector('%s/%s' % (dirName, fileNameStr))  # 加载数据集
    return trainingMat, hwLabels




def testDigits(kTup=('rbf', 10)):
    '''
    手写数字分类测试
    Parameters:
        kTup - 包含核函数信息的元组
    Returns:
        无
    '''
    dataArr, labelArr = loadImages('trainingDigits')  # 加载训练集
    b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 100, kTup)  # 计算b和alphas
    datMat = np.mat(dataArr)
    labelMat = np.mat(labelArr).transpose()
    svInd = np.nonzero(alphas.A > 0)[0]  # 获得支持向量
    sVs = datMat[svInd]  # 支持向量的特征数据
    labelSV = labelMat[svInd]  # 支持向量的类别
    print("there are %d Support Vectors" % np.shape(sVs)[0])
    m, n = np.shape(datMat)
    errorCount = 0
    for i in range(m):
        # 计算各个点的核
        kernelEval = kernelTrans(sVs, datMat[i, :], kTup)
        # 根据支持向量的点,计算超平面,返回预测结果
        predict = kernelEval.T * np.multiply(labelSV, alphas[svInd]) + b
        # 返回数组中各元素的正负符号,用1和-1表示,并统计错误个数
        if np.sign(predict) != np.sign(labelArr[i]):
            errorCount += 1
    # 打印错误率
    print("the training error rate is: %f" % (float(errorCount) / m))
    # 加载测试集
    dataArr, labelArr = loadImages('testDigits')
    errorCount = 0
    datMat = np.mat(dataArr)
    labelMat = np.mat(labelArr).transpose()
    m, n = np.shape(datMat)
    # 循环测试集
    for i in range(m):
        # 计算各个点的核
        kernelEval = kernelTrans(sVs, datMat[i, :], kTup)
        # 根据支持向量的点,计算超平面,返回预测结果
        predict = kernelEval.T * np.multiply(labelSV, alphas[svInd]) + b
        # 返回数组中各元素的正负符号,用1和-1表示,并统计错误个数
        if np.sign(predict) != np.sign(labelArr[i]):
            errorCount += 1
    # 打印错误率
    print("the test error rate is: %f" % (float(errorCount) / m))


运行结果如下:
在这里插入图片描述
可以看到准确率很高,同时书上说修改径向核函数的参数会对准确率有很大的影响,但是经本人测试,并没有发现准确率变化很大。

完整代码

简化版SMO

from numpy import *
import matplotlib.pyplot as plt


def loadDataSet(fileName):
    dataMat = []
    labelMat = []
    fr = open(fileName)
    for line in fr.readlines():  # 逐行读取
        lineArr = line.strip().split('\t')  # 去除首尾空白字符,按制表符分割
        dataMat.append([float(lineArr[0]), float(lineArr[1])])  # 添加数据
        labelMat.append(float(lineArr[2]))  # 添加标签
    return dataMat, labelMat


def selectJrand(i, m):
    j = i
    while j == i:  # 随机选择一个与i不同的j
        j = int(random.uniform(0, m))
    return j


def clipAlpha(aj, H, L):  # 调整aj的值,使其处于 L<=aj<=H
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj


# dataMatIn: 数据集,classLabels: 类别标签,C: 松弛变量,toler: 容错率,maxIter: 最大迭代次数
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    dataMatrix = mat(dataMatIn)
    labelMat = mat(classLabels).transpose()
    b = 0
    m, n = shape(dataMatrix)
    # 初始化alpha参数为0
    alphas = mat(zeros((m, 1)))
    iter = 0
    while iter < maxIter:
        alphasPairsChanged = 0  # 用于记录alpha是否已经进行优化
        for i in range(m):
            # ❶ 计算预测值与真实值的误差
            # fXi是预测的类别,Ei是误差。最后计算出来的Ei要拿去和toler比较,如果误差很大,就要对alpha进行优化
            fXi = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[i, :].T)) + b
            Ei = fXi - float(labelMat[i])
            # 优化alpha,更设定一定的容错率。若误差很大,且alpha值不在边界0或C上,则对alpha进行优化
            if ((labelMat[i] * Ei < -toler) and (alphas[i] < C)) or ((labelMat[i] * Ei > toler) and (alphas[i] > 0)):
                # ❷ 随机选择第⼆个alpha
                # 计算出误差,用于比较,看是否需要对alpha进行优化
                j = selectJrand(i, m)
                fXj = float(multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[j, :].T)) + b
                Ej = fXj - float(labelMat[j])

                # 保存更新前的alpha值,后续需要用到来计算新的alpha值
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy()

                # ❸ 计算上下界L和H,用于将alpha[j]调整到0到C之间
                if labelMat[i] != labelMat[j]:
                    L = max(0, alphas[j] - alphas[i])
                    H = min(C, C + alphas[j] - alphas[i])
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, alphas[j] + alphas[i])

                # 如果L==H,则说明alpha[j]已经在边界上,不需要再优化了,直接执行下一次循环
                if L == H:
                    print("L==H")
                    continue

                # ❹ 计算eta,用于计算新的alpha[j],这里的eta是alpha[j]的最优修改量,可以用数学上的推导来证明
                # 前文中提到的分解为一系列步骤的优化问题,每一步都可以通过求导来求解,这里就是求解alpha[j]的最优修改量
                eta = 2.0 * dataMatrix[i, :] * dataMatrix[j, :].T - \
                      dataMatrix[i, :] * dataMatrix[i, :].T - \
                      dataMatrix[j, :] * dataMatrix[j, :].T

                # 如果eta>=0,则不做任何改变,直接执行下一次循环,因为eta是alpha[j]的最优修改量,如果eta>=0,说明alpha[j]的改变量太⼩,不值得优化
                if eta >= 0:
                    print("eta>=0")
                    continue

                # ❺ 更新alpha[j],由于更新后的alpha[j]必须满⾜0<=alpha[j]<=C,所以紧接着要对alpha[j]进行调整
                alphas[j] -= labelMat[j] * (Ei - Ej) / eta
                alphas[j] = clipAlpha(alphas[j], H, L)
                # 如果alpha[j]的改变量太⼩,则不做任何改变,直接执行下一次循环
                if abs(alphas[j] - alphaJold) < 0.00001:
                    print("j not moving enough")
                    continue

                # ❻ 更新alpha[i],并同时更新b,这里的更新alpha[i]的公式是由前文提到的alpha[i]和alpha[j]的约束条件推导出来的,即他们进行一些运算后,和为常数
                alphas[i] += labelMat[j] * labelMat[i] * (alphaJold - alphas[j])

                # 更新b,但是要先计算出两个常数项b1和b2,这两个常数项的意义是:如果alpha[i]和alpha[j]都满⾜0<alpha<C,那么b1和b2是相等的
                # 如果alpha[i]或alpha[j]有⼀个等于0或C,那么b1和b2是不相等的,但是他们的平均值是相等的
                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

                # 如果0<alpha[i]<C,则b=b1;如果0<alpha[j]<C,则b=b2;如果alpha[i]和alpha[j]都等于0或者都等于C,则b=(b1+b2)/2
                if (0 < alphas[i]) and (C > alphas[i]):  # 如果0<alpha[i]<C,则b=b1
                    b = b1
                elif (0 < alphas[j]) and (C > alphas[j]):  # 如果0<alpha[j]<C,则b=b2
                    b = b2
                else:  # 如果alpha[i]和alpha[j]都等于0或者都等于C,则b=(b1+b2)/2
                    b = (b1 + b2) / 2.0
                alphasPairsChanged += 1  # 记录alpha是否已经进行优化
                print("iter: %d i:%d, pairs changed %d" % (iter, i, alphasPairsChanged))

        # ❼ 如果alpha没有更新,即alphasPairsChanged==0,则进行下一次循环,如果alpha有更新,则将iter置0后继续循环
        if alphasPairsChanged == 0:  # 如果alpha没有更新,则迭代次数+1,继续下一次循环
            iter += 1
        else:  # 如果alpha有更新,则迭代次数置0,继续循环
            iter = 0
        print("iteration number: %d" % iter)
    return b, alphas  # 返回常数项b和alpha向量


# ❸ 计算w
def calcWs(alphas, dataArr, classLabels):
    X = mat(dataArr)
    labelMat = mat(classLabels).transpose()
    m, n = shape(X)
    w = zeros((n, 1))  # 初始化w
    for i in range(m):  # 计算w
        w += multiply(alphas[i] * labelMat[i], X[i, :].T)
    return w


#  ❹ 画出数据集和SVM分类直线的函数
def showClassifer(dataMat, w, b, labelMat):

    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    data_plus = []  # 正样本
    data_minus = []  # 负样本
    for i in range(len(dataMat)):  # 分类
        if labelMat[i] > 0:
            data_plus.append(dataMat[i])
        else:
            data_minus.append(dataMat[i])
    data_plus_np = array(data_plus)  # 转换为numpy矩阵
    data_minus_np = array(data_minus)  # 转换为numpy矩阵
    plt.scatter(transpose(data_plus_np)[0], transpose(data_plus_np)[1], s=30, alpha=0.7)  # 正样本散点图
    plt.scatter(transpose(data_minus_np)[0], transpose(data_minus_np)[1], s=30, alpha=0.7)  # 负样本散点图
    x1 = max(dataMat)[0]  # x轴最大值
    x2 = min(dataMat)[0]  # x轴最小值
    a1, a2 = w  # 获得w的两个参数
    b = float(b)
    a1 = float(a1[0])
    a2 = float(a2[0])
    y1, y2 = (-b - a1 * x1) / a2, (-b - a1 * x2) / a2  # 根据x1、x2计算y1、y2
    plt.plot([x1, x2], [y1, y2])  # 画出直线
    for i, alpha in enumerate(alphas):  # 找出支持向量点
        if abs(alpha) > 0:
            x, y = dataMat[i]
            plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
    plt.show()


dataArr, classLabels = loadDataSet('testSet.txt')

b, alphas = smoSimple(dataArr, classLabels, 0.6, 0.001, 40)

w = calcWs(alphas, dataArr, classLabels)

showClassifer(dataArr, w, b, classLabels)

核函数版SMO

import matplotlib.pyplot as plt
import numpy as np
import random


def kernelTrans(X, A, kTup):
    m, n = np.shape(X)
    K = np.mat(np.zeros((m, 1)))
    if kTup[0] == 'lin':
        K = X * A.T
    elif kTup[0] == 'rbf':
        for j in range(m):
            deltaRow = X[j, :] - A
            K[j] = deltaRow * deltaRow.T
        # ❶ 元素间的除法
        K = np.exp(K / (-1 * kTup[1] ** 2))
    else:
        raise NameError('Houston We Have a Problem -- That Kernel is not recognized')
    return K


class optStruct:
    def __init__(self, dataMatIn, classLabels, C, toler, kTup):
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = np.shape(dataMatIn)[0]
        self.alphas = np.mat(np.zeros((self.m, 1)))
        self.b = 0
        self.eCache = np.mat(np.zeros((self.m, 2)))
        self.K = np.mat(np.zeros((self.m, self.m)))
        for i in range(self.m):
            self.K[:, i] = kernelTrans(self.X, self.X[i, :], kTup)


def loadDataSet(fileName):
    """
    读取数据
    Parameters:
        fileName - 文件名
    Returns:
        dataMat - 数据矩阵
        labelMat - 数据标签
    """
    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 calcEk(oS, k):
    """
    函数说明:计算误差

    Parameters:
        oS - 数据结构
        k - 标号为k的数据
    Returns:
        Ek - 标号为k的数据误差
    """
    fXk = float(np.multiply(oS.alphas, oS.labelMat).T * oS.K[:, k] + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek


def selectJrand(i, m):
    """
    函数说明:随机选择alpha_j的索引值

    Parameters:
        i - alpha_i的索引值
        m - alpha参数个数
    Returns:
        j - alpha_j的索引值
    """
    j = i  # 选择一个不等于i的j
    while (j == i):
        j = int(random.uniform(0, m))
    return j


def selectJ(i, oS, Ei):
    '''
     选择第二个alpha,即内循环的alpha值
     依据最大步长(max(abs(Ei - Ej)))选择
     '''
    maxK = -1;
    maxDeltaE = 0;
    Ej = 0
    # 存储Ei,第一位为有效标记
    oS.eCache[i] = [1, Ei]
    # # 构建一个非零表,返回非零E值所对应的index
    validEcacheList = np.nonzero(oS.eCache[:, 0].A)[0]

    if (len(validEcacheList)) > 1:
        for k in validEcacheList:
            if k == i:
                continue
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k
                maxDeltaE = deltaE
                Ej = Ek
        return maxK, Ej
    else:  # 第一次循环随机选一个alpha值
        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]


def clipAlpha(aj, H, L):
    """
    修剪alpha_j
    Parameters:
        aj - alpha_j的值
        H - alpha上限
        L - alpha下限
    Returns:
        aj - 修剪后的alpah_j的值
    """
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj



def innerL(i, oS):
    '''
    内循环
    :param i: 选择第i个alpha进行优化
    :param oS: 数据结构
    :return:
    '''
    Ei = calcEk(oS, i)
    # 违背KKT条件
    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)):
        # 启发式选择第二个alpha
        j, Ej = selectJ(i, oS, Ei)
        # 保存更新前的alpha值,使用深拷贝
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()
        # 保证alpha在0和C之间
        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])
        # 如果上下界相等,不做任何改变,直接返回0
        if L == H:
            print("L==H")
            return 0
        # 计算eta
        eta = 2.0 * oS.K[i, j] - oS.K[i, i] - oS.K[j, j]  # changed for kernel
        if eta >= 0:
            print("eta>=0")
            return 0
        # 更新alpha[j]
        oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej) / eta
        # 修剪alpha[j]
        oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
        # 更新Ej至误差缓存
        updateEk(oS, j)  # added this for the Ecache
        # alpha[j]变化太小
        if (abs(oS.alphas[j] - alphaJold) < 0.00001):
            print("j not moving enough")
            return 0
        # 更新alpha[i]
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j])
        # 更新Ei至误差缓存
        updateEk(oS, i)  # added this for the Ecache
        # 更新b1和b2
        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]
        # 根据b1和b2更新b
        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
        # alpha值发生变动返回1,否则返回0
        return 1
    else:
        return 0


def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):  # full Platt SMO
    '''
    主函数--外循环
    '''
    # 用构建的数据结构容纳所有数据
    oS = optStruct(np.mat(dataMatIn), np.mat(classLabels).transpose(), C, toler, kTup)
    iter = 0
    entireSet = True
    alphaPairsChanged = 0
    # 第一个alpha值的选择会在两种方式之间进行交替
    # 一种方式是在所有数据集上进行单遍扫描,另一种方式则是在非边界alpha中实现单遍扫描
    # 这里非边界指的是那些不等于边界0或C的alpha值
    # 同时,这里会跳过那些已知的不会改变的alpha值
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:
            # 在数据集上遍历任意可能的alpha
            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:  # 遍历非边界alpha值
            nonBoundIs = np.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  # 控制循环的flag
        elif (alphaPairsChanged == 0):
            entireSet = True
        print("iteration number: %d" % iter)
    return oS.b, oS.alphas


def showClassifer(dataMat, classLabels, alphas):
    """
    分类结果可视化
    Parameters:
        dataMat - 数据矩阵
        w - 直线法向量
        b - 直线解决
    Returns:
        无
    """
    # 绘制样本点
    data_plus = []  # 正样本
    data_minus = []  # 负样本
    for i in range(len(dataMat)):
        if classLabels[i] > 0:
            data_plus.append(dataMat[i])
        else:
            data_minus.append(dataMat[i])
    data_plus_np = np.array(data_plus)  # 转换为numpy矩阵
    data_minus_np = np.array(data_minus)  # 转换为numpy矩阵
    plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1], s=30, alpha=0.7)  # 正样本散点图
    plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1], s=30, alpha=0.7)  # 负样本散点图
    # 找出支持向量点
    for i, alpha in enumerate(alphas):
        if abs(alpha) > 0:
            print(alpha)
            x, y = dataMat[i]
            plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
    plt.show()



def testRBF(k1=1.3):
    '''
    测试函数
    Parameters:
        k1 - 使用高斯核函数的时候表示到达率
    Returns:
        无
    '''
    dataArr, classLabels = loadDataSet('testSetRBF.txt')  # 加载训练集
    b, alphas = smoP(dataArr, classLabels, 200, 0.0001, 100, ('rbf', k1))  # 根据训练集计算b和alphas
    showClassifer(dataArr, classLabels, alphas)  # 画出分类图
    datMat = np.mat(dataArr)
    labelMat = np.mat(classLabels).transpose()
    # 获得支持向量
    svInd = np.nonzero(alphas.A > 0)[0]
    sVs = datMat[svInd]  # 支持向量的特征数据
    labelSV = labelMat[svInd]  # 支持向量的类别
    print("there are %d Support Vectors" % np.shape(sVs)[0])
    m, n = np.shape(datMat)
    errorCount = 0
    for i in range(m):
        # 计算各个点的核
        kernelEval = kernelTrans(sVs, datMat[i, :], ('rbf', k1))
        # 根据支持向量的点,计算超平面,返回预测结果
        predict = kernelEval.T * np.multiply(labelSV, alphas[svInd]) + b
        # 返回数组中各元素的正负符号,用1和-1表示,并统计错误个数
        if np.sign(predict) != np.sign(classLabels[i]):
            errorCount += 1
    # 打印错误率
    print("the training error rate is: %f" % (float(errorCount) / m))
    # 加载测试集
    dataArr, classLabels = loadDataSet('testSetRBF2.txt')
    errorCount = 0
    datMat = np.mat(dataArr)
    labelMat = np.mat(classLabels).transpose()
    m, n = np.shape(datMat)
    # 循环测试集
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], ('rbf', k1))
        # 根据支持向量的点,计算超平面,返回预测结果
        predict = kernelEval.T * np.multiply(labelSV, alphas[svInd]) + b
        # 返回数组中各元素的正负符号,用1和-1表示,并统计错误个数
        if np.sign(predict) != np.sign(classLabels[i]):
            errorCount += 1
        # 打印错误率
    print("the test error rate is: %f" % (float(errorCount) / m))


if __name__ == '__main__':
    testRBF()

  • 1
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值