機械学習ノート——ベクトルマシンSMOアルゴリズムの完全版コード(核関数)をサポートします.

55432 ワード

機械学習ノート——ベクトルマシンSMOアルゴリズムの完全版コード(核関数)をサポートします.
\Qquadの前のコードに使用される「間隔」は、ベクトル内積(K i j=(x i,x j)K_{ij}=(xuui,xuj)Kij=(xi,x xj)です.
<Qquadがサンプルデータの線形不可分の場合、核関数によって変換され、低次元データをより上位の中にマッピングして分類する必要があるので、アルゴリズムのコードを修正します.
\qquadは「間隔」計算方法を修正しましたので、つまりコードにK i j K_が含まれています.{ij}Kijの計算は全部修正しなければなりません.η 、 b、u i\eta、b、u u u uiη、b、ui、k e r n e l T r a n s kersnelTrans核変換関数を追加しました.ここにガウス核のみを追加しました.他の核関数を使って自分で追加する必要があります.以下は修正後のコードです.
# -*- coding: utf-8 -*-
"""
Created on Thu Aug 22 15:07:35 2019

@author:wangtao_zuel

E-mail:[email protected]

"""

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

def svmSmo(dataPath,outerPath,C=0.6,toler=0.001,maxIters=40,kTup=('lin',0)):
    """
       、   
        kTup      ,            ,                
            'lin':  
            'rbf':   
    """
    #     ,     xlsx    ,              
    dataMat,labelMat = loadData(dataPath)
    #          ,             
    oS = parameter(dataMat,labelMat,C,toler,kTup)
    #        
    iters = 0
    #       、     
    entireSet = True
    #      、             ,               
    alphaPairsChanged =0
    #       :                 (   、             )
    while (iters < maxIters) and ((alphaPairsChanged > 0) or (entireSet)):
        #         
        alphaPairsChanged =0
        #       alpha    0,          
        if entireSet:
            for i in range(oS.m):
                alphaPairsChanged += innerL(oS,i)
            print("fullSet. iters:%d. alphaPairsChanged:%d"%(iters,alphaPairsChanged))
            iters += 1
        else:
            #  nonzero          : alpha!=0oralpha!=C   
            nonBond = np.nonzero((oS.alpha.A > 0)*(oS.alpha.A < oS.C))[0]
            for i in nonBond:
                alphaPairsChanged += innerL(oS,i)
            print("nonBond. iters:%d. alphaPairsChanged:%d"%(iters,alphaPairsChanged))
            iters += 1
        #     、     ,                
        #       ,         ;                 ,       ,     True       False
        if entireSet:
            entireSet = False
        #                        ,         ,    entireSet   True,          、        
        elif alphaPairsChanged == 0:
            entireSet = True
    print("      !")
#    print(oS.alpha)
    #     
    w = calcW(oS)
#    print(w)
#    print(oS.b)
    #    ,        
#    draw(oS,w)
    #      ,        
#    parameterAnalyze(oS,w)
    #         
    predict(w,oS.b,outerPath)
#    return w[0,0],w[1,0],oS.b[0,0]

def loadData(datapath):
    """
        
    """
    data = pd.read_excel(datapath)
    #               
    features = np.mat(data.iloc[:,:-1])
    labels = np.mat(data.iloc[:,-1]).T
    
    return features,labels

class parameter:
    """
        :toler           ;C     ;eCache    Ei,     j      
    """
    def __init__(self,dataMat,labelMat,C,toler,kTup):
        self.x = dataMat
        self.y = labelMat
        self.m = dataMat.shape[0]
        self.alpha = np.mat(np.zeros((self.m,1)))
        self.b = 0
        self.C = C
        self.toler = toler 
        self.eCache = np.mat(np.zeros((self.m,2)))
        self.K = np.mat(np.zeros((self.m,self.m)))
        for ii in range(self.m):
            self.K[:,ii] = kernelTrans(self.x,self.x[ii,:],kTup)

def kernelTrans(X,A,kTup):
    """
         :kTup       ('lin':  ,'rbf':   )
    """
    m,n = X.shape
    K = np.mat(np.zeros((m,1)))
    #   
    if kTup[0] == 'lin':
        K = X*A.T
    #    
    elif kTup[0] == 'rbf':
        for jj in range(m):
            deltaRow = X[jj,:] - A
            K[jj,0] = deltaRow*deltaRow.T
        K = np.exp(K/(-2*kTup[1]**2))
    else:
        raise NameError('       !')
    
    return K

def innerL(oS,i):
    """
          ,       1;  KKT  、      0
    """
    Ei = calcEi(oS,i)
    #       KKT  ,         
    if ((oS.alpha[i,0] < oS.C) and (oS.y[i,0]*Ei < -oS.toler)) or ((oS.alpha[i,0] > 0) and (oS.y[i,0]*Ei > oS.toler)):
        #        j
        j,Ej = selectJ(oS,i,Ei)
        #         alpha,  alpha       
        alphaIOld = oS.alpha[i,0].copy()
        alphaJOld = oS.alpha[j,0].copy()
        #   alpha   
        if oS.y[i,0] != oS.y[j,0]:
            L = max(0,oS.alpha[j,0]-oS.alpha[i,0])
            H = min(oS.C,oS.C+oS.alpha[j,0]-oS.alpha[i,0])
        else:
            L = max(0,oS.alpha[j,0]+oS.alpha[i,0]-oS.C)
            H = min(oS.C,oS.alpha[j,0]+oS.alpha[i,0])
        #  L=H, alpha      ,       ,     0 
        if L == H:
            return 0
        eta = oS.K[i,i] + oS.K[j,j] - 2*oS.K[i,j]
        #  eta 0,   0,       0,  eta      
        if eta == 0:
            return 0
        #      ,       ,         alpha  0 ,         
        oS.alpha[j,0] += oS.y[j,0]*(Ei-Ej)/eta
        #             
        oS.alpha[j,0] = clipAlpha(oS.alpha[j,0],H,L)
        #   eChache
        updateEi(oS,j)
        #          ,   0
        if abs(oS.alpha[j,0]-alphaJOld) < 0.00001:
#            print("      ,      !")
            return 0
        #   alpha_i,  alpha_i*y_i alpha_j*y_j               
        oS.alpha[i,0] += oS.y[i,0]*oS.y[j,0]*(alphaJOld-oS.alpha[j,0])
        updateEi(oS,i)
        #     b
        bi = oS.b - Ei - oS.y[i,0]*(oS.alpha[i,0]-alphaIOld)*oS.K[i,i] - oS.y[j,0]*(oS.alpha[j,0]-alphaJOld)*oS.K[i,j]
        bj = oS.b - Ej - oS.y[i,0]*(oS.alpha[i,0]-alphaIOld)*oS.K[i,j] - oS.y[j,0]*(oS.alpha[j,0]-alphaJOld)*oS.K[j,j]
        #   b,     b            
        if (oS.alpha[i,0] > 0) and (oS.alpha[i,0] < oS.C):
            oS.b = bi
        elif (oS.alpha[j,0] > 0) and (oS.alpha[j,0] < oS.C):
            oS.b = bj
        else:
            oS.b = (bi+bj)/2
        #            1
        return 1
    #    KKT  ,   0
    else:
        return 0

def updateEi(oS,i):
    """
      eChache
    """
    Ei = calcEi(oS,i)
    oS.eCache[i,:] = [1,Ei]
        
def clipAlpha(alpha,H,L):
    """
    alpha      
    """
    if alpha < L:
        alpha = L
    elif alpha > H:
        alpha = H
    
    return alpha        

def selectJ(oS,i,Ei):
    """
       i       j
    """
    maxJ = 0
    maxdeltaE = 0
    oS.eCache[i,:] = [1,Ei]
    validEcacheList = np.nonzero(oS.eCache[:,0].A)[0]
    #     j        j
    if len(validEcacheList) > 1:
        for j in validEcacheList:
            if j == i:
                continue
            Ej = calcEi(oS,j)
            deltaE = abs(Ei-Ej)
            if deltaE > maxdeltaE:
                maxJ = j
                maxdeltaE = deltaE
                best_Ej = Ej
        
        return maxJ,best_Ej
    #        j,         j
    else:
        j = randomJ(i,oS.m)
        Ej = calcEi(oS,j)
        
        return j,Ej

def randomJ(i,m):
    """
        j
    """
    j = i
    while j == i:
        j = np.random.randint(0,m+1)
    
    return j

def calcEi(oS,i):
    """
      Ei,            j
    """
    ui = float(np.multiply(oS.alpha,oS.y).T*oS.K[:,i]+oS.b)
    Ei = ui - float(oS.y[i,0])
    
    return Ei

def calcW(oS):   
    """
        w
    """
    w = oS.x.T*np.multiply(oS.alpha,oS.y)
    
    return w

def draw(oS,w):
    """
           :                
    """
    x1 = []
    y1 = []
    x2 = []
    y2 = []
    for i in range(oS.m):
        if oS.y[i,0] == -1:
            x1.append(oS.x[i,0])
            y1.append(oS.x[i,1])
        else:
            x2.append(oS.x[i,0])
            y2.append(oS.x[i,1])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(x1,y1,marker='*')
    ax.scatter(x2,y2)
    x = np.arange(3,22,0.5)
    y = -(w[0,0]*x+oS.b[0,0])/w[1,0]
    ax.plot(x,y)
    plt.show()
    print(w)
    print(oS.b)

def parameterAnalyze(k,b,c):
    """
        :           
    """
    fig = plt.figure()
    ax1 = fig.add_subplot(211)
    ax1.plot(c,k)
    ax2 = fig.add_subplot(212)
    ax2.plot(c,b)
    plt.tight_layout()
    plt.show()

def predict(w,b,outerPath):
    """
            
    """
    data = pd.read_excel(outerPath)
    #        
    dataMat = np.mat(data)
    result = dataMat*w + b
    #       0    1,  0    -1 
    result[result>0] = 1
    result[result<0] = -1
    #     
    data['classLabel'] = result
    data.to_excel(outerPath,index=False)
    print("    !")