支持向量机书写体识注释，含文件

基于支持向量机的手写体识别

代码全是《机器学习》上的，只是将其整合到了一起，能够运行手写体识别。
内容大部分进行了注释，可能有些注释不够精准或者不容理解，见谅！

代码

from numpy import * #导入数据处理库
from time import sleep #导入时间间隙函数

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

def selectJrand(i,m): # 从0~m个alpha中选择除i外的另一个alpha
    j=i 
    while (j==i): #选择一个不等于i的j
        j = int(random.uniform(0,m))
    return j

def clipAlpha(aj,H,L): # 用来调整，使得它在H—L之间
    if aj > H: 
        aj = H
    if L > aj:
        aj = L
    return aj

def kernelTrans(X, A, kTup): # kTup元组中的第一个参数是描述所用核函数类型的一个字符串，其他2个参数是核函数可能需要的可选参数，数据矩阵、单个数据向量、包含和函数信息的元组
m,n = shape(X)
    K = mat(zeros((m,1)))
    if kTup[0]=='lin': K = X * A.T   #linear kernel
    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)) #在numpy中除法是按元素而不是像matlab那样的矩阵    
else: raise NameError('Houston We Have a Problem -- \
    That Kernel is not recognized')
    return K

class optStruct:
    def __init__(self,dataMatIn, classLabels, C, toler, kTup):  # Initialize the structure with the parameters 
        self.X = dataMatIn  # X:训练的数据集
        self.labelMat = classLabels  # labelMat:X对应的类别标签
        self.C = C  # C:松弛变量系数
        self.tol = toler  # tol:容错率
        self.m = shape(dataMatIn)[0]  # m:样本的个数
        self.alphas = mat(zeros((self.m,1)))  # alphas:拉格朗日系数，需要优化项
        self.b = 0  # b:阈值
        self.eCache = mat(zeros((self.m,2))) #first column is valid flag
        self.K = mat(zeros((self.m,self.m)))
        for i in range(self.m):  # 第一列 标志位,标志Ek是否有效，1为有效，0为无效 第二列 错误率Ek
            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)
            
def calcEk(oS, k): #计算预测值与真是值的误差
    fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b) #(alpha*label).t*k[:,k]+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]  #将Ei在缓存中设置成为已经计算好的
validEcacheList = nonzero(oS.eCache[:,0].A)[0] #构建一非零Ei值对应的alpha值 
    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
    else:   #否则，我们将没有有效的误差值
        j = selectJrand(i, oS.m)#随机选择一个j值
        Ej = calcEk(oS, j))#j值对应的误差值Ej
    return j, Ej

def updateEk(oS, k):#将计算出的误差存入缓存中 ，数据集、alpha
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1,Ek]
    
def innerL(i, oS): #完整的内循环优化进程， 输入选出的第一个alpha和相关的参数，输出是否进行了参数的优化
    Ei = calcEk(oS, i) #计算第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) #选择有最大的步长的alpha
        alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy();#保存更新前/后的aplpha值，使用深拷贝
        if (oS.labelMat[i] != oS.labelMat[j]): #如果yi和yj的标签不一样，alphas[j]new的取值范围
            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*Kij+Kii+Kjj，而这儿eta=2*Kij-Kii-Kjj,所以下面公式中用的是减号,也可以看做关于oS.alphas[j]的二阶导数
        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#目标函数是求最小值，这儿eta是负的二阶导数，因此二阶导数不能小于等于0
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta#更新oS.alphas[j]
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)#消减到规定范围内
        updateEk(oS, j) #added this for the Ecache#更新Ej至误差缓存
        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


def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):    #完整的线性SMO，数据矩阵、数据标签、松弛变量、容错率、最大迭代次数
    oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler, kTup)#初始化数据结构
    iter = 0 #初始化当前迭代次数
entireSet = True; alphaPairsChanged = 0
# 遍历整个数据集都alpha也没有更新或者超过最大迭代次数,则退出循环
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:  #遍历整个数据集#首先进行完整遍历，过程和简化版的SMO一样
            for i in range(oS.m):        
                alphaPairsChanged += innerL(i,oS)#使用优化的SMO算法，如果oS.alphas[j]和oS.alphas[i]更新，则返回1
                print ("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
            iter += 1#循环完一次，迭代次数加1
        else:#非边界遍历，挑选其中alpha值在0和C之间非边界alpha进行优化
            nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]#numpy.nonzeros返回非零元素的下标位置
            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 #如果alpha没有更新,计算全样本遍历 
        print ("iteration number: %d" % iter)
    return oS.b,oS.alphas

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 loadImages(dirName): #将文件夹下面的数字文件 转换成数据特征（每个数字1/0就是一个特征，每个文件有1024个特征）和类别标签
    from os import listdir  #方法用于返回指定的文件夹包含的文件或文件夹的名字的列表
    hwLabels = []
    trainingFileList = listdir(dirName)           #加载训练集
    m = len(trainingFileList)
    trainingMat = zeros((m,1024))
    for i in range(m):
        fileNameStr = trainingFileList[i]
        fileStr = fileNameStr.split('.')[0]     #对文件名进行分割。例如2_45.txt，从.那个地方开始分割文件名，就得到2_45和txt两部分
        classNumStr = int(fileStr.split('_')[0])#对文件名进行分割。例如2_45，从_那个地方开始分割文件名
        if classNumStr == 9: hwLabels.append(-1) #二分类，数字9标签为-1 其他的标签为+1
        else: hwLabels.append(1)
        trainingMat[i,:] = img2vector('%s/%s' % (dirName, fileNameStr))#读取数据到矩阵
    return trainingMat, hwLabels    

def testDigits(kTup=('rbf', 10)):
    dataArr,labelArr = loadImages('trainingDigits')  #加载训练数据
    b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, kTup)
    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,:],kTup)  #计算各个点的核
        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 = loadImages('testDigits')  #加载测试数据
    errorCount = 0
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    m,n = shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],kTup)
        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)) 

def img2vector(filename):  # 读取数据到矩阵
    returnVect = zeros((1,1024))  # 创建向量
    fr = open(filename)
    for i in range(32):
        lineStr = fr.readline()
        for j in range(32):  # 将每行前32字符转成int存入向量
            returnVect[0,32*i+j] = int(lineStr[j])
    return returnVect

运行结果

如果labelMat[i]>0，其属于1类，若<0，属于-1类，验证结果：

手写数字的识别：
首先是alpha=20

接着alpha=10：

最后alpha=5：

最后附上testDigits测试文件

地址：百度云盘
链接：https://pan.baidu.com/s/1qeOBsij1DxssLaU14tMh_g
提取码：pvpg

注意：由于文件太多，运行时间会比较长，请耐心等待。