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

最新推荐文章于 2022-10-30 11:24:58 发布

weixin_44275584

最新推荐文章于 2022-10-30 11:24:58 发布

阅读量168

点赞数

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

本文链接：https://blog.csdn.net/weixin_44275584/article/details/110499135

版权

机器学习专栏收录该内容

5 篇文章 1 订阅

订阅专栏

支持向量机就是希望找到到支持向量距离最小的那一个超平面。

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

线性可分：可以使用一条线（超平面）完全把‘0’，‘1’两种类型的数据分割开——该线称为分割超平面，也就是分类的决策边界
支持向量机就是希望能采取数据点离决策边界最远的方法来构建分类器。

分类器工作原理：
输入数据给分类器，会输出一个分类标签，这相当于类似sigmoid函数的作用。

超平面确定后，各个点到超平面的距离称为间隔。距离分割超平面最近的几个点我们称之为支持向量（support vector）。分割超平面可能会有很多，需要取间隔最大的那一个，也就是使得支持向量的间隔最大。
分割超平面用公式可以表示为:
$\vec{w}^{T} \vec{x} + b$
每个点到超平面的距离可以表示为：
$\frac{\left | \vec{w}^{T} \vec{x} + b \right |}{\begin{Vmatrix} \vec{w} \end{Vmatrix}} = \frac{ y *(\vec{w}^{T} \vec{x} + b) }{\begin{Vmatrix} \vec{w} \end{Vmatrix}}$

我们要找的是所有点中距离超平面最远的，用公式表示为：
$\underset{n}{min}(d)= \underset{n}{min}(\frac{\left | \vec{w}^{T} \vec{x} + b \right |}{\begin{Vmatrix} \vec{w} \end{Vmatrix}}) = \underset{n}{min}( \frac{ y * (\vec{w}^{T} \vec{x} + b) }{\begin{Vmatrix} \vec{w} \end{Vmatrix}} )$

接下来的问题在于找到这样的分类器，就需要找到能够使得上述最小距离最大化的 $\vec{w}$ 和 $b$ 所确定的那个超平面:
$\underset{ \vec{w},b}{max} \{\underset{n}{min}(d)\}= \underset{ \vec{w},b}{max} \{\underset{n}{min}(\frac{\left | \vec{w}^{T} \vec{x} + b \right |}{\begin{Vmatrix} \vec{w} \end{Vmatrix}}) \} =\underset{ \vec{w},b}{max} \{\underset{n}{min}(\frac{ y *(\vec{w}^{T} \vec{x} + b) }{\begin{Vmatrix} \vec{w} \end{Vmatrix}}) \}$

其中 $y$ 为类别标签，取值为-1 和1

上述式子的求解非常复杂，为了更容易求解，令支持向量的 $*(\vec{w}^{T} \vec{x} + b)$ 都为1，那么就可以通过求 $\begin{Vmatrix} \vec{w} \end{Vmatrix}^{-1}$ 的最大值来得到最终的解，但并非所有数据点的 $*(\vec{w}^{T} \vec{x} + b)$ 都为1，因此多了一个约束条件 $*(\vec{w}^{T} \vec{x} + b) ≥1$ ，对于此类问题，可以使用拉格朗日乘子法求解，最终优化目标函数为：
$max[\sum_{i=1}^{m}\alpha -\frac{1}{2}\sum_{i,j=1}^{m}y_{i} \cdot y_{j} \cdot a_{i} \cdot a_{j} <x_{i},x_{j}> ]$
其中 $x_{i},x_{j}>$ 为两个向量的内积。
约束条件为：
$\alpha \geq 0 , \sum_{i=1}^{m}\alpha_{i} \cdot y_{i} = 0$
至此，我们解决了线性可分下的分割，对于线性不可分的情况，我们可以引入松弛变量（slack variable），约束条件变为：
$\geq \alpha \geq 0 , \sum_{i=1}^{m}\alpha_{i} \cdot y_{i} = 0$
这里的常数C用于控制“最大化间隔”和“保证大部分点的函数间隔小于1.0”这两个目标的权重。

SMO算法

SMO表示序列最小优化（Sequential minimal optimiza），该算法将大优化问题分解为多个小优化问题来求解。这些小优化问题往往很容易求解，并且对他们进行顺序求解的结果与将他们整体求解得到的结果一致。
SOM算法的目标是求出一系列 $\alpha$ 和b，一旦求出这些 $\alpha$ ，就很容易计算出权重向量并得到分割超平面
SMO算法的工作原理：每次循环中选择两个 $\alpha$ 进行优化处理。一点找到一对合适的 $\alpha$ ，那么就增大其中一个同时减小另一个。这里的合适指要满足两个条件，条件一是这两个 $\alpha$ 必须要在间隔边界之外，第二个条件是这两个 $\alpha$ 还没有进行过区间化处理或者不在边界上。

简化版SMO算法

SMO算法的完整实现需要大量代码。下述代码会对算法进行简化处理，以便了解算法的基本工作思路。简化版代码虽然少，但是执行速度慢。简化版会跳过在外循环中确定需要优化的最佳alpha对这一步骤，首先在数据及上遍历每一个alpha，然后在剩下的alpha集合中随机选择另一个alpha，从而建立alpha对。在这里我们要改变两个alpha，之所以要这样，是因为我们有一个约束条件：
$\sum_{i=1}^{m}\alpha_{i} \cdot y_{i} = 0$
由于只改变一个alpha可能会导致该约束条件失效，因此需要同时改变两个alpha

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

测试数据

import random
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
data = []
for i in range(50):
    x1 = random.randint(20,50)
    x2 = random.uniform(2,3)*x1
    if x1 > x2 :
        z = 1
    else:
        z = -1
    data.append([x1,x2,z])
for i in range(50):
    x1 = random.randint(20,50)
    x2 = random.uniform(0.5,0.8)*x1
    if x1 > x2 :
        z = 1
    else:
        z = -1
    data.append([x1,x2,z])
   
    
data = np.array(data)  
plt.scatter(data[:,0],data[:,1]) 
plt.show()

在这里插入图片描述

def loadDataSet(fileName):
    '''
    读取文件数据，得到数据集和标签
    '''
    dateMat = []
    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):
    import random
    '''
    i: 第一个alpha的下标
    m: 所有alpha的数目
    在0-m（不含）内随机进行选择不为i的随机数
    '''
    j = i
    while (j==i):
        j = random.randint(0,m-1)
    return j

def clipAlpha(aj, H, L):
    '''
    使得aj（alpha）在L到H 的范围内
    '''
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    '''
    dataMatIn:数据集
    classLabels ： 类别标签
    C： 松弛变量-常量
    toler：容错率
    maxIter：最大循环次数
    '''
    dataMatrix = np.mat(dataMatIn)  # 转换为矩阵
    labelMat = np.mat(classLabels).transpose()   #转换为矩阵并转置，得到列向量
    b = 0
    m, n = np.shape(dataMatrix)
    alphas = np.mat(np.zeros((m,1))) #列向量
    iter = 0  # 存放在没有alpha改变时，遍历数据集的次数；当iter>=maxIter时，循环结束
    while (iter < maxIter):
        alphaPairsChanged = 0   # 用于记录alpha是否已经优化
        for i in range(m):
            fXi = float(np.multiply(alphas,labelMat).T * (dataMatrix * dataMatrix[i,:].T)) + b  #预测类别 np.multiply‘对应位置’相乘
            Ei = fXi - float(labelMat[i])  #误差，若果误差大，则优化alpha
            if ((labelMat[i] * Ei < -toler) and (alphas[i] < C)) or ((labelMat[i] * Ei > toler) and (alphas[i] > 0)):
                j = selectJrand(i,m)   # 任意选取另一个alpha
                fXj = float(np.multiply(alphas,labelMat).T * (dataMatrix * dataMatrix[j,:].T)) + b
                Ej = fXj - float(labelMat[j])
                alphaIold= alphas[i].copy()  #存旧的alpha
                alphaJold= alphas[j].copy()
                # 用于将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])
                if L==H :
                    print("L==H")
                    continue 
                # 最优修改量，如果>=0，则退出当前迭代
                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)    # 调整至L-H之间
                if (abs(alphas[j] - alphaJold) < 0.00001) :
                    print("j not moving enough ")
                    continue
                alphas[i] += labelMat[j] * labelMat[i] * (alphaJold - alphas[j])
                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
                print("iter:%d i:%d , paris changed %d" %(iter, i, alphaPairsChanged))
        if (alphaPairsChanged == 0):
            iter += 1
        else :
            iter = 0
        print("iteration number : %d" %iter)
    return b, alphas

# 调用
dataArray = data[:,0:2]
labelArray = data[:,2]

b,alphas = smoSimple(dataArray,labelArray,0.6,0.001,40)

结果

b = 0.91825168
alphas[alphas>0] = [3.07209585e-03, 5.75965503e-05, 2.36928839e-03, 6.36523659e-04]

for i in range(100):
    if alphas[i] > 0.0:
        print(dataArray[i],labelArray[i])
# 支持向量
[24.         51.69105293] -1.0
[30.         23.44589062] 1.0
[46.         36.67181744] 1.0
[20.         15.62208869] 1.0
 
# 画图
plt.scatter(data[:50,0],data[:50,1],c = 'pink') 
plt.scatter(data[50:,0],data[50:,1],c = 'green') 

x = []
for i in range(100):
    if alphas[i] > 0.0:
        print(dataArray[i],labelArray[i])
        x.append(dataArray[i])
x = np.array(x)
# 画散点图
plt.scatter(x[:,0], x[:,1],c ='black' , s=20, alpha=0.5)


plt.plot()
plt.show()

在这里插入图片描述

利用完整 Platt SMO算法加速优化

简化版的SMO算法在计算小规模数据集上没有问题，但是在大型数据集上运行速度就会变慢。完整版和简化版之间，在实现alpha的更改和代数优化环节是一致的，唯一不同之处在于选择alpha的方式。
Platt SMO 算法是通过一个外循环来选择第一个alpha值的，并且选择过程会在两种方法之间进行交替：1.在所有数据集上进行单遍扫描；2.在非边界alpha中实现单遍扫描（非边界：不等于边界0或C的alpha值）。该步骤会跳过哪些已知的不会改变的alpha值。
在选择第一个alpha值后，算法会通过一个内循环来选择第二个alpha值。在优化过程中，会通过最大化步长的方式来获取第二个alpha值。在简化版SMO算法中，会在选择j之后计算错误率Ej。而在完整版中，会建立一个全局的缓存用于保存误差值，并从中选择使得步长或者Ei-Ej最大的alpha值。

import numpy as np
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)))  # m*2 误差缓存。第一列：是否有效；第二列：实际的E值
    
def calcEk(oS,k):
    '''
    计算E值并返回
    '''
    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 selectJ(i, oS, Ei):
    '''
    选择第二个alpha值，
    '''
    maxK = -1
    maxDeltaE = 0
    Ej = 0
    oS.eCache[i] = [1, Ei]
    validEcacheList = np.nonzero(oS.eCache[:,0].A)[0]
    # 选择具有最大步长行的j
    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

def updateEk(oS, k):
    '''
    计算误差并存储，在对alpha值进行优化之后会用到
    '''
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1, Ek]

def innerL(i, oS):
    '''
    完整Platt SMO 算法中的优化例程
    有任意一对alpha发生改变 返回1
    '''
    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): 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
        return 1
    else: return 0

def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):    
    '''
    完整版Platte SMO 算法的外循环代码
    '''
    oS = optStruct(np.mat(dataMatIn), np.mat(classLabels).transpose(), C, toler) # 实例化 
    iter = 0
    entireSet = True
    alphaPairsChanged = 0
    # 当迭代超过maxIter或者整个数据集都未对任意alpha进行修改时，退出循环
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:   
            # 遍历所有值
            for i in range(oS.m): 
                alphaPairsChanged += innerL(i,oS) #选择第二个alpha
                print ("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
            iter += 1
        else:
            nonBoundIs = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
            for i in nonBoundIs: # 遍历所有非边界alpha值
                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

# 调用

import random
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
data = []
for i in range(50):
    x1 = random.randint(20,50)
    x2 = random.uniform(2,3)*x1
    if x1 > x2 :
        z = 1
    else:
        z = -1
    data.append([x1,x2,z])
for i in range(50):
    x1 = random.randint(20,50)
    x2 = random.uniform(0.5,0.8)*x1
    if x1 > x2 :
        z = 1
    else:
        z = -1
    data.append([x1,x2,z])
   
dataArray = data[:,0:2]
labelArray = data[:,2]

b,alphas = smoP(dataArray,labelArray,0.6,0.001,40)

运行结果

fullSet, iter: 0 i:0, pairs changed 1
fullSet, iter: 0 i:1, pairs changed 2
fullSet, iter: 0 i:2, pairs changed 2
fullSet, iter: 0 i:3, pairs changed 2
fullSet, iter: 0 i:4, pairs changed 3
fullSet, iter: 0 i:5, pairs changed 3
fullSet, iter: 0 i:6, pairs changed 3
L==H
fullSet, iter: 0 i:7, pairs changed 3
fullSet, iter: 0 i:8, pairs changed 3
fullSet, iter: 0 i:9, pairs changed 3
L==H
fullSet, iter: 0 i:10, pairs changed 3
fullSet, iter: 0 i:11, pairs changed 3
fullSet, iter: 0 i:12, pairs changed 3
……
fullSet, iter: 2 i:95, pairs changed 0
fullSet, iter: 2 i:96, pairs changed 0
fullSet, iter: 2 i:97, pairs changed 0
j not moving enough
fullSet, iter: 2 i:98, pairs changed 0
j not moving enough
fullSet, iter: 2 i:99, pairs changed 0
iteration number: 3

查看效果，绘制出支持向量

plt.scatter(data[:50,0],data[:50,1],c = 'pink') 
plt.scatter(data[50:,0],data[50:,1],c = 'green') 

x = []
for i in range(100):
    if alphas[i] > 0.0:
        print(dataArray[i],labelArray[i])
        x.append(dataArray[i])
x = np.array(x)
# 画散点图
plt.scatter(x[:,0], x[:,1],c ='black' , s=20, alpha=0.5)


plt.plot()
plt.show()

结果

[25.         67.37676254] -1.0
[29.         58.55461893] -1.0
[21.        45.3533745] -1.0
[47.         35.94429139] 1.0
[29.         20.26931489] 1.0

在这里插入图片描述为了最后预测，计算w值

def calcWs(alphas, dataArr, classLabels):
    """
    分类器，计算W
    alphas：数据向量
    dataArr：数据集
    classLabels：类别标签

    输出w：w*x+b中的w
    """ 
    X = np.mat(dataArr)
    labelMat = np.mat(classLabels).transpose()  # 转置
    m,n = shape(X)
    w = zeros((n,1))
    #遍历alpha，更新w
    for i in range(m):
        w += multiply(alphas[i]*labelMat[i],X[i,:].T)
    return w

利用核函数来实现更加复杂的情况

核函数的主要思想为将数据映射到更高维空间去，变线性不可分为线性可分
SVM优化中所有的运算都可以写成向量内积（点积）的形式。可以将这种向量内积替换成核函数，而不必再做其他处理。

径向基核函数——常用核函数

径向基函数采用向量作为自变量，能够基于向量距离运算出一个标量。
径向基高斯版本公式：
$e^{(\frac{- ||x-y||^{2}}{2\sigma ^{2}})}$
其中， $\sigma$ 是到达率，即函数值跌落到0的速度参数

def kernelTrans(X, A, kTup):
    '''
    核转换函数
    X：数据集矩阵
    A：数据集的一行
    kTup：核函数信息，第1个参数（字符串）代表不同类型的核函数，其他两个参数是可选的参数
    '''
    m,n = shape(X)
    K = mat(zeros((m,1))) 
    if kTup[0] == 'lin': K = X * A.T   # 线性核函数
    elif kTup[0] == 'rbf':  # 径向基函数rbf
        for j in range(m):
            deltaRow = X[j,:] - A
            K[j] = deltaRow*deltaRow.T
        K = exp(K/(-1*kTup[1]**2)) # 第2个参数是可选参数
    else: raise NameError('Houston We Have a Problem -- That Kernel is not recognized')  # 抛出异常用raise
    return K


class optStruct:
    '''
    创建类，用来保存重要值--仅作为数据结构使用
    较之前，引入了新变量kTup 包含核函数信息的元组 
    '''
    def __init__(self,dataMatIn, classLabels, C, toler, 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))) #first column is valid flag
        self.K = mat(zeros((self.m,self.m)))
        # 只需计算一次K
        for i in range(self.m):
            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)

            
    
def calcEk(oS, k):
    '''
    计算误差E值并返回
    '''   
    fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek
 
      
def selectJ(i, oS, Ei):         #this is the second choice -heurstic, and calcs Ej
    maxK = -1; maxDeltaE = 0; Ej = 0
    oS.eCache[i] = [1,Ei]  #set valid #choose the alpha that gives the maximum delta E
    validEcacheList = nonzero(oS.eCache[:,0].A)[0]
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:   #loop through valid Ecache values and find the one that maximizes delta E
            if k == i: continue #don't calc for i, waste of time
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k; maxDeltaE = deltaE; Ej = Ek
        return maxK, Ej
    else:   #in this case (first time around) we don't have any valid eCache values
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej

'''
计算误差值，并存入缓存中
'''
def updateEk(oS, k):#after any alpha has changed update the new value in the cache
    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) #this has been changed from selectJrand
        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] #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) #added this for the Ecache
        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])#update i by the same amount as j
        updateEk(oS, i) #added this for the Ecache                    #the update is in the oppostie direction
        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

'''
完整版Platt SMO的外循环代码，是内核版本
'''
def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):    #full Platt SMO
    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:   #go over all
            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:#go over non-bound (railed) alphas
            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 #toggle entire set loop
        elif (alphaPairsChanged == 0): entireSet = True  
        print("iteration number: %d" % iter)
    return oS.b,oS.alphas

'''
# 测试简化版SMO算法
dataArr, labelArr = loadDataSet('testSet.txt')
b, alphas = smoP(dataArr, labelArr, 0.6, 0.001, 40)
'''

'''
用计算出的alpha值对数据进行分类
@param alphas：计算出的alphas值向量，大多数为0，非零alpha即为支持向量
@param dataArr：
@param classLabels：
'''
def calcWs(alphas,dataArr,classLabels):
    X = mat(dataArr); labelMat = mat(classLabels).transpose()
    m,n = shape(X)
    w = zeros((n,1)) # w向量初始化为1
    for i in range(m):
        w += multiply(alphas[i]*labelMat[i],X[i,:].T)  # 计算w
    return w

'''
# 测试分类
dataArr, labelArr = loadDataSet('testSet.txt')
b, alphas = smoP(dataArr, labelArr, 0.6, 0.001, 40)
ws = calcWs(alphas,dataArr,labelArr)
print(ws)  # 打印w1、w2值
datMat = mat(dataArr)
label = datMat[0] * ws + b  # 对第1个数据点进行分类
print(label) # -0.925... -> -1
label = datMat[2] * ws + b  # 对第3个数据点进行分类
print(label) # 2.304... -> 1
'''

'''
在测试中使用核函数rbf
@param k1：用户定义的可选参数，即径向基函数中的σ，σ太小，支持向量的点越多，错误率会上升
'''
def testRbf(k1=1.3):
    dataArr,labelArr = loadDataSet('testSetRBF.txt') # 训练集
    b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) # C=200 important
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    svInd=nonzero(alphas.A>0)[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,:],('rbf', k1))
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1   # sign函数计算各元素的符号值，1（+），0，-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))
        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))    

# 测试径向基函数作用于训练集和测试集上的分类错误率
testRbf();  # there are 30 Support Vectors; the training error rate is: 0.130000; the test error rate is: 0.150000

weixin_44275584

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
《机器学习实战》第六章——支持向量机

支持向量机就是希望找到到支持向量距离最小的那一个超平面。优点：泛化错误率低，计算开销不大，容易得到结果缺点：对参数调解和核函数的选择敏感，原始分类器不加修改仅适用于处理二分类问题适用于：线性可分：可以使用一条线（超平面）完全把‘0’，‘1’两种类型的数据分割开——该线称为分割超平面，也就是分类的决策边界支持向量机就是希望能采取数据点离决策边界最远的方法来构建分类器。分类器工作原理：输入数据给分类器，会输出一个分类标签，这相当于类似sigmoid函数的作用。超平面确定后，各个点到超平面
复制链接

扫一扫