【机器学习算法原理及其实现剖析】SVM篇

往往我们只是知道了简单的证明，应付了考试，但是却没有把所学的东西，真正的去实践和应用一番，又可能苦于代码或理论难理解，而选择直接放弃。

本主题【机器学习算法原理及其实现剖析】就是，把常用的机器学习算法的原理实现和对应理论支撑相结合，希望能提供一种不同的学习思路。

SVM

如果想明白一个东西，第一步就是知道它用来干什么的？也就是解决了什么问题。

1、解决什么问题？

最基本的应用是数据分类，特别是对于非线性不可分数据集。支持向量机不仅能对非线性可分数据集进行分类，对于非线性不可分数据集的也可以分类（这也是SVM应用广泛的重要原因）

现实场景一 ：样本数据大部分是线性可分的，但是只是在样本中含有少量噪声或特异点，去掉这些噪声或特异点后线性可分 => 用支持向量机的软间隔方法进行分类；

现实场景二：样本数据完全线性不可分 => 引入核函数， 将低维不可分的非线性数据集转化为高维可分的数据集，用支持向量机的软间隔方法进行分类；

对分类器准确性有影响只有样本中的支持向量，支持向量机的名字由来也正是如此。
所以，问题就转化为获取离边界最近样本点（支持向量）到超平面最远距离。在代码中就是一个entireSet标志，地磁全部遍历，之后就变遍历支持向量，此处以SMO代码举例。

    # 全量遍历标志
    entireSet = True
    # alpha对是否优化标志
    alphaPairsChanged = 0
    # 外循环 终止条件：1.达到最大次数 或者 2.alpha对没有优化
    while (iter < maxInter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        # 全量遍历 ，遍历每一行数据 alpha对有修改，alphaPairsChanged累加
        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:
            # 获取(0，C)范围内数据索引列表，也就是只遍历属于支持向量的数据
            nonBounds = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]  # .A将mat转换为array
            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 = True
        print("iteation number: %d" % iter)
        print("entireSet :%s" % entireSet)
        print("alphaPairsChanged :%d" % alphaPairsChanged)
    return oS.b, oS.alphas```

1.1 软间隔问题

但是这样就有一个问题（软间隔问题），如果支持向量是离群点的时候，就会影响最终结果，于是引入松弛变量。以hinge损失为例，

这个公式不好懂，没关系，我们可以可视化成下图这样，

因此，最终的优化目标函数变成了两个部分组成：距离函数和松弛变量误差
上式中，我们加入权重参数C，将其与目标函数中的松弛变量误差相乘，这样，就可以通过调整C来对二者的系数进行调和。

简单概括一下规律
（1) $\alpha$ = C， 松弛系数在0-1之间，是分类没有合并的
（2）松弛系数为0，表示分类正确
（3）$\alpha$ = C ，松弛系数大于1 表示分类错误。

1.2 非线性不可分问题

思想：对于在N维空间中线性不可分的数据，在N+1维以上的空间会有更大到可能变成线性可分的。
这样，我们将原问题变成了如何对原始数据进行映射，才能使其在新空间中线性可分。
为了实现这种映射，引入了核函数的方法。
那为什么核函数一定起作用呢？这里给出一段引用：

RBF核函数可以将维度扩展到无穷维的空间，因此，理论上讲可以满足一切映射的需求。为什么会是无穷维呢？我以前都不太明白这一点。后来老师讲到，RBF对应的是泰勒级数展开，在泰勒级数中，一个函数可以分解为无穷多个项的加和，其中，每一个项可以看做是对应的一个维度，这样，原函数就可以看做是映射到了无穷维的空间中。这样，在实际应用中，RBF是相对最好的一个选择。当然，如果有研究的话，还可以选用其他核函数，可能会在某些问题上表现更好。但是，RBF是在对问题不了解的情况下，对最广泛问题效果都很不错的核函数。因此，使用范围也最广。

# 核转换函数（一个特征空间映射到另一个特征空间，低维空间映射到高维空间）
# 高维空间解决线性问题，低维空间解决非线性问题
# 线性内核 = 原始数据矩阵（100*2）与原始数据第一行矩阵转秩乘积（2*1） =>（100*1）
# 非线性内核公式：k(x,y) = exp(-||x - y||**2/2*(e**2))
# 1.原始数据每一行与原始数据第一行作差，
# 2.平方
'''
X: dataMat 原始特征数据
Y: rowDataMat 标签数据
kTup 指定核方式
'''
def kernelTrans(dataMat, rowDataMat, kTup):
    m , n =np.shape(dataMat)
    # 初始化核矩阵 m*1
    K = np.mat(np.zeros((m ,1)))
    # 线性核
    if kTup[0] == 'lin':
        K = dataMat * rowDataMat.T
    # 高斯核
    elif kTup[0] == 'rbf':  # 非线性核
        for j in range(m):
            # xi - xj
            deltaRow = dataMat[j ,:] - rowDataMat
            K[j] = deltaRow *deltaRow.T
        # 1*m m*1 => 1*1
        K = np.exp( K /(- 2 *kTup[1]**2))  # m个样本1列

    # 若需要不同核，下面可以继续添加

    else: raise NameError('Houston We Have a Problem -- That Kernel is not recognized')
    return K

至于拉格朗日对偶法转化目标函数，这里就不再做解释，直接给出优化目标。

2、实现

大致明白SVM解决什么问题之后，我们以其中一种简单实现SMO来看理论与实现的结合。
一句话，由于拉格朗日变换后的$\alpha$比较多,SMO算法则采用了一种启发式的方法。它每次只优化两个变量，将其他的变量都视为常数。
换句话说，由于${\sum }_{i=1}^m{\alpha }_i{y}_i=0$.假如将${\alpha }_3$，${\alpha }_4$，${\alpha }_m$　固定，那么${\alpha }_1$,${\alpha }_2$之间的关系也确定了。这样SMO算法将一个复杂的优化算法转化为一个比较简单的两变量优化问题。

下面是重头戏：
首先，我们在上面目标优化式子（下图）基础上，单独考虑如何计算变量${\alpha }_2$

这里我们就需要引入误差和迭代的思想，首先引入对于变量${\alpha }_2$误差的计算：

相信上面一部很好理解，最后就是利用迭代的思想更新${\alpha }_2$：

更新$\alpha$的代码如下：

# 第一次随机选取不等于i的数据项，其后根据误差最大选取数据项
        j, Ej = selectJ(i, oS, Ei)   # 这一步运用了启发式
        # 初始化，开辟新的内存，为了不将原来内存占用而改变值
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()
        # 通过 a1y1 + a2y2 = 常量
        #    0 <= a1,a2 <= C 求出L,H
        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
        # 内核分母 K11 + K22 - 2K12
        eta = oS.K[i, i] + oS.K[j, j] - 2.0 * oS.K[i, j]
        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], L, H)
        # alpha有改变，就需要更新缓存数据
        updateEk(oS, j)  # 主要更改self.eCache值
        # 如果优化后的alpha 与之前的alpha变化很小，则舍弃，并重新选择数据项的alpha
        if (abs(oS.alphas[j] - alphaJold) < 0.00001):   # 判断是否改变较小，若太小，则不改变，因改变太小没有意义
            print("j not moving enough, abandon it.")
            return 0
        # 计算alpha对的另一个alpha i
        # ai_new*yi + aj_new*yj = 常量
        # ai_old*yi + ai_old*yj = 常量
        # 作差=> ai = ai_old + yi*yj*(aj_old - aj_new)
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j])
        # alpha有改变，就需要更新缓存数据
        updateEk(oS, i) # 主要更改self.eCache值

这里简单说下，两个变量的选择：

最后，就是更新bias和误差便于下次迭代的变量选择：

        bi = 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]
        bj = 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]
        # 首选alpha i，相对alpha j 更准确
        if (0 < oS.alphas[i]) and (oS.alphas[i] < oS.C):
            oS.b = bi
        elif (0 < oS.alphas[j]) and (oS.alphas[j] < oS.C):
            oS.b = bj
        else:
            oS.b = (bi + bj) / 2.0
        return 1
    else:
        return 0

最后感谢提供的代码：机器学习之支持向量机—SVM原理代码实现

3、小结：

贴出最终版的代码，供学习与交流之用：
算法步骤：

是不是很容易就明白了呢？~
那就点个赞加个关注吧~

import numpy as np


# 核转换函数（一个特征空间映射到另一个特征空间，低维空间映射到高维空间）
# 高维空间解决线性问题，低维空间解决非线性问题
# 线性内核 = 原始数据矩阵（100*2）与原始数据第一行矩阵转秩乘积（2*1） =>（100*1）
# 非线性内核公式：k(x,y) = exp(-||x - y||**2/2*(e**2))
# 1.原始数据每一行与原始数据第一行作差，
# 2.平方
def kernelTrans(dataMat, rowDataMat, kTup):
    m , n =np.shape(dataMat)
    # 初始化核矩阵 m*1
    K = np.mat(np.zeros((m ,1)))
    # 线性核
    if kTup[0] == 'lin':
        K = dataMat * rowDataMat.T
    # 高斯核
    elif kTup[0] == 'rbf':  # 非线性核
        for j in range(m):
            # xi - xj
            deltaRow = dataMat[j ,:] - rowDataMat
            K[j] = deltaRow *deltaRow.T
        # 1*m m*1 => 1*1
        K = np.exp( K /(- 2 *kTup[1]**2))  # m个样本1列

    # 若需要不同核，下面可以继续添加


    else: raise NameError('Houston We Have a Problem -- That Kernel is not recognized')
    return K


# 定义数据结构体，用于缓存，提高运行速度
class optStruct:
    def __init__(self, dataSet, labelSet, C, toler, kTup):
        self.dataMat = np.mat(dataSet)  # 原始数据，转换成m*n矩阵
        self.labelMat = np.mat(labelSet).T  # 标签数据 m*1矩阵
        self.C = C  # 惩罚参数，C越大，容忍噪声度小，需要优化；反之，容忍噪声度高，不需要优化；
        # 所有的拉格朗日乘子都被限制在了以C为边长的矩形里，  此处的C 就是  0<alpha(i)<C
        self.toler = toler  # 容忍度
        self.m = np.shape(self.dataMat)[0]  # 原始数据行长度
        self.alphas = np.mat(np.zeros((self.m ,1))) # alpha系数，m*1矩阵
        self.b = 0  # 偏置
        self.eCache = np.mat(np.zeros((self.m ,2))) # 保存原始数据每行的预测值
        self.K = np.mat(np.zeros((self.m ,self.m))) # 核转换矩阵 m*m
        for i in range(self.m):
            self.K[: ,i] = kernelTrans(self.dataMat, self.dataMat[i ,:], kTup)  # 核函数的转换
            # K = dataMat * rowDataMat.T



# 从文件获取特征数据，标签数据
def loadDataSet(fileName):
    dataSet = []; labelSet = []
    fr = open(fileName)
    for line in fr.readlines():
        # 分割
        lineArr = line.strip().split()
        n=len(lineArr)-1
        dataSet.append([float(lineArr[i])for i in range(n)])
        labelSet.append(float(lineArr[2]))
    m,n=np.shape(dataSet)
    print(m,n)
    dataSet=np.reshape(dataSet,(m,n))
    labelSet=np.reshape(labelSet,(1,m))
    return dataSet, labelSet



# 计算原始数据第k项对应的预测误差  1*m m*1 =>1*1
# oS：结构数据
# k： 原始数据行索引
def calEk(oS, k):
    # 用此处函数的原因可以看推导公式，将会给你详细的解释
    fXk = float(np.multiply(oS.alphas ,oS.labelMat). T *oS.K[: ,k] + oS.b)  # f(x) = w*x + b
    Ek = fXk - float(oS.labelMat[k])  # EK=fX(i)-label(i)
    return Ek


# 在alpha有改变都要更新缓存
def updateEk(oS, k):
    Ek = calEk(oS, k)
    oS.eCache[k] = [1, Ek]



# 启发式选择另一个j
def selectJ(i, oS, Ei):
    # 初始化
    maxK = -1  # 误差最大时对应索引
    maxDeltaE = 0  # 最大误差,这里类似启发式算法的比较
    Ej = 0   # j索引对应预测误差
    # 保存每一行的预测误差值 1相对于初始化为0的更改
    oS.eCache[i] = [1 ,Ei]
    # 获取数据缓存结构中非0的索引列表(先将矩阵第0列转化为数组)
    validEcacheList = np.nonzero(oS.eCache[: ,0].A)[0]   # 选择非零值，得到索引
    # 遍历索引列表，寻找最大误差对应索引
    if len(validEcacheList) > 1:
        for k in validEcacheList:
            if k == i:
                continue
            Ek = calEk(oS, k)
            deltaE = abs(Ei - Ek)
            if(deltaE > maxDeltaE):  # 找到的值大于该值就替换，以此得到最优的值
                maxK = k
                maxDeltaE = deltaE
                Ej = Ek
        return maxK, Ej
    else:
        # 随机选取一个不等于i的j
        j = i
        while (j == i):
            j = int(np.random.uniform(0, oS.m))  # 随机选取一个不等于i的索引
        Ej = calEk(oS, j)
    return j, Ej


# alpha范围剪辑
def clipAlpha(aj, L, H):
    if aj > H:
        aj = H
    if aj < L:
        aj = L
    return aj

# 计算 w 权重系数，这样就有了   y = wx + b
def calWs(alphas, dataSet, labelSet):
    '''
    m 表示样本数量，n表示一个样本的所有特征
    :param alphas: [m,1]
    :param dataSet: [m,n]
    :param labelSet:[m,1]
    :return: w
    '''
    dataMat = np.mat(dataSet)
    # 1*100 => 100*1
    labelMat = np.mat(labelSet)
    m, n = np.shape(dataMat)
    w = np.zeros((n, 1))
    alphas=np.reshape(alphas, (m, 1))
    for i in range(m):
        w += np.multiply(alphas[i] * labelMat[i], dataMat[i,:].T)  # 可以查看推导公式
    return w


# 计算原始数据每一行alpha,b，保存到数据结构中，有变化及时更新
def innerL(i, oS):
    # 该函数每次通过改变类中的self.alphas和更新echache值，完成alphas的更新和对应的EI更新
    # 计算预测误差
    Ei = calEk(oS, i)
    # 选择第一个alpha，违背KKT条件2
    # 正间隔，负间隔
    if ((oS.labelMat[i] * Ei < -oS.toler) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i] * Ei > oS.toler) and (oS.alphas[i] > 0)):
        # 第一次随机选取不等于i的数据项，其后根据误差最大选取数据项
        j, Ej = selectJ(i, oS, Ei)   # 这一步运用了启发式
        # 初始化，开辟新的内存，为了不将原来内存占用而改变值
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()
        # 通过 a1y1 + a2y2 = 常量
        #    0 <= a1,a2 <= C 求出L,H
        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
        # 内核分母 K11 + K22 - 2K12
        eta = oS.K[i, i] + oS.K[j, j] - 2.0 * oS.K[i, j]
        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], L, H)
        # alpha有改变，就需要更新缓存数据
        updateEk(oS, j)  # 主要更改self.eCache值
        # 如果优化后的alpha 与之前的alpha变化很小，则舍弃，并重新选择数据项的alpha
        if (abs(oS.alphas[j] - alphaJold) < 0.00001):   # 判断是否改变较小，若太小，则不改变，因改变太小没有意义
            print("j not moving enough, abandon it.")
            return 0
        # 计算alpha对的另一个alpha i
        # ai_new*yi + aj_new*yj = 常量
        # ai_old*yi + ai_old*yj = 常量
        # 作差=> ai = ai_old + yi*yj*(aj_old - aj_new)
        oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j])
        # alpha有改变，就需要更新缓存数据
        updateEk(oS, i) # 主要更改self.eCache值
        ''' 
        计算b1,b2
        y(x) = w*x + b => b = y(x) - w*x
        w = aiyixi(i= 1->N求和)
        b1_new = y1_new - (a1_new*y1*k11 + a2_new*y2*k21 + ai*yi*ki1(i = 3 ->N求和 常量))
        b1_old = y1_old - (a1_old*y1*k11 + a2_old*y2*k21 + ai*yi*ki1 (i = 3 ->N求和 常量))
        作差=> b1_new = b1_old + (y1_new - y1_old) - y1*k11*(a1_new - a1_old) - y2*k21*(a2_new - a2_old)
        => b1_new = b1_old + Ei - yi*(ai_new - ai_old)*kii - yj*(aj_new - aj_old)*kij
        同样可推得 b2_new = b2_old + Ej - yi*(ai_new - ai_old)*kij - yj*(aj_new - aj_old)*kjj
        kij实际为核函数，xi*xj相当于将1行n列与转置n行1列相乘，得到一个数  
        
        为什么Ei=y1_new-y1_old，上面不是公式Ei=yi-label[i]吗？
        答：试想每次更改得到的y是与标签作比较，为此y1_old就可看成label[i]
        '''
        bi = 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]
        bj = 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]
        # 首选alpha i，相对alpha j 更准确
        if (0 < oS.alphas[i]) and (oS.alphas[i] < oS.C):
            oS.b = bi
        elif (0 < oS.alphas[j]) and (oS.alphas[j] < oS.C):
            oS.b = bj
        else:
            oS.b = (bi + bj) / 2.0
        return 1
    else:
        return 0


# 完整SMO核心算法，包含线性核核非线性核，返回alpha,b
# dataSet 原始特征数据
# labelSet 标签数据
# C 凸二次规划参数
# toler 容忍度
# maxInter 循环次数
# kTup 指定核方式
# 程序逻辑：
# 第一次全部遍历，遍历后根据alpha对是否有修改判断，
# 如果alpha对没有修改，外循环终止；如果alpha对有修改，则继续遍历属于支持向量的数据。
# 直至外循环次数达到maxIter
# 相比简单SMO算法，运行速度更快，原因是：
# 1.不是每一次都全量遍历原始数据，第一次遍历原始数据，
# 如果alpha有优化，就遍历支持向量数据，直至alpha没有优化，然后再转全量遍历，这是如果alpha没有优化，循环结束；
# 2.外循环不需要达到maxInter次数就终止；

def smoP(dataSet, labelSet, C, toler, maxInter, kTup=('lin', 0)):
    '''
    :param dataSet:  输入dataX
    :param labelSet: 输入label
    :param C:  限制0<α<C
    :param toler: yi*Ei的限制
    :param maxInter: 最大的迭代数
    :param kTup: 选用什么核
    :return: oS.b, oS.alphas
    '''
    oS = optStruct(dataSet, labelSet, C, toler, kTup)  # 初始化结构体类，获取实例
    iter = 0
    # 全量遍历标志
    entireSet = True
    # alpha对是否优化标志
    alphaPairsChanged = 0
    # 外循环 终止条件：1.达到最大次数 或者 2.alpha对没有优化
    while (iter < maxInter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        # 全量遍历 ，遍历每一行数据 alpha对有修改，alphaPairsChanged累加
        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:
            # 获取(0，C)范围内数据索引列表，也就是只遍历属于支持向量的数据
            nonBounds = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]  # .A将mat转换为array
            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 = True
        print("iteation number: %d" % iter)
        print("entireSet :%s" % entireSet)
        print("alphaPairsChanged :%d" % alphaPairsChanged)
    return oS.b, oS.alphas





# 绘制支持向量
def drawDataMap(dataArr, labelArr, b, alphas):
    import matplotlib.pyplot as plt
    # alphas.A>0 获取大于0的索引列表，只有>0的alpha才对分类起作用
    svInd = np.nonzero(alphas.A > 0)[0]
    # 分类数据点
    classified_pts = {'+1': [], '-1': []}       # 记录分类结果
    m, n = np.shape(dataArr)
    labelArr=np.reshape(labelArr,(m,1))

    for point, label in zip(dataArr, labelArr):
        if label == 1.0:
            classified_pts['+1'].append(point)
        else:
            classified_pts['-1'].append(point)


    fig = plt.figure()
    ax = fig.add_subplot(111)
    # 绘制数据点
    for label, pts in classified_pts.items():
        pts = np.array(pts)
        ax.scatter(pts[:, 0], pts[:, 1], label=label)
    # 绘制分割线
    w = calWs(alphas, dataArr, labelArr)
    # 函数形式：max( x ,key=lambda a : b )        #    x可以是任何数值，可以有多个x值
    # 先把x值带入lambda函数转换成b值，然后再将b值进行比较
    x1, _ = max(dataArr, key=lambda x: x[0])  # 第一维度最大值
    x2, _ = min(dataArr, key=lambda x: x[0])  # 第一维度最小值
    a1, a2 = w
    y1, y2 = (-b - a1 * x1) / a2, (-b - a1 * x2) / a2
    # 矩阵转化为数组.A
    ax.plot([x1, x2], [y1.A[0][0], y2.A[0][0]])

    # 绘制支持向量
    for i in svInd:
        x, y = dataArr[i]
        ax.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='#AB3319')
    plt.show()

    # alpha>0对应的数据才是支持向量，过滤不是支持向量的数据
    sVs = np.mat(dataArr)[svInd]  # get matrix of only support vectors
    print("there are %d Support Vectors.\n" % np.shape(sVs)[0])





# 训练结果
def getTrainingDataResult(dataSet, labelSet, b, alphas, k1=1.3):
    datMat = np.mat(dataSet)
    # 100*1
    labelMat = np.mat(labelSet).T
    # alphas.A>0 获取大于0的索引列表，只有>0的alpha才对分类起作用
    svInd = np.nonzero(alphas.A > 0)[0]
    sVs = datMat[svInd]
    labelSV = labelMat[svInd]

    m, n = np.shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], ('rbf', k1))
        # y(x) = w*x + b => b = y(x) - w*x
        # w = aiyixi          (i= 1->N求和)
        predict = kernelEval.T * np.multiply(labelSV, alphas[svInd]) + b

        if np.sign(predict) != np.sign(labelSet[i]): errorCount += 1
    print("the training error rate is: %f" % (float(errorCount) / m))



def getTestDataResult(dataSet, labelSet, b, alphas, k1=1.3):
    datMat = np.mat(dataSet)
    # 100*1
    labelMat = np.mat(labelSet).T
    # alphas.A>0 获取大于0的索引列表，只有>0的alpha才对分类起作用
    svInd = np.nonzero(alphas.A > 0)[0]
    sVs = datMat[svInd]
    labelSV = labelMat[svInd]
    m, n = np.shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs, datMat[i, :], ('rbf', k1))
        # y(x) = w*x + b => b = y(x) - w*x
        # w = aiyixi(i= 1->N求和)
        predict = kernelEval.T * np.multiply(labelSV, alphas[svInd]) + b
        if np.sign(predict) != np.sign(labelSet[i]): errorCount += 1
    print("the test error rate is: %f" % (float(errorCount) / m))




if __name__=='__main__':
    path='D:\\CODE\\SVM\\SVM.txt'
    dataSet, labelSet=loadDataSet(path)
    print(dataSet.shape,labelSet.shape)

    # b,alphas=smoP(dataSet, labelSet, C, toler, maxInter, kTup=('lin', 0) )
    b, alphas = smoP(dataSet, labelSet, 4, 0.5, 50, kTup=('rbf', 0.05))
    print('b=',b)
    print('alphas=',alphas)
    drawDataMap(dataSet, labelSet, b, alphas)