机器学习理论学习：基于python的SVM实践

通过前面两篇文章的学习，这里我们简单的用python代码实现SVM,通过前面两篇文章，我们知道在求解SVM模型之前我们需要知道alpha（拉格朗日乘子），然后计算w和b，这样才能得到模型，其相应的公式如下：

计算alpha：

计算w:

计算b：

（1）定义存储所有的数据的数据结构

 

'''
功能：计算核函数
输入：X-整个数据集
    A-第i行数据
    kTup-核信息元组
输出：k-核矩阵
'''        
def kernelTrans(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))
    else:
        raise NameError('Houston we Have a problem--that kernel is not recognized')
    return k;

class optStructK:
    def __init__(self,dataMatIn,classLabels,C,toler,kTup):
        self.X=dataMatIn#样本
        self.labelMat=classLabels#标签
        self.C=C#alpha的上边界
        self.tol=toler#容错率
        self.m=shape(dataMatIn)[0]#样本的数量
        self.alphas=mat(zeros((self.m,1)))#alphas-拉格朗日乘子
        self.b=0#模型的b值
        self.eCache=mat(zeros((self.m,2)))#误差缓存
        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）SMO优化

 

 

'''
功能：计算E值
输入：oS-输入数据
    k-第k个样本
输出：Ek-预测值与标签的差值
'''
def calcEk(oS,k):
    fXk=float(multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T))+oS.b
    Ek=fXk-float(oS.labelMat[k])
    return Ek

'''
功能：用于调整alphas的大小
'''
def clipAlpha(aj,H,L):
    if aj>H:
        aj=H
    if aj<L:
        aj=L
    return aj

'''
功能：更新Ek
'''
def updateEk(oS,k):
    Ek=calcEk(oS,k)
    oS.eCache[k]=[1,Ek]

'''
功能：内循环中启发式方法选择第二个alphas的下标
输入：i-第一个alphas的下标
    oS-输入数据
    Ei-第i个样本的误差
输出：maxK-具有最大变化的第二个alphas的下标
    Ej-第二个alphas的预测值与标签的差值
'''
def selectJ(i,oS,Ei):
    #初始化
    maxK=-1;maxDeltaE=0;Ej=0
    oS.eCache[i]=[1,Ei]
    
    #将变化量转化为非零表
    validEcacheList=nonzero(oS.eCache[:,0].A)[0]
    
    #如果存在第二个alphas
    if(len(validEcacheList))>1:
        #遍历所有alphas使得alphas变化最大的下标
        for k in validEcacheList:
            if k==i:continue
            Ek=calcEk(oS,k)
            deltaE=abs(Ei-Ek)
            #选择具有最大步长的j
            if(deltaE>maxDeltaE):
                maxK=k;maxDeltaE=deltaE;Ej=Ek
        return maxK,Ej
    #如果不存在alphas，则采用随机选取法
    else:
        j=selectJrand(i,oS.m)
        Ej=calcEk(oS,j)
    return j,Ej

    
'''
功能：随机选择第二个alpha j的下标
输入：i-第一个alphas的下标
    m-样本数量
输出：j-第二个alphas的下标
'''
def selectJrand(i,m):
    j=i
    while(j==i):
        j=int(random.uniform(0,m))
    return j

'''
功能：计算E值
输入：oS-输入数据
    k-第k个样本
输出：Ek-预测值与标签的差值
'''
def calcEkK(oS,k):
    fXk=float(multiply(oS.alphas,oS.labelMat).T*(oS.K[:,k])+oS.b)
    Ek=fXk-float(oS.labelMat[k])
    return Ek

'''
功能：计算两个变量的二次规划
输入：i-第一个变量alphas
    oS-输入数据
'''
def innerLK(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)#使用启发式方法选取第二个alphas
        alphaIold=oS.alphas[i].copy();alphaJold=oS.alphas[j].copy()#初始化两个alphas
        
        #计算第二个alphas的剪切边界
        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
        
        #计算2K12-K11-K22
        eta=2.0*oS.K[i,j]-oS.K[i,i]-oS.K[j,j]
        if eta>=0:
            print 'eta>=0';return 0
        #计算第二个alphas，并对其进行约束
        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 enought';return 0
        #计算第一个alphas
        oS.alphas[i]+=oS.labelMat[j]*oS.labelMat[i]*(alphaJold-oS.alphas[j])
        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]
        #如果0<a1<C，则b=b1
        if(0<oS.alphas[i])and(oS.C>oS.alphas[i]):
            oS.b=b1
        #如果0<a2<C，则b=b2
        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

'''
功能：SMO算法
输入：dataMatIn-输入样本
    classLabels-输入标签
    C-alphas的上边界
    toler-误差
    maxIter-最大迭代次数
输出：b-b
    alphas-alphas
'''    
def smoPK(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)): 
    oS = optStructK(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 += innerLK(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 += innerLK(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

（3）使用SVM

 

 

'''
功能：计算w
'''              
def calcWs(alphas,dataArr,classLabels):
    X=mat(dataArr)
    labelMat=mat(classLabels).transpose()
    m,n=shape(X)
    w=zeros((n,1))
    for i in range(m):
        w+=multiply(alphas[i]*labelMat[i],X[i,:].T)
    return w   

def testRbf(k=1.3):
    dataArr,labelArr=loadDataSet('./testSetRBF.txt')
    
    b,alphas=smo.smoPK(dataArr,labelArr,200,0.001,10000,('rbf',k))
    dataMat=mat(dataArr);labelMat=mat(labelArr).transpose()
    svInd=nonzero(alphas.A>0)[0]
    sVs=dataMat[svInd]
    labelSV=labelMat[svInd]
    print 'there are %d Support Vectors'%shape(sVs)[0]
    m,n=shape(dataMat)
    errorCount=0
    for i in range(m):
        kernelEval = smo.kernelTrans(sVs,dataMat[i,:],('rbf', k))
        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)
    
    b,alphas=smo.smoPK(dataArr,labelArr,200,0.001,10000,('rbf',k))
    dataMat=mat(dataArr);labelMat=mat(labelArr).transpose()
    svInd=nonzero(alphas.A>0)[0]
    sVs=dataMat[svInd]
    labelSV=labelMat[svInd]
    print 'there are %d Support Vectors'%shape(sVs)[0]
    m,n=shape(dataMat)
    errorCount=0
    for i in range(m):
        kernelEval = smo.kernelTrans(sVs,dataMat[i,:],('rbf', k))
        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)
    
    dataArr1,labelArr1 = loadDataSet('testSetRBF2.txt')
    errorCount = 0
    datMat=mat(dataArr1); labelMat = mat(labelArr1).transpose()
    m,n = shape(datMat)
    for i in range(m):
        kernelEval = smo.kernelTrans(sVs,datMat[i,:],('rbf', k))
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr1[i]): errorCount += 1    
    print "the test error rate is: %f" % (float(errorCount)/m)    
    
    #显示数据集
    dataMat0 = array(dataArr)
    n = shape(dataMat0)[0] 
    xcord1 = []; ycord1 = []
    xcord2 = []; ycord2 = []
    for i in range(n):
        if int(labelArr[i])== 1:
            xcord1.append(dataMat0[i,0]); ycord1.append(dataMat0[i,1])
        else:
            xcord2.append(dataMat0[i,0]); ycord2.append(dataMat0[i,1])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
    ax.scatter(xcord2, ycord2, s=30, c='green')
    for i in range(n):
        if alphas[i]>0.0:
            circle = Circle((float(dataArr[i][0]), float(dataArr[i][1])), 0.1, facecolor='blue', \
                edgecolor=(0,0.8,0.8), linewidth=3, alpha=0.5)
            print dataArr[i],labelArr[i]
            ax.add_patch(circle)
    plt.show()