支持向量机（四）：支持向量机的Python语言实现

最新推荐文章于 2023-07-08 19:44:24 发布

0x12A2A7F

最新推荐文章于 2023-07-08 19:44:24 发布

阅读量595

点赞数

分类专栏： Machine Learning Charms of Mathematics

Machine Learning 同时被 2 个专栏收录

53 篇文章 0 订阅

订阅专栏

Charms of Mathematics

13 篇文章 1 订阅

订阅专栏

转：http://www.cnblogs.com/pursued-deer/p/7892342.html

1 数据样本集的介绍

这篇文章是根据《机器学习实战》一书的实例进行代码的详细解读，我在查找这方面的资料没有人对支持向量机算法 python 实现的详细说明，我就把我在看代码时的思路和代码详细注解。如果存在不足，欢迎给我留言相互探讨。好了，废话不多说，正文开始。。。

首先我们使用的数据是二维的坐标点，还有对应的类标号（1 或 -1）。数据集以 “testSet.txt” 命名，如下代码段中：

数据集

然后，数据集中对象在下图显示，我们的工作就是找到最佳的超平面（这里是直线）。

2 Python 代码的实现

2.1 准备数据和分析数据

既然我们有了数据，那么就要把数据转换成可以被程序调用和处理的形式。写一个加载数据的函数，把文本中的数据转化成列表类型。代码如下：

# 将文本中的样本数据添加到列表中
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

由于在实现支持向量机的过程中要多次调用样本数据、类标号、对偶因子、对象经过核函数映射之后的值和其他数据，因此利用创造一个类，初始化这些数据利于保存和调用。代码如下：

class optStruct:
    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))) # 差值矩阵 m * 2，第一列是对象的标志位 1 表示存在不为零的差值 0 表示差值为零，第二列是实际的差值 E
        self.K = mat(zeros((self.m, self.m))) # 对象经过核函数映射之后的值
        for i in range(self.m): # 遍历全部样本集
            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup ) # 调用径向基核函数

2.2 对偶因子的求解函数

我们通过支持向量机（二）：SMO算法，可知对偶因子是我们求得超平面的关键，只有解得对偶因子才能得到超平面的法向量和截距。代码如下：

# 遍历所有能优化的 alpha
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup=('lin', 0)):
    oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler, kTup) # 创建一个类对象 oS ,类对象 oS 存放所有数据
    iter = 0 # 迭代次数的初始化
    entireSet = True # 违反 KKT 条件的标志符
    alphaPairsChanged = 0 # 迭代中优化的次数

    # 从选择第一个 alpha 开始，优化所有alpha
    while(iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet )): # 优化的终止条件：在规定迭代次数下，是否遍历了整个样本或 alpha 是否优化
        alphaPairsChanged  = 0
        if entireSet: #
            for i in range(oS.m): # 遍历所有对象
                alphaPairsChanged += innerL(i ,oS) # 调用优化函数（不一定优化）
            print "fullSet , iter: %d i %d, pairs changed %d" % (iter, i , alphaPairsChanged )
            iter += 1 # 迭代次数加 1
        else:
            nonBoundIs = nonzero((oS.alphas .A > 0) * (oS.alphas.A < C))[0]
            for i in nonBoundIs : # 遍历所有非边界样本集
                alphaPairsChanged += innerL(i, oS) # 调用优化函数（不一定优化）
                print "non-bound, iter: %d i :%d, pairs changed %d" % (iter, i, alphaPairsChanged )
            iter += 1 # 迭代次数加 1
        if entireSet : # 没有违反KKT 条件的alpha ，终止迭代
            entireSet = False
        elif (alphaPairsChanged == 0): # 存在违反 KKT 的alpha
            entireSet = True
        print "iteration number: %d" % iter
    return oS.b, oS.alphas # 返回截距值和 alphas

上述代码中，第 13 行调用了 innerL() 函数，它是对对偶因子向量的单个元素求解。代码如下：

# 优化选取两个 alpha ，并计算截距 b
def innerL(i, oS):
    Ei = calcEk(oS, i) # 计算 对象 i 的差值
    # 第一个 alpha 符合选择条件进入优化
    if ((oS.labelMat [i]*Ei <- oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat [i]*Ei > oS.tol) and (oS.alphas[i] > 0)):
        j,Ej =selectJ(i, oS, Ei) # 选择第二个 alpha
        alphaIold = oS.alphas[i].copy() # 浅拷贝
        alphaJold = oS.alphas[j].copy() # 浅拷贝

        # 根据对象 i 、j 的类标号（相等或不等）确定KKT条件的上界和下界
        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 # 不符合优化条件（第二个 alpha）
        eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j]  # 计算公式的eta ,是公式的相反数
        if eta >= 0:print "eta>=0";return 0 # 不考虑eta 大于等于 0 的情况（这种情况对 alpha 的解是另外一种方式，即临界情况的求解）
        # 优化之后的第二个 alpha 值
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
        updateEk(oS, j) # 更新差值矩阵
        if (abs(oS.alphas[j] - alphaJold) < 0.00001): # 优化之后的 alpha 值与之前的值改变量太小，步长不足
            print "j not moving enough"
            return 0
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j]) # 优化第二个 alpha
        updateEk(oS, i) # 更新差值矩阵
        # 计算截距 b
        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

这里第 3 行代码调用了calcEk() 函数，它是计算预测值与真实值的差值。代码如下：

 # 预测的类标号值与真实值的差值，参数 oS 是类对象，k 是样本的对象的标号
 def calcEk(oS, k):
     fXk = float(multiply(oS.alphas, oS.labelMat).T * oS.K[:, k] + oS.b)  # 公式（1）
     Ek = fXk - float(oS.labelMat[k]) # 差值
     return Ek

其中，第 3 行中指到的公式（1）是超平面方程。

在 innerL（）函数中的第 6 行，调用了 selectJ （）函数，此函数是用来选择优化的第二个对偶因子元素。代码如下：

# 由启发式选取第二个 alpha，以最大步长为标准
def selectJ(i, oS, Ei): # 函数的参数是选取的第一个 alpha 的对象号、类对象和对象差值
    maxK = -1; maxDeltaE = 0; Ej = 0 # 第二个 alpha 的初始化
    oS.eCache[i] = [1,Ei] # 更新差值矩阵的第 i 行
    validEcacheList = nonzero(oS.eCache[:,0].A)[0] # 取差值矩阵中第一列不为 0 的所有行数（标志位为 1 ），以元组类型返回
    if (len(validEcacheList)) > 1 : #
        for k in validEcacheList : # 遍历所有标志位为 1 的对象的差值
            if k == i: continue
            Ek = calcEk(oS, k) # 计算对象 k 的差值
            deltaE = abs(Ei - Ek) # 取两个差值之差的绝对值
            if (deltaE > maxDeltaE): # 选取最大的绝对值deltaE
                maxK = k; maxDeltaE = deltaE; Ej = Ek
        return maxK, Ej # 返回选取的第二个 alpha
    else: # 随机选取第二个 alpha
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS,j)
    return j, Ej # 返回选取的第二个 alpha

这个函数的第 15 行调用了 selectJrand() 函数，是为了没有合适的第二个元素，就用随机的方式选择。代码如下：

 # 随机选取对偶因子alpha ,参数i 是alpha 的下标，m 是alpha 的总数
 def selectJrand(i,m):
     j = i
     while (j==i):
         j = int(random.uniform(0,m))
     return j

在innerL（）函数中的第 23 行调用 clipAlpha() 函数，是根据KKT 条件对对偶因子的修剪。代码如下：

# 对所求的对偶因子按约束条件的修剪
def clipAlpha(aj, H, L): # H 为上界，L为下界
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

innerL（）函数的第 29 行调用了 updateEk() 函数，更新差值矩阵的数据。代码如下：

 # 更新差值矩阵的数据
 def updateEk(oS, k):
     Ek = calcEk(oS, k) # 调用计算差值的函数
     oS.eCache [k] = [1,Ek]

2.3 核函数

这里我们使用的径向基核函数，把数据从低维映射到高维，把不可分数据变成可分。当然，我们给出的数据集是线性可分，现在给出另外数据集，分为训练集“testSetRBF.txt”和测试集"testSetRBF2.txt"。

径向基核函数公式如下：

数据集

训练集
-0.214824    0.662756    -1.000000
-0.061569    -0.091875    1.000000
0.406933    0.648055    -1.000000
0.223650    0.130142    1.000000
0.231317    0.766906    -1.000000
-0.748800    -0.531637    -1.000000
-0.557789    0.375797    -1.000000
0.207123    -0.019463    1.000000
0.286462    0.719470    -1.000000
0.195300    -0.179039    1.000000
-0.152696    -0.153030    1.000000
0.384471    0.653336    -1.000000
-0.117280    -0.153217    1.000000
-0.238076    0.000583    1.000000
-0.413576    0.145681    1.000000
0.490767    -0.680029    -1.000000
0.199894    -0.199381    1.000000
-0.356048    0.537960    -1.000000
-0.392868    -0.125261    1.000000
0.353588    -0.070617    1.000000
0.020984    0.925720    -1.000000
-0.475167    -0.346247    -1.000000
0.074952    0.042783    1.000000
0.394164    -0.058217    1.000000
0.663418    0.436525    -1.000000
0.402158    0.577744    -1.000000
-0.449349    -0.038074    1.000000
0.619080    -0.088188    -1.000000
0.268066    -0.071621    1.000000
-0.015165    0.359326    1.000000
0.539368    -0.374972    -1.000000
-0.319153    0.629673    -1.000000
0.694424    0.641180    -1.000000
0.079522    0.193198    1.000000
0.253289    -0.285861    1.000000
-0.035558    -0.010086    1.000000
-0.403483    0.474466    -1.000000
-0.034312    0.995685    -1.000000
-0.590657    0.438051    -1.000000
-0.098871    -0.023953    1.000000
-0.250001    0.141621    1.000000
-0.012998    0.525985    -1.000000
0.153738    0.491531    -1.000000
0.388215    -0.656567    -1.000000
0.049008    0.013499    1.000000
0.068286    0.392741    1.000000
0.747800    -0.066630    -1.000000
0.004621    -0.042932    1.000000
-0.701600    0.190983    -1.000000
0.055413    -0.024380    1.000000
0.035398    -0.333682    1.000000
0.211795    0.024689    1.000000
-0.045677    0.172907    1.000000
0.595222    0.209570    -1.000000
0.229465    0.250409    1.000000
-0.089293    0.068198    1.000000
0.384300    -0.176570    1.000000
0.834912    -0.110321    -1.000000
-0.307768    0.503038    -1.000000
-0.777063    -0.348066    -1.000000
0.017390    0.152441    1.000000
-0.293382    -0.139778    1.000000
-0.203272    0.286855    1.000000
0.957812    -0.152444    -1.000000
0.004609    -0.070617    1.000000
-0.755431    0.096711    -1.000000
-0.526487    0.547282    -1.000000
-0.246873    0.833713    -1.000000
0.185639    -0.066162    1.000000
0.851934    0.456603    -1.000000
-0.827912    0.117122    -1.000000
0.233512    -0.106274    1.000000
0.583671    -0.709033    -1.000000
-0.487023    0.625140    -1.000000
-0.448939    0.176725    1.000000
0.155907    -0.166371    1.000000
0.334204    0.381237    -1.000000
0.081536    -0.106212    1.000000
0.227222    0.527437    -1.000000
0.759290    0.330720    -1.000000
0.204177    -0.023516    1.000000
0.577939    0.403784    -1.000000
-0.568534    0.442948    -1.000000
-0.011520    0.021165    1.000000
0.875720    0.422476    -1.000000
0.297885    -0.632874    -1.000000
-0.015821    0.031226    1.000000
0.541359    -0.205969    -1.000000
-0.689946    -0.508674    -1.000000
-0.343049    0.841653    -1.000000
0.523902    -0.436156    -1.000000
0.249281    -0.711840    -1.000000
0.193449    0.574598    -1.000000
-0.257542    -0.753885    -1.000000
-0.021605    0.158080    1.000000
0.601559    -0.727041    -1.000000
-0.791603    0.095651    -1.000000
-0.908298    -0.053376    -1.000000
0.122020    0.850966    -1.000000
-0.725568    -0.292022    -1.000000

测试集
0.676771    -0.486687    -1.000000
0.008473    0.186070    1.000000
-0.727789    0.594062    -1.000000
0.112367    0.287852    1.000000
0.383633    -0.038068    1.000000
-0.927138    -0.032633    -1.000000
-0.842803    -0.423115    -1.000000
-0.003677    -0.367338    1.000000
0.443211    -0.698469    -1.000000
-0.473835    0.005233    1.000000
0.616741    0.590841    -1.000000
0.557463    -0.373461    -1.000000
-0.498535    -0.223231    -1.000000
-0.246744    0.276413    1.000000
-0.761980    -0.244188    -1.000000
0.641594    -0.479861    -1.000000
-0.659140    0.529830    -1.000000
-0.054873    -0.238900    1.000000
-0.089644    -0.244683    1.000000
-0.431576    -0.481538    -1.000000
-0.099535    0.728679    -1.000000
-0.188428    0.156443    1.000000
0.267051    0.318101    1.000000
0.222114    -0.528887    -1.000000
0.030369    0.113317    1.000000
0.392321    0.026089    1.000000
0.298871    -0.915427    -1.000000
-0.034581    -0.133887    1.000000
0.405956    0.206980    1.000000
0.144902    -0.605762    -1.000000
0.274362    -0.401338    1.000000
0.397998    -0.780144    -1.000000
0.037863    0.155137    1.000000
-0.010363    -0.004170    1.000000
0.506519    0.486619    -1.000000
0.000082    -0.020625    1.000000
0.057761    -0.155140    1.000000
0.027748    -0.553763    -1.000000
-0.413363    -0.746830    -1.000000
0.081500    -0.014264    1.000000
0.047137    -0.491271    1.000000
-0.267459    0.024770    1.000000
-0.148288    -0.532471    -1.000000
-0.225559    -0.201622    1.000000
0.772360    -0.518986    -1.000000
-0.440670    0.688739    -1.000000
0.329064    -0.095349    1.000000
0.970170    -0.010671    -1.000000
-0.689447    -0.318722    -1.000000
-0.465493    -0.227468    -1.000000
-0.049370    0.405711    1.000000
-0.166117    0.274807    1.000000
0.054483    0.012643    1.000000
0.021389    0.076125    1.000000
-0.104404    -0.914042    -1.000000
0.294487    0.440886    -1.000000
0.107915    -0.493703    -1.000000
0.076311    0.438860    1.000000
0.370593    -0.728737    -1.000000
0.409890    0.306851    -1.000000
0.285445    0.474399    -1.000000
-0.870134    -0.161685    -1.000000
-0.654144    -0.675129    -1.000000
0.285278    -0.767310    -1.000000
0.049548    -0.000907    1.000000
0.030014    -0.093265    1.000000
-0.128859    0.278865    1.000000
0.307463    0.085667    1.000000
0.023440    0.298638    1.000000
0.053920    0.235344    1.000000
0.059675    0.533339    -1.000000
0.817125    0.016536    -1.000000
-0.108771    0.477254    1.000000
-0.118106    0.017284    1.000000
0.288339    0.195457    1.000000
0.567309    -0.200203    -1.000000
-0.202446    0.409387    1.000000
-0.330769    -0.240797    1.000000
-0.422377    0.480683    -1.000000
-0.295269    0.326017    1.000000
0.261132    0.046478    1.000000
-0.492244    -0.319998    -1.000000
-0.384419    0.099170    1.000000
0.101882    -0.781145    -1.000000
0.234592    -0.383446    1.000000
-0.020478    -0.901833    -1.000000
0.328449    0.186633    1.000000
-0.150059    -0.409158    1.000000
-0.155876    -0.843413    -1.000000
-0.098134    -0.136786    1.000000
0.110575    -0.197205    1.000000
0.219021    0.054347    1.000000
0.030152    0.251682    1.000000
0.033447    -0.122824    1.000000
-0.686225    -0.020779    -1.000000
-0.911211    -0.262011    -1.000000
0.572557    0.377526    -1.000000
-0.073647    -0.519163    -1.000000
-0.281830    -0.797236    -1.000000
-0.555263    0.126232    -1.000000

核函数代码如下：

# 径向基核函数（高斯函数）
def kernelTrans(X, A, kTup): # X 是样本集矩阵，A 是样本对象（矩阵的行向量） ， kTup 元组
    m,n = shape(X)
    K = mat(zeros((m,1)))
    # 数据不用核函数计算
    if kTup [0] == 'lin': K = X * A.T

    # 用径向基核函数计算
    elif kTup[0] == 'rbf':
        for j in range(m):
            deltaRow = X[j,:] - A
            K[j] = deltaRow * deltaRow.T
        K = exp(K/(-1*kTup[1]**2))
    # kTup 元组值异常，抛出异常信息
    else:raise NameError('Houston We Have a Problem --That Kernel is not recognized')
    return K

3 支持向量机预测分类测试结果

测试函数代码如下：

# 训练样本集的错误率和测试样本集的错误率
def testRbf(k1=1.3):
    dataArr,labelArr = loadDataSet('testSetRBF.txt') # 训练样本的提取
    b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) # 计算得到截距和对偶因子
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    svInd=nonzero(alphas.A>0)[0] # 对偶因子大于零的值，支持向量的点对应对偶因子
    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)) # 对象 i 的映射值
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b # 预测值
        if sign(predict)!=sign(labelArr[i]): errorCount += 1
    print "the training error rate is: %f" % (float(errorCount)/m)
    dataArr,labelArr = loadDataSet('testSetRBF2.txt') # 测试样本集的提取
    errorCount = 0
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    m,n = shape(datMat)
    # 对测试样本集的测试
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1)) # 测试样本对象 i 的映射值
        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)

运行结果：

从结果可知，训练集的错误率在3%，测试集的错误率在4%，说明支持向量机的分类效果很好。当然在每次运行的结果与图中的结果不一定相同，因为代码段中存在随机选择对偶因子的情况，这是无法避免的，但结果的错误率是基本相同的。

到这里支持向量机算法是讲完了，期间花费了我大量时间学习支持向量机的原理和代码的实现。虽然这部分内容存在一定的难度，但认真地去看每步的由来还是能够理解的。像我这么不怎么聪明的人都能对它说出 1,2,3,4 来，相信你也可以。我总是这么认为，你看懂支持向量机是一回事，能够说出向量机是一回事，把向量机讲出来又是另外一回事。我现在对支持向量机还有很多半只不懂的地方，希望以后对它有了新的理解再来修改支持向量机系列的文章。

0x12A2A7F

关注

0
点赞
踩
3

收藏

觉得还不错? 一键收藏
0
评论
支持向量机（四）：支持向量机的Python语言实现

转：http://www.cnblogs.com/pursued-deer/p/7892342.html1 数据样本集的介绍这篇文章是根据《机器学习实战》一书的实例进行代码的详细解读，我在查找这方面的资料没有人对支持向量机算法 python 实现的详细说明，我就把我在看代码时的思路和代码详细注解。如果存在不足，欢迎给我留言相互探讨。好了，废话不多说，正文开始。。。首先我们使用的数据是二...
复制链接

扫一扫