支持向量机核函数的实现

一：回顾SVM中的SMO算法

https://www.cnblogs.com/ssyfj/p/13363526.html

二：核函数的了解

（一）西瓜书（粗略了解）

（二）统计学习方法（详细）

（三）推文：支持向量机原理(三)线性不可分支持向量机与核函数

（四）推文：核函数和核矩阵

（五）机器学习实战（代码加强理解）

三：重点知识

对于在低维线性不可分的数据，在映射到了高维以后，就变成线性可分。这个映射过程分为隐式转换和显示转换两种情况。

（一）显式映射（x->显示构造高维特征φ(x)->内积计算φ(x).T@φ(x))

线性不可分的低维特征数据，我们可以将其映射到高维，就能线性可分。现在我们将它运用到我们的SVM的算法上。回顾线性可分SVM的优化目标函数：

注意到上式低维特征仅仅以内积xi∙xj的形式出现，如果我们定义一个低维特征空间到高维特征空间的映射ϕ(x)（比如从2维到5维的映射），将所有特征映射到一个更高的维度，让数据线性可分，我们就可以继续按前面的SMO方法来优化目标函数，求出分离超平面和分类决策函数了。

也就是说现在的SVM的优化目标函数变成：

可以看到，和线性可分SVM的优化目标函数的区别仅仅是将内积x_i∙x_j替换为ϕ(x_i)∙ϕ(x_j)。

这就是显式映射，我们需要真正意义上将数据映射到高维特征空间中去，然后对这个转换后的高维向量进行内积运算。

看起来似乎这样我们就已经完美解决了线性不可分SVM的问题了，但是事实是不是这样呢？

我们看看，假如是一个2维特征的数据，我们可以将其映射到5维来做特征的内积，如果原始空间是三维，可以映射到到19维空间，似乎还可以处理。

但是如果我们的低维特征是100个维度，1000个维度呢？

那么我们要将其映射到超级高的维度来计算特征的内积。

这时候映射成的高维维度是爆炸性增长的，这个计算量实在是太大了，而且如果遇到无穷维的情况，就根本无从计算了。

显示转换存在以下两个问题：

1.通常我们不知道映射ϕ的具体形式；
2.映射后的空间维数可能非常高，甚至无限维，直接计算开销太大，十分困难.

（二）隐式映射（核函数）---见统计学习方法p136

利用核函数，不需要将输入空间中的样本映射到新的空间中去，在输入空间中就可以直接计算出内积了。

对输入空间向高维空间的一种隐式映射（注意，这也是低维空间到高维空间的一种映射<理解意义上的>）。不需要给出那个映射，而在输入空间中就可以计算内积，即核函数方法。【核函数可以简化映射空间中的内积运算】

（三）显式映射与隐式映射案例对比（加深理解）---见西瓜书笔记

显示映射：需要构建高维空间，需要大量计算

隐式映射：使用核函数，直接将X作为输入，避免中间构造高维特征空间的计算。所有的计算都是在原始特征空间中进行的，并没有构建高维空间。

（四）核技巧---统计学习方法

在学习核预测中只定义核函数K(x,z),而不是显式的定义映射函数。因为直接计算K(x,z)比较容易，而通过φ(x)核φ(z)计算K(x,z)并不容易（也就是说明我们即便是想要通过显式映射到高维特征空间进行划分数据集并不容易）：

因为对于给定的核K(x,z)或者我们想要的高维中划分超平面选取，高维特征空间和映射函数的取法并不唯一，可以选取不同的特征空间，即便是在同一特征空间里，也可以取不同的映射

四：核函数代码实现

（一）实现核函数转换函数

#实现核函数（核转换函数）
def kernelTrans(data_X,data_Xi,kTup):
    """
    因为我们的核函数是对数据集进行的隐式映射，后面求解α等，并没有对数据集本身修改过，
    所以我们只需要调用一次核转换函数即可求得全局的核矩阵解
    :param data_X:  表示传入的全部数据集
    :param data_Xi: 表示传入的第i个样本数据
    :param kTup: 是一个包含核函数信息的元组
    :return: 返回核矩阵的第i列
    """
    m,n = data_X.shape
    K = np.zeros((m,1)) #这就是我们要返回的核矩阵的第i列初始状态

    #根据我们输入的核函数信息判断选取的是哪一个核函数，并获取传入参数
    if kTup[0] == 'lin':    #线性核函数，即没有使用核函数
        K = data_X@data_Xi  #(m,1)
    elif kTup[0] == 'rbf':     #高斯核函数
        for j in range(m):  #对数据集每一个样本进行高斯处理
            deltaRow = data_X[j,:] - data_Xi    #获取差值
            K[j] = deltaRow@deltaRow.T  #获取高斯函数分子部分，无论上面求解的deltaRow正负，这里内积以后，全是正值
        K = np.exp(K/(-1*kTup[1]**2))   #对整体K进行除方差操作，获取高斯处理后的值
    else:   #如果还有其他的核函数，接着使用elif即可
        raise NameError("the kernel is not recognized")#对于没有的核函数，我们抛出异常

    return K

（二）修改数据结构，获取全局核矩阵

#首先我们定义一个数据结构（类），来保存所有的重要值
class optStruct:
    def __init__(self,data_X,data_Y,C,toler,kTup):   #输入参数分别是数据集、类别标签、常数C用于软间隔、和容错率toler
        self.X = data_X
        self.label = data_Y
        self.C = C
        self.toler = toler  #就是软间隔中的ε，调节最大间隔大小
        self.m = data_X.shape[0]
        self.alphas = np.zeros(self.m)    #存放每个样本点的α值
        self.b = 0  #存放阈值
        self.eCache = np.zeros((self.m,2))  #用于缓存误差，每个样本点对应一个Ei值，第一列为标识符，标志是否为有效值，第二列存放有效值
        self.K = np.zeros((self.m,self.m))
        #这里获取全局核矩阵即可
        for i in range(self.m):
            self.K[:,i] = kernelTrans(data_X,data_X[i],kTup).flatten()

注意：核矩阵是对称阵，由（一）可以知道，我们求解的最终K核矩阵是对称阵，即半正定矩阵（实对称矩阵A称为半正定矩阵），所以我们是符合核矩阵要求的。

（三）修改内循环函数---因为这里面使用了核函数

#三：实现内循环函数，相比于外循环，这里包含了主要的更新操作
def innerL(i,oS):   #由外循环提供i值（具体选取要违背kkT<这里实现>,使用交替遍历<外循环中实现>）---提供α_1的索引
    Ei = calcEk(oS,i)   #计算E1值，主要是为了下面KKT条件需要使用到

    #如果下面违背了KKT条件，则正常进行α、Ek、b的更新,重点：后面单独说明下面是否满足违反KKT条件
    if ((oS.label[i]*Ei < -oS.toler) and (oS.alphas[i] < oS.C)) or\
       ((oS.label[i]*Ei > oS.toler) and (oS.alphas[i] > 0)):  #注意：对于硬间隔，我们直接和1对比，对于软间隔，我们要和1 +或- ε对比
        #开始在内循环中，选取差值最大的α_2下标索引
        j,Ej = selectJ(i,oS,Ei)
        #因为后面要修改α_1与α_2的值，但是后面修改阈值b的时候需要用到新旧两个值，所以我们需要在更新α值之前进行保存旧值
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()

        #分析约束条件（是对所有α都适用），一会对我们新的α_2进行截取纠正，注意：α_1是由α_2推出的，所以不需要进行验证了。
        #如果y_1!=y_2异号时：
        if oS.label[i] != oS.label[j]:
            L = max(0,alphaJold-alphaIold)
            H = min(oS.C,oS.C+alphaJold-alphaIold)
        else:   #如果y_1==y_2同号时
            L = max(0,alphaJold+alphaIold-oS.C)
            H = min(oS.C,alphaJold+alphaIold)
        #上面就是将α_j调整到L,H之间
        if L==H:    #如果L==H，之间返回0，跳出这次循环，不进行改变(单值选择，没必要）
            return 0

        #计算η值=k_11+k_22-2k_12
        eta = oS.K[i,i] + oS.K[j,j] - 2.0*oS.K[i,j]  #eta性质可以知道是>=0的，所以我们只需要判断是否为0即可
        if eta <= 0:
            print("eta <= 0")
            return 0

        #当上面所有条件都满足以后，我们开始正式修改α_2值，并更新对应的Ek值
        oS.alphas[j] += oS.label[j]*(Ei-Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        updateEk(oS,j)

        #查看α_2是否有足够的变化量，如果没有足够变化量，我们直接返回，不进行下面更新α_1,注意：因为α_2变化量较小，所以我们没有必要非得把值变回原来的旧值
        if abs(oS.alphas[j] - alphaJold) < 0.00001:
            print("J not move enough")
            return 0

        #开始更新α_1值,和Ek值
        oS.alphas[i] += oS.label[i]*oS.label[j]*(alphaJold-oS.alphas[j])
        updateEk(oS,i)

        #开始更新阈值b,正好使用到了上面更新的Ek值
        b1 = oS.b - Ei - oS.label[i] * (oS.alphas[i] - alphaIold) * oS.K[i,i] - oS.label[j] * (
                    oS.alphas[j] - alphaJold) * oS.K[i,j]

        b2 = oS.b - Ej - oS.label[i] * (oS.alphas[i] - alphaIold) * oS.K[i,j] - oS.label[j] * (
                    oS.alphas[j] - alphaJold) * oS.K[j,j]

        #根据统计学习方法中阈值b在每一步中都会进行更新，
        #1.当新值alpha_1不在界上时(0<alpha_1<C)，b_new的计算规则为：b_new=b1
        #2.当新值alpha_2不在界上时(0 < alpha_2 < C)，b_new的计算规则为：b_new = b2
        #3.否则当alpha_1和alpha_2都不在界上时，b_new = 1/2(b1+b2)
        if oS.alphas[i] > 0 and oS.alphas[i] < oS.C:
            oS.b = b1
        elif oS.alphas[j] > 0 and oS.alphas[j] < oS.C:
            oS.b = b2
        else:
            oS.b = 1/2*(b1+b2)

        #注意：这里我们应该根据b_new更新一次Ei，但是我们这里没有写，因为我们将这一步提前到了最开始，即selectJ中

        #以上全部更新完毕，开始返回标识
        return 1
    return 0    #没有违背KKT条件

（四）修改计算误差函数---含核函数

#计算每个样本点k的Ek值，就是计算误差值=预测值-标签值
def calcEk(oS,k):
    # 根据西瓜书6.24，我们可以知道预测值如何使用α值进行求解
    fxk = np.multiply(oS.alphas,oS.label).T@(oS.K[:,k])+oS.b #np.multiply之后还是(m,1),(oS.X@oS.X[k,:])之后是(m,1),通过转置(1,m)@(m,1)-->实数后+b即可得到预测值fx
    #获取误差值Ek
    Ek = fxk - oS.label[k]
    return Ek

（五）全部代码

import numpy as np
import matplotlib.pyplot as plt
from numpy import *

#一：SMO算法中的辅助函数
#加载数据集
def loadDataSet(filename):
    dataSet = np.loadtxt(filename)
    m,n = dataSet.shape
    data_X = dataSet[:,0:n-1]
    data_Y = dataSet[:,n-1]

    return data_X,data_Y

#随机选取一个数J，为后面内循环选取α_2做辅助（如果α选取不满足条件，就选择这个方法随机选取）
def selectJrand(i,m):   #主要就是根据α_1的索引i,从所有数据集索引中随机选取一个作为α_2的索引
    j = i
    while j==i:
        j = np.int(np.random.uniform(0,m))  #从0~m中随机选取一个数，是进行整数化的
    print("random choose index for α_2:%d"%(j))
    return j    #由于这里返回随机数，所以后面结果 可能导致不同

def clipAlpha(aj,H,L):  #根据我们的SVM算法中的约束条件的分析，我们对获取的aj,进行了截取操作
    if aj > H:
        aj = H
    if aj < L:
        aj = L
    return aj

#二：SMO的支持函数
#实现核函数（核转换函数）
def kernelTrans(data_X,data_Xi,kTup):
    """
    因为我们的核函数是对数据集进行的隐式映射，后面求解α等，并没有对数据集本身修改过，
    所以我们只需要调用一次核转换函数即可求得全局的核矩阵解
    :param data_X:  表示传入的全部数据集
    :param data_Xi: 表示传入的第i个样本数据
    :param kTup: 是一个包含核函数信息的元组
    :return: 返回核矩阵的第i列
    """
    m,n = data_X.shape
    K = np.zeros((m,1)) #这就是我们要返回的核矩阵的第i列初始状态

    #根据我们输入的核函数信息判断选取的是哪一个核函数，并获取传入参数
    if kTup[0] == 'lin':    #线性核函数，即没有使用核函数
        K = data_X@data_Xi  #(m,1)
    elif kTup[0] == 'rbf':     #高斯核函数
        for j in range(m):  #对数据集每一个样本进行高斯处理
            deltaRow = data_X[j,:] - data_Xi    #获取差值
            K[j] = deltaRow@deltaRow.T  #获取高斯函数分子部分
        K = np.exp(K/(-1*kTup[1]**2))   #对整体K进行除方差操作，获取高斯处理后的值
    else:   #如果还有其他的核函数，接着使用elif即可
        raise NameError("the kernel is not recognized")#对于没有的核函数，我们抛出异常

    return K

#首先我们定义一个数据结构（类），来保存所有的重要值
class optStruct:
    def __init__(self,data_X,data_Y,C,toler,kTup):   #输入参数分别是数据集、类别标签、常数C用于软间隔、和容错率toler
        self.X = data_X
        self.label = data_Y
        self.C = C
        self.toler = toler  #就是软间隔中的ε，调节最大间隔大小
        self.m = data_X.shape[0]
        self.alphas = np.zeros(self.m)    #存放每个样本点的α值
        self.b = 0  #存放阈值
        self.eCache = np.zeros((self.m,2))  #用于缓存误差，每个样本点对应一个Ei值，第一列为标识符，标志是否为有效值，第二列存放有效值
        self.K = np.zeros((self.m,self.m))
        #这里获取全局核矩阵即可
        for i in range(self.m):
            self.K[:,i] = kernelTrans(data_X,data_X[i],kTup).flatten()

#计算每个样本点k的Ek值，就是计算误差值=预测值-标签值
def calcEk(oS,k):
    # 根据西瓜书6.24，我们可以知道预测值如何使用α值进行求解
    fxk = np.multiply(oS.alphas,oS.label).T@(oS.K[:,k])+oS.b #np.multiply之后还是(m,1),(oS.X@oS.X[k,:])之后是(m,1),通过转置(1,m)@(m,1)-->实数后+b即可得到预测值fx
    #获取误差值Ek
    Ek = fxk - oS.label[k]
    return Ek

#内循环的启发式方法，获取最大差值|Ei-Ej|对应的Ej的索引J
def selectJ(i,oS,Ei):   #注意我们要传入第一个α对应的索引i和误差值Ei,后面会用到
    maxK = -1   #用于保存临时最大索引
    maxDeltaE = 0   #用于保存临时最大差值--->|Ei-Ej|
    Ej = 0  #保存我们需要的Ej误差值

    #重点：这里我们是把SMO最后一步（根据最新阈值b，来更新Ei)提到第一步来进行了，所以这一步是非常重要的
    oS.eCache[i] = [1,Ei]

    #开始获取各个Ek值，比较|Ei-Ej|获取Ej的所有
    #获取所有有效的Ek值对应的索引
    validECacheList = np.where(oS.eCache[:,0]!=0)[0]    #根据误差缓存中第一列非0，获取对应的有效误差值
    if len(validECacheList) > 1:   #如果有效误差缓存长度大于1(因为包括Ei),则正常进行获取j值，否则使用selectJradn方法选取一个随机J值
        for k in validECacheList:
            if k == i:  #相同则不处理
                continue
            #开始计算Ek值，进行对比，获取最大差值
            Ek = calcEk(oS,k)
            deltaE = abs(Ei - Ek)
            if deltaE > maxDeltaE:  #更新Ej及其索引位置
                maxK = k
                maxDeltaE = deltaE
                Ej = Ek
        return maxK,Ej  #返回我们找到的第二个变量α_2的位置
    else:   #没有有效误差缓存，则随机选取一个索引，进行返回
        j = selectJrand(i,oS.m)
        Ej = calcEk(oS,j)
        return j,Ej

#实现更新Ek操作,因为除了最后我们需要更新Ei之外，我们在内循环中计算α_1与α_2时还是需要用到E1与E2，
#因为每次的E1与E2由于上一次循环中更新了α值，所以这一次也是需要更新E1与E2值，所以单独实现一个更新Ek值的方法还是有必要的
def updateEk(oS,k):
    Ek = calcEk(oS,k)
    oS.eCache[k] = [1,Ek]   #第一列1，表示为有效标识

#三：实现内循环函数，相比于外循环，这里包含了主要的更新操作
def innerL(i,oS):   #由外循环提供i值（具体选取要违背kkT<这里实现>,使用交替遍历<外循环中实现>）---提供α_1的索引
    Ei = calcEk(oS,i)   #计算E1值，主要是为了下面KKT条件需要使用到

    #如果下面违背了KKT条件，则正常进行α、Ek、b的更新,重点：后面单独说明下面是否满足违反KKT条件
    if ((oS.label[i]*Ei < -oS.toler) and (oS.alphas[i] < oS.C)) or\
       ((oS.label[i]*Ei > oS.toler) and (oS.alphas[i] > 0)):  #注意：对于硬间隔，我们直接和1对比，对于软间隔，我们要和1 +或- ε对比
        #开始在内循环中，选取差值最大的α_2下标索引
        j,Ej = selectJ(i,oS,Ei)
        #因为后面要修改α_1与α_2的值，但是后面修改阈值b的时候需要用到新旧两个值，所以我们需要在更新α值之前进行保存旧值
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()

        #分析约束条件（是对所有α都适用），一会对我们新的α_2进行截取纠正，注意：α_1是由α_2推出的，所以不需要进行验证了。
        #如果y_1!=y_2异号时：
        if oS.label[i] != oS.label[j]:
            L = max(0,alphaJold-alphaIold)
            H = min(oS.C,oS.C+alphaJold-alphaIold)
        else:   #如果y_1==y_2同号时
            L = max(0,alphaJold+alphaIold-oS.C)
            H = min(oS.C,alphaJold+alphaIold)
        #上面就是将α_j调整到L,H之间
        if L==H:    #如果L==H，之间返回0，跳出这次循环，不进行改变(单值选择，没必要）
            return 0

        #计算η值=k_11+k_22-2k_12
        eta = oS.K[i,i] + oS.K[j,j] - 2.0*oS.K[i,j]  #eta性质可以知道是>=0的，所以我们只需要判断是否为0即可
        if eta <= 0:
            print("eta <= 0")
            return 0

        #当上面所有条件都满足以后，我们开始正式修改α_2值，并更新对应的Ek值
        oS.alphas[j] += oS.label[j]*(Ei-Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        updateEk(oS,j)

        #查看α_2是否有足够的变化量，如果没有足够变化量，我们直接返回，不进行下面更新α_1,注意：因为α_2变化量较小，所以我们没有必要非得把值变回原来的旧值
        if abs(oS.alphas[j] - alphaJold) < 0.00001:
            print("J not move enough")
            return 0

        #开始更新α_1值,和Ek值
        oS.alphas[i] += oS.label[i]*oS.label[j]*(alphaJold-oS.alphas[j])
        updateEk(oS,i)

        #开始更新阈值b,正好使用到了上面更新的Ek值
        b1 = oS.b - Ei - oS.label[i] * (oS.alphas[i] - alphaIold) * oS.K[i,i] - oS.label[j] * (
                    oS.alphas[j] - alphaJold) * oS.K[i,j]

        b2 = oS.b - Ej - oS.label[i] * (oS.alphas[i] - alphaIold) * oS.K[i,j] - oS.label[j] * (
                    oS.alphas[j] - alphaJold) * oS.K[j,j]

        #根据统计学习方法中阈值b在每一步中都会进行更新，
        #1.当新值alpha_1不在界上时(0<alpha_1<C)，b_new的计算规则为：b_new=b1
        #2.当新值alpha_2不在界上时(0 < alpha_2 < C)，b_new的计算规则为：b_new = b2
        #3.否则当alpha_1和alpha_2都不在界上时，b_new = 1/2(b1+b2)
        if oS.alphas[i] > 0 and oS.alphas[i] < oS.C:
            oS.b = b1
        elif oS.alphas[j] > 0 and oS.alphas[j] < oS.C:
            oS.b = b2
        else:
            oS.b = 1/2*(b1+b2)

        #注意：这里我们应该根据b_new更新一次Ei，但是我们这里没有写，因为我们将这一步提前到了最开始，即selectJ中

        #以上全部更新完毕，开始返回标识
        return 1
    return 0    #没有违背KKT条件

#四：开始外循环，由于我们在内循环中实现了KKT条件的判断，所以这里我们只需要进行交替遍历即可
#交替遍历一种方式是在所有的数据集上进行单遍扫描，另一种是在非边界上（不在边界0或C上的值）进行单遍扫描
# 交替遍历：
# 交替是通过一个外循环来选择第一个alpha值的，并且其选择过程会在两种方式之间交替：
# 一种方式是在所有数据集上进行单遍扫描，
# 另一种方式则是在非边界alpha中实现单遍扫描，所谓非边界alpha指的是那些不等于边界0或C的alpha值。
# 对整个数据集的扫描相当容易，
# 而实现非边界alpha值的扫描时，首先需要建立这些alpha值的列表，然后对这个表进行遍历。
# 同时，该步骤会跳过那些已知不变的alpha值。
def smoP(data_X,data_Y,C,toler,maxIter,kTup):
    oS = optStruct(data_X,data_Y,C,toler,kTup)
    iter = 0
    entireSet = True    #标志是否应该遍历整个数据集
    alphaPairsChanged = 0   #标志一次循环中α更新的次数
    #开始进行迭代
    #当iter >= maxIter或者((alphaPairsChanged == 0) and not entireSet)退出循环
    #前半个判断条件很好理解，后面的判断条件中，表示上一次循环中，是在整个数据集中遍历，并且没有α值更新过，则退出
    while iter < maxIter and ((alphaPairsChanged > 0) or entireSet):
        alphaPairsChanged = 0
        if entireSet:   #entireSet是true，则在整个数据集上进行遍历
            for i in range(oS.m):
                alphaPairsChanged += innerL(i,oS)   #调用内循环
            print("full dataset, iter: %d i:%d,pairs changed:%d"%(iter,i,alphaPairsChanged))
            iter += 1   #无论是否更新过，我们都计算迭代一次
        else:   #遍历非边界值
            nonBounds = np.where((oS.alphas>0) & (oS.alphas<C))[0]  #获取非边界值中的索引
            for i in nonBounds: #开始遍历
                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还是False,下一次还是在非边界上进行遍历。可以认为这里是倾向于非边界遍历，因为非边界遍历的样本更符合内循环中的违反KKT条件
            entireSet = True

        print("iteration number: %d"%iter)

    return oS.b,oS.alphas

def calcWs(alphas,data_X,data_Y):
    #根据西瓜书6.37求W
    m,n = data_X.shape
    w = np.zeros(n)
    for i in range(m):
        w += alphas[i]*data_Y[i]*data_X[i].T

    return w

#绘制图像
def plotFigure(data_X,data_Y,alphas):
    m,n = data_X.shape
    # 进行数据集分类操作
    cls_1x = data_X[np.where(data_Y==1)]
    cls_2x = data_X[np.where(data_Y!=1)]

    plt.scatter(cls_1x[:,0].flatten(), cls_1x[:,1].flatten(), s=30, c='b')
    plt.scatter(cls_2x[:,0].flatten(), cls_2x[:,1].flatten(), s=30, c='g', marker='s')

    # 画出支持向量点
    for i in range(m):
        if alphas[i] > 0.0:
            plt.scatter(data_X[i, 0], data_X[i, 1], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='r')

    plt.xlim((-1.5, 1.5))
    plt.ylim((-1, 1.5))
    plt.show()

def testRbf(data_X,data_Y,C,toler,maxIter,kTup):
    #获取阈值b和α值
    b,alphas = smoP(data_X,data_Y,C,toler,maxIter,kTup)
    plotFigure(data_X,data_Y,alphas)
    # print(b)
    # print(alphas)
    #获取支持向量，即α>0点
    svIdx = np.where(alphas>0)[0]
    sVs = data_X[svIdx]
    labelSV = data_Y[svIdx]
    print("there are %d support vectors"%sVs.shape[0])
    m,n = data_X.shape
    errorCount = 0

    #获取数据异常错误率
    for i in range(m):
        #先进行预测，预测公式见西瓜书6.24
        kernelEval = kernelTrans(sVs,data_X[i],kTup)    #使用支持向量点，和对应的预测数据点，进行预测
        pred = kernelEval.T@np.multiply(labelSV,alphas[svIdx]) + b #进行预测
        if np.sign(pred) != np.sign(data_Y[i]):
            errorCount += 1
    print("the trainning set error rate is: %f"%(errorCount/m))

    #测试集预测
    data_X,data_Y = loadDataSet("testSetRBF2.txt")
    m,n = data_X.shape
    errorCount = 0

    #获取数据异常错误率
    for i in range(m):
        #先进行预测，预测公式见西瓜书6.24
        kernelEval = kernelTrans(sVs,data_X[i],kTup)    #使用支持向量点，和对应的预测数据点，进行预测
        pred = kernelEval.T@np.multiply(labelSV,alphas[svIdx]) + b #进行预测
        if np.sign(pred) != np.sign(data_Y[i]):
            errorCount += 1
    print("the testing set error rate is: %f"%(errorCount/m))

data_X,data_Y = loadDataSet("testSetRBF.txt")
C = 200
toler = 0.00001
maxIter = 10000
kTup = ('rbf',1.3)

testRbf(data_X,data_Y,C,toler,maxIter,kTup)

（六）结果显示

kTup=('rbf',1.3)时：

 

kTup=('rbf',0.1)时：

 

五：手写数字识别问题

（一）与KNN对比

与KNN算法对比，KNN算法效果不错，但是想要保留所有的训练样本，在预测的时候也是需要所有的样本，但是我们使用SVM时，我们只需要保留支持向量即可，保留的样本少了很多，效果不会太差。

注意:我们主要实现了SVM对数字识别的二分类问题,主要对数字1和9进行了分类

（二）代码实现获取数据集

#将每一个数字文件转换为矩阵向量
def img2vector(filename):
    data = []
    with codecs.open(filename,'r') as fp:
        for i in range(32):
            linestr = fp.readline() #读取一行数据
            for j in range(32):
                data.append(int(linestr[j])) #添加数据
        fp.close()
    return np.array(data)

def loadImages(path):
    # 读取数据
    hwLabels = []
    filelist = listdir(path)  # 获取所有文件目录
    m = len(filelist)
    data = np.zeros((m, 1024))

    # 先获取标签值
    for i in range(m):
        filename = filelist[i]
        classNumStr = int(filename.split('_')[0])
        if classNumStr == 9: hwLabels.append(-1)
        else: hwLabels.append(1)
        data[i, :] = img2vector("%s/%s" % (path, filename))

    return data, hwLabels

（三）测试数据集错误率

def testDights(data_X,data_Y,C,toler,maxIter,kTup):
    #获取阈值b和α值
    b,alphas = smoP(data_X,data_Y,C,toler,maxIter,kTup)
    #获取支持向量，即α>0点
    svIdx = np.where(alphas>0)[0]
    sVs = data_X[svIdx]
    labelSV = data_Y[svIdx]
    print("there are %d support vectors"%sVs.shape[0])
    m,n = data_X.shape
    errorCount = 0

    #获取数据异常错误率
    for i in range(m):
        #先进行预测，预测公式见西瓜书6.24
        kernelEval = kernelTrans(sVs,data_X[i],kTup)    #使用支持向量点，和对应的预测数据点，进行预测
        pred = kernelEval.T@np.multiply(labelSV,alphas[svIdx]) + b #进行预测
        if np.sign(pred) != np.sign(data_Y[i]):
            errorCount += 1
    print("the trainning set error rate is: %f"%(errorCount/m))

    #测试集预测
    data_X, data_Y = loadImages("digits/testDigits")
    m,n = data_X.shape
    errorCount = 0

    #获取数据异常错误率
    for i in range(m):
        #先进行预测，预测公式见西瓜书6.24
        kernelEval = kernelTrans(sVs,data_X[i],kTup)    #使用支持向量点，和对应的预测数据点，进行预测
        pred = kernelEval.T@np.multiply(labelSV,alphas[svIdx]) + b #进行预测
        if np.sign(pred) != np.sign(data_Y[i]):
            errorCount += 1
    print("the testing set error rate is: %f"%(errorCount/m))

（四）进行测试

data_X,data_Y = loadImages("digits/trainingDigits")
data_Y = np.array(data_Y)   #注意我们loadImages返回的data_Y是列表，我们要转化为numpy数组类型
C = 200
toler = 0.0001
maxIter = 10000
kTup = ('rbf',10)

testDights(data_X,data_Y,C,toler,maxIter,kTup)