参考《机器学习实战》
一、KNN
原理：
步骤：
数据处理归一化(x-min)/(max-min)计算已知类别数据集中的点与当前点之间的距离选取与当前点距离最小的k个点确定k个点所在类别出现的频率，返回频率最高的类别作为当前点的预测分类
def classify0(inX,dataSet,labels,k):dataSetSize=dataSet,shape[0]##tile函数用来复制一份inX数组的副本##dataSet是数据集，我们要求输入数据与数据集的欧氏距离diffMat=tile(inX,(dataSetSize,1))-dataSetsqDiffMat=diffMat**2##axis=1表示用列相加sqDistance=sqDiffMat.sum(axis=1)distance=sqDistance**0.5##对数据按从小到大的次序排序,选择距离最小的k个点sortedDistIndicies=distances.argsort()classCount={}for i in range(k):voteIlabel=labels[sortedDistIndicies[i]]##get的赋值语句，建立新的字典键对，统计每个label出现的频率##dict.get(key, default=None)##如果之前classCount[voteIlabel]有值就返回这个值##否则新建一个字典键对，默认返回0classCount[voteIlabel]=sorted(classCount.get(voteIlabel,0)+1##对ClassCount的内容进行排序，就可以找到距离最近的k个目标中##出现频率最高的标签啦sortedClassCount=sorted(classCount.iteritems(),key=operator.itemgetter(1),reverse=True)return sortedClassCount[0][0]
测试部分
已有数据的一部分数据作为训练样本，另一部分作为测试样本。
这里我们用错误率作为测试的依据。
对数据的预处理是很重要的
例子：手写数字识别系统
def handwritingClassTest():hwLabels=[]##listdir用于返回指定绝对地址中目录的条目##存储进trainingFileList这个列表中trainingFileList=listdir('trainingDigits')m=len(trainingFileList)##创建一个m行1024列的训练矩阵，每行存储一个数字的图像trainingMat=zeros({m,1024})for i in range(m):##区分当前行是哪个数字的图像fileNameStr=trainingFileList[i]##split函数以符号为基准切分字符串，返回切分结果列表fileStr=fileNameStr.split('.')[0]classNumStr=int(fileStr.split('_')[0])hwLabels.append(classNumStr)##img2vector:首先创建1×1024的NumPy数组，然后打开给定的文件，##循环读出文件的前32行，并将每行的头32个字符值存储在NumPy数组trainingMat[i,:]=img2vector('trainingDigits/%s'%fileNameStr')##下面用同样的方式处理测试集testFileList=listdir('testDigits')errorCount=0.0mTest=len(testFileList)for i in range(mTest)fileNameStr.testFileList[i]fileStr=fileNameStr.split('.')[0]classNumStr=int(fileStr.split('_')[0])vectorUnderTest=img2vector('testDigits/%s' % fileNameStr)##用上面的classify0去跑测试集和训练集classifierResult=classify0(vectorUnderTest,\trainingMat,hwLabels,3)##打印每一次测试集跑出的结果print" the classifier came back with: %d,the real answer is:%d"\ %(classifierResult,classNumStr)##统计计算错误率if(classifierResult != classNumStr):errorCount+=1.0print"\n the total number of errors is :%d" %errorCountprint"\n the total error rate is %f"% (errorCount/float(mTest))
二、决策树
步骤：
检测是否属于一个分类用信息论寻找划分数据集的最好特征划分数据集对每个划分子集继续1—3的操作
不同的划分数据集的算法是重点
划分数据集之前之后信息发生的变化称为信息增益。获得信息增益最高的特征就是最好的选择。
信息期望值计算公式：
−∑i=1np(xi)log2p(xi)-\sum_{i=1}^np(x_i)log_2p(x_i) −i=1∑np(xi)log2p(xi)
计算概率后计算信息熵的代码：
from math import logdef calcShannonEnt(dataSet)numEntries=len(dataSet)labelCounts={}##做一个键是类别，值是出现次数的字典for featVec in dataSet:currentLabel=featVec[-1]if currentLabel not in labelCounts.keys():labelCounts[currentLabel]=0labelCounts[currentLabel]+=1shannonEnt=0.0for key in labelCounts:prob=float(labelCounts[key])/numEntriesshannonEnt-=prob*log(prob,2)return shannonEnt
ID3
按照给定特征划分数据集
##list.append(object) 向列表中添加一个对象object##list.extend(sequence) 把一个序列seq的内容添加到列表中def splitDataSet(dataSet,axis,value):retDataSet=[]for featVec in dataSet:if featVec[axis]==value##把已经用来分类的这个属性去掉reducedFeatVec=featVec[:axis]reducedFeatVec.extend(featVec[axis+1:])##添加到分类结果列表retDataSet.append(reducedFeatVec)return retDataSet
选择最好的数据集划分方式
def chooseBeatFeatureToSplit(dataSet):numFeatures=len(dataSet[0])-1baseEntropy=calcShannonEnt(dataSet)bestInfoGain=0.0bestFeature=-1for i in range(numFeatures):##将列表中第i个特征值有可能存在的值放入新的列表featList=[example[i] for examples in dataSet]##featlist去重uniqueVals=set(featList)newEntropy=0.0##对每一个唯一属性值划分数据集for value in uniqueVals：##按照给定特征划分数据集subDataSet=splitDataSet(data,i,value)##计算数据集的新熵值prob=len(subDataSet)/float(len(dataSet))newEntropy+=prob*calcShannonEnt(SubDataSet)##计算这样划分得到的信息增益，返回最好特征划分的索引值infoGain=baseEntropy-newEntropyif(infoGain>bestInfoGain):bestInfoGain=infoGainbestFeature=-1return bestFeature
在一棵决策树上，我们需要
得到数据划分数据在划分的子集上再次划分数据如果程序遍历完所有划分数据集的属性（之后采用多数表决方法确定该叶子节点的分类）
或者每个分支下的每个实例都有相同的分类
则停止递归
def createTree(dataSet,labels):classList=[example[-1] for example in dataSet]##类别完全相同则停止划分if classList.count(classList[0])==len(classList)return classList[0]##遍历完所有特征，返回出现次数最多的类别if len(dataSet[0]==-1)return majorityCnt(classList)##选出最好的特征用于分割bestFeat=chooseBestFeatureToSplit(dataSet)bestFeatLabel=labels[bestFeat]myTree={bestFeatLabel:{}}##在集合里删除这个特征del(labels[bestFeat])##得到列表包含的所有属性值featValues=[example[bestFeat] for example in dataSet]uniqueVals=set(featValues)##对属性值中的每个属性创建子树for value in uniqueVals:subLabels=labels[:]myTree[bestFeatLabel][value]=createTree(splitDataSet\(dataSet,bestFeat,value),subLabels)return myTree
最终得到的结果：
一个子树：一个字典
一个叶子结点：一个类标签
绘图
获取叶节点的数目
def getNumLeafs(myTree):numLeafs=0firstStr=myTree.keys()[0]##找出树种第一个关键字secondDict=myTree[firstStr]for key in secondDict.keys()##type()判定测试节点的数据是否为字典##name属性表示类的名称if type(secondDict[key]).__name__=='dict':##递归调用getNumLeafsnumLeafs+=getNumLeafs(secondDict[key])elsenumLeafs+=1return numLeafs
获取树的层数
def getTreeDepth(myTree):maxDepth=0firstStr=myTree.keys()[0]secondDict=myTree[firstStr]for key in secondDict.keys():if type(secondDict[key]).__name__=='dict':thisDepth=1+getTreeDepth(secondDict[key])else:thisDepth=1if thisDepth>maxDepth:maxDepth=thisDeothreturn maxDepth
绘制树节点
import matplotlib.pyplot as plt##sawtooth是将文本框的类型定义为锯齿形decisionNode=dict(boxstyle="sawtooth",fc="0.8")leafNode=dict(boxstyle="round4",fc="0.8")arrow_args=dict(arrowstyle="<-")def plotNode(nodeTxt,centerPt,parentPt.nodeType):##定义一个绘图区域,绘制一个节点createPlot.ax1.annotate(nodeTxt,xy=parentPt,\xycoords='axes fraction',xytext=centerPt,textcoords='axes fraction',\va="center",bbox=nodeType,arrowprops=arrow_args)def createPlot():fig=plt.figure(1,facecolor='white')fig.clf()##subplot(pos, **kwargs)##pos是子块的索引##frameon表示是否绘制矩形边框createPlot.ax1=plt.subplot(111,frameon=False)##绘制带有标签的节点plotNode('决策节点',(0.5,0.1),(0.1,0.5),decisionNode)plotNode('叶节点',(0.8,0.1),(0.3,0.8),leafNode)
绘制整棵决策树
##在父子节点间填写文本信息def plotMidText(cntrPt,parentPt,txtString):##找到一个中间位置xMid=(parentPt[0]-cntrPt[0])/2.0+cntrPt[0]yMid=(parentPt[1]-cntrPt[1])/2.0+cntrPt[1]createPlot,ax1,text(xMid,yMid,txtString)def plotTree(myTree,parentPt,nodeTxt):##用之前的函数找出叶子节点的数量和深度##以此计算树的宽和高numLeafs=getNumLeafs(myTree)depth=getTreeDepth(myTree)firstStr=myTree.keys()[0]##用plotTree.xOff和plotTree.yOff将绘制的节点放在##相对合适的位置cntrPt=(PlotTree.xOff+(1.0+float(numLeafs))/2.0/plotTree.totalW,\plotTree.yOff)##标记该节点属性值,由此分支向下的节点必包含这一属性plotMidText(cntrPt,parentPt,nodeTxt)##绘制该节点plotNode(firstStr,cntrPt,parentPt,decisionNode)secondDict=myTree[firstStr]##因为树的层数增加了，按比例减小yOff值plotTree.yOff=plotTree.yOff-1.0/plotTree.totalD##绘制子节点，如果子节点是字典，就递归绘制##如果子节点是叶子，就直接绘制一个节点即可for key in secondDict.keys():if type(secondDict[key]).__name__=='dict':plotTree(secondDict[key],cntrPt,str(key))else:plotTree.xOff=plotTree.xOff+1.0/plotTree.totalWplotNode(secondDict[key],(plotTree.xOff,plotTree.yOff),cntrPt,leafNode)plotMidText((plotTree.xOff,plotTree.yOff),cntrPt,str(key))plotTree.yOff=plotTree.yOff+1.0/plotTree.totalD
createPlot应用于决策树的2.0版本
用整棵树的宽和高确定画树的位置
def createPlot(inTree)fig=plt.figures(1,facecolor='white')fig.clf()axprops=dict(xticks=[],yticks=[])createPlot.ax1=plt.subplot(111,frameon=False,**axprops)##计算总宽度，总高度，xOff,yOffplotTree.totalW=float(getNumLeafs(inTree))plotTree.totalD=float(getTreeDepth(inTree))plotTree.xOff=-0.5/plotTree.totalWplotTree.yOff=1.0##设置好参数后调用plotTreeplotTree(inTree,(0.5,1.0),'')plt.show()
使用构建的决策树执行分类
递归比较测试数据和决策树上的数值，递归执行该过程直到进入叶子节点
def classify(inputTree,featLabels,testVec):firstStr=inputTree.keys()[0]secondDict=inputTree[firstStr]##将标签字符串转化为索引featIndex=featLabels.index(firstStr)for key in secondDict.keys():if testVec[featIndex]==key:if type(secondDict[key]).__name__=='dict'##继续递归分类classLbel=classify(secondDict[key],featLabels,testVec)else:##找到合适的分类啦classLabel=secondDict[key]return classLabel
用pickle序列化决策树，存储在磁盘中，下次用就不用再构造啦
##存储def storeTree(inputTree,filename):import picklefw=open(filename,'w')pickle.dump(inputTree,fw)fw.close()##调用def grabTree(filename)import picklefr=open(filename)return pickle.load(fr)
三、朴素贝叶斯
数学基础：条件概率，贝叶斯公式
两个假设：
所有特征之间相互独立每个特征同等重要
p(ci∣w)=p(w0∣ci)...p(wn∣ci)p(ci)p(w)p(c_{i}|w)=\frac{p(w_{0}|c_{i})...p(w_{n}|c_{i})p(c_{i})}{p(w)} p(ci∣w)=p(w)p(w0∣ci)...p(wn∣ci)p(ci)
准备数据：从文本间构建词向量
def loadDataSet():postingList=[['my','dog','has','flea','problems','help','please'],\['maybe','not','take','him','to','dog','park','stupid'],\['my','dalmation','is','so','cute','I','love','him'],\['stop','posting','stupid','worthless','garbage'],\['mr','licks','ate','my','steak','how','to','stop','him'],\['quit','buying','worthless','dog','food','stupid']]classVec=[0,1,0,1,0,1]return postingList,classVecdef createVocabList(dataSet):##创建一个空集合vocabSet=set([])for document in dataSet:##求没有重复元素的并集vocabSet=vocabSet|set(document)return list(vocabSet)##这个函数用来检测vocabList里是否有一段话里面的单词##返回一个01词向量def setOfWords2Vec(vocabList,inputSet):##创建元素全为0的向量returnVec=[0]*len(vocabList)for word in inputSet:if word in vocabList:returnVec[vocabList.index(word)]=1else:print"The word:%s is not in my volcabulary" % wordreturn returnVec
朴素贝叶斯训练器
返回的概率将直接用于分类中
def trainNB0(trainMatrix,trainCategory):numTrainDocs=len(trainMatrix)numWords=len(trainMatrix[0])##求p(c)pAbusive=sum(trainCategory)/float(numTrainDocs)##初始化p(w0|c)和p(w1|c)p0Num=zeros(numWords)p1Num=zeros(numWords)p0Denom=0.0p1Denom=0.0##统计所有训练集中的词汇，一旦在文档中出现##该词对应的0，则pNum0++,该词对应的1,则pNum1++for i in range(numTrainDocs):if trainCategory[i]==1:p1Num+=trainMatrix[i]p1Denom+=sum(trainMatrix[i])p1Vect=p1Num/p1Denomp0Vect=p0Num/p0Denomreturn p1Vect,p0Vect,pAbusive
朴素贝叶斯分类函数
def classifyNB(vec2Classify,p0Vec):listOPosts,listClasses=loadDataSet()myVocabList=createVocabList(listOPosts)##numpy中对应相乘##词汇表的所有词的概率相加，并映射到log函数p1=sum(vec2Classify*p1Vec)+log(pClass1)p0=sum(vec2Classify*p0Vec)+log(pClass0)##分类为骂人句子和非骂人句子if p1>p0:return 1elsereturn 0
测试一下
def testingNB()listOPosts,listClasses=loadDataSet()myVocabList=createVocabList(listOPosts)trainMat=[]for postinDoc in listOPosts:trainMat.append(setOfWords2Vec(myVocabList,postinDoc))p0V,p1V,pAb=trainNB0(array(trainMat),array(listClasses))testEntry=['love','my','dalmation']##利用setOfWords2Vec进行词表到向量的转换thisDoc=array(setOfWords2Vec(myVocabList,testEntry))##进行classify并输出结果print testEntry,'classified as:',classifyNB(thisDoc,p0V,p1V,pAb)
将每个词出现与否作为一个特征，可以被描述为词集模型
对应之前的函数setOfWords2Vec
如果一个词在文档中出现不止一次，意味着会传达更多的信息，不能单纯用0和1衡量。于是引入了词袋模型。
def bagOfWords2VecMN(vocabList,inputSet):returnVec=[0]*len(vocabList)for word in inputSet:if word in vocabList:returnVec[vocabList.index(word)]+=1return returnVec
使用朴素贝叶斯进行交叉验证
def textParse(bigString):import re##把字符串分割成词listOfTokens=re.split(r'\W*,bigString)##返回纯小写的词条,去掉少于两个字符的字符串return [tok.lower() for tok in listOfTokens if len(tok)>2]def spamTest():docList=[]classList=[]fullText=[]for i in range(1,26):##打开要切分的文本wordList=textParse(open('email/spam/%d.txt'%i).read())docList.append(wordList)fullText.extend(wordList)##这些词被分类为abusiveclassList.append(1)wordList=textParse(open('email/ham/%d.txt'%i).read())docList.append(wordList)fullList.extend(wordList)classList.append(0)##创建词条模型vocabList=createVocabList(docList)trainingSet=range(50)testSet=[]##接下来随机构造训练集,留存交叉验证for i in range(10):##trainingSet是0到49的整数列表randIndex=int(random.uniform(0,len(trainingSet)))testSet.append(trainingSet[randIndex])del(trainingSet[randIndex])trainMat=[]trainClasses=[]for docIndex in trainingSet:##对每封邮件构建词向量trainMat.append(setOfWords2Vec(vocabList,docList[docIndex]))trainClasses.append(classList[docIndex])##计算概率p0V,p1V,pSpam=trainNB0(array(trainMat),array(trainClasses))errorCount=0##遍历测试集，对每封邮件进行验证for docIndex in testSet:wordVector=setOfWords2Vec(vocabList,docList[docIndex])if classifyNB(array(wordVector),p0V,p1Vm,pSpam)!=classList:errorCount+=1print'the error rate is',float(errorCount)/len(testSet)
使用RSS源
##计算最经常出现的30个单词def calcMostFreq(vocabList,fullText)import operatorfreqDict=[]for token in vocabList:freqDict[token]=fullText.count(token)sortedFreq=sorted(freqDict.iteritems(),key=operator.itemgetter(1),reverse=True)return sortedFreq[:30]##运用RSS源的spamtextdef localWords(feed1,feed0):import feedparserdocList=[]classList=[]fullText=[]minLen=min(len(feed1['entries'][i]['summary']),len(feed0['entries'][i]['summary']))for i in range(minLen):wordList=textParse(feed0['entries'][i]['summary'])docList.append(wordList)fullText.extend(wordList)classList.append(0)wordList=textParse(feed1['entries'][i]['summary'])docList.append(wordList)fullText.extend(wordList)classList.append(1)##创建vocabListvocabList=createVocabList(docList)top30Words=calcMostFreq(vocabList,fullText)##去掉频率最高的30个单词##（出于实际应用考虑，因为广告词中重复的词是没用的）for pairW in top30Words:if pairW[0] in vocabList:vocabList.remove(pairW[0])trainingSet=range(2*minLen);##随机选出训练集for i in range(20):##uniform(x,y):随机生成一个在x,y之间的实数randIndex=int(random.uniform(0,len(trainingSet)))testSet.append(trainingSet[randIndex])del(trainingSet[randIndex])trainMat=[]trainClasses=[]for docIndex in trainingSet:trainMat.append(bagOfWords2VecMN(vocabList,docList[docIndex]))trainClasses.append(classList[docIndex])##利用trainNB0求classify所需的概率p0V,p1V,pSpam=trainNB0(array(trainMat),array(trainClasses))errorCount=0for docIndex in testSet:wordVector=bagOfWords2VecMN(vocabList,docList[docIndex])if classifyNB(array(wordVector),p0V,p1V,pSpam)!=classList[docIndex]:errorCount+=1print'the error rate is:',float(errorCount)/len(testSet)return vocabList,p0V,p1V
四、Logistic回归
回归：寻找最佳的拟合系数
单位阶跃函数不能满足连续性。但如果我们把每个特征乘上一个回归系数，然后把所有结果值相加，总和代入sigmoid函数中。
大于0.5分为1类，小于0.5分为0类
sigmoid函数式和图像
梯度上升
z=wtxz=w^{t}x z=wtx
一般用梯度上升法求局部最大值，梯度下降法求局部最小值
在每一点处朝着梯度方向迭代（梯度方向就是偏导数方向）
w:=w+α∇w+f(w)w:=w+\alpha\nabla_{w}+f(w)w:=w+α∇w+f(w)
α\alphaα是向目标移动的步长
直到某个停止条件为止
def loadDataSet():dataMat=[]labelMat=[]fr=open('testSet.txt')for line in fr.readlines():##数据中，每行前两个是X1和X2，第三个是数据标签lineArr=line.strip().split()dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])labelMat.append(int(lineArr[2]))return dataMat,labelMat##计算sigmoid函数def sigmoid(inX):return 1.0/(1+exp(-inX))def gradAscent(dataMatIn,classLabels):##转化为numpy矩阵数据类型dataMatrix=mat(dataMatIn)##转置为列向量labelMat=mat(classLabels).transpose()m,n=shape(dataMatrix)alpha=0.001maxCycles=500weights=ones(n,1)for k in range(maxCycles):h=sigmoid(dataMatrix*weights)##用sigmoid函数很好计算errorerror=(labelMat-h)##梯度上升迭代weights=weights+alpha*dataMatrix.transpose()*errorreturn weights
matplotlib画出最佳拟合直线函数w2y=w1x+b
def plotBestFit(weights):import matplotlib.pyplot as pltdataMat,labelMat=loadDataSet()dataArr=array(dataMat)n=shape(dataArr)[0]xcord1=[]ycord1=[]xcord2=[]ycord2=[]for i in range(n):if int(labelMat[i])==1:xcord1.append(dataArr[i,1])ycord1.append(dataArr[i,2])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')x=arrange(-3.0,3.0,0.1)##反解出x2##0是两个类别的分界处##0=w0x0+w1x1+w2x2##这里w0=1（因为是一条直线）y=(-weights[0]-weights[1]*x)/weights[2]ax.plot(x,y)plt.xlabel('X1')plt.ylabel('X2')plt.show()
随机梯度上升算法
随机取数据集中的数据进行梯度上升算法，降低时间复杂度，是一个在线算法。
def stocGradAscent1(dataMatrix,classLabels):m,n=shape(dataMatrix)weights=ones(n)##j是迭代次数，i是样本点的下标for j in range(numIter):dataIndex=range(m)for i in range(m):##alpha每次迭代都会减小##有利于最后迭代次数降低，稳定在一个值，收敛的更快alpha=4/(1.0+j+i)+0.01##随机选择一个点更新参数randIndex=int(random.uniform(0,len(dataIndex))h=sigmoid(sum(dataMatrix[randIndex]*weights))error=classLabels[randIndex]-hweights=weights+alpha*error*dataMatrix[randIndex]##将用过的index删掉del(dataIndex[randIndex])return weights
示例：预测病马的死亡率
处理数据中缺失值的方法：
使用可用特征的均值填补缺失值使用特殊值填补缺失值，如-1忽略有缺失值的样本使用相似样本的均值填补缺失值使用另外的机器学习算法预测缺失值
这里运用logistic回归，因此补0就好(sigmoid(0)=0.5)，不会影响到分类函数
另外，分类的时候把测试集上的每个特征和回归系数相乘再相加，输入sigmoid函数，即可得到分类结果。
##分类函数def classifyVector(inX,weights):prob=sigmoid(sum(inX*weights))if prob>0.5:return 1.0else:return 0.0def colicTest():frTrain=open('horseColicTraining.txt')frTest=open('horseColicTest.txt')trainingSet=[]trainingLabels=[]##对训练集的处理for line in frTrain.readlines():currLine=line.strip().split('\t')lineArr=[]for i in range(21):lineArr.append(float(currLine[i]))trainingSet.append(lineArr)trainingLabels.append(float(currLine[21]))trainWeights=stocGradienAscent1(array(trainingSet),trainingLabels,500)##对测试集的处理errorCount=0numTestvec=0.0for line in frTest.readlines():numTestVec+=1.0currLine=line.strip().split('\t')lineArr=[]for i in range(21):lineArr.append(float(currLine[i]))if int(classifyVector(array(lineArr),trainWeights))!=int(currLine[21]):errorCount+=1errorRate=(float(errorCount)/numTestVec)print"the error rate of this test is %f" % errorRatereturn errorRatedef multiTest():numTests=10errorSum=0.0for k in range(numTests):errorSum+=colicTest()print"after %d iterations the average error rate is:%f" %(numTests,errorSum/float(numTests))
五、支持向量机
基本推导：
是个泛化错误率低的二分类问题
支持向量就是离分隔超平面最近的点。要最大化支持向量到分隔面的距离。
argmax{min(label(wTx+b))1∣∣w∣∣}argmax \{min(label(w^{T}x+b))\frac{1}{||w||}\}argmax{min(label(wTx+b))∣∣w∣∣1}
其中label是分类标签（-1、1）
约束条件：（共有i个）
label∗(wTx+b)≥1.0label*(w^{T}x+b)\geq1.0label∗(wTx+b)≥1.0
下面引入拉格朗日乘数法求条件极值，解决一个凸优化问题
对w和b求偏导代入，得到：（这个不用讲了吧，高等数学）
L=∑i=1mα−12∑i,j=1mlabelilabeljαiαj<xixj>L=\sum_{i=1}^{m}\alpha-\frac{1}{2}\sum_{i,j=1}^{m}label^{i}label^{j}\alpha^{i}\alpha{j}<x^{i}x^{j}>L=i=1∑mα−21i,j=1∑mlabelilabeljαiαj<xixj>
但是我们的约束条件是不等号，于是引入对偶问题，这里是一个强对偶关系
一般来说，minmax大于等于maxmin，只有在满足kkt条件的时候才会有两者相等。
条件为：
这篇博客讲的比较清晰：
我是对偶问题推导
minmax[∑i=1mα−12∑i,j=1mlabelilabeljαiαj<xixj>]minmax[\sum_{i=1}^{m}\alpha-\frac{1}{2}\sum_{i,j=1}^{m}label^{i}label^{j}\alpha^{i}\alpha{j}<x^{i}x^{j}>]minmax[i=1∑mα−21i,j=1∑mlabelilabeljαiαj<xixj>]
st.α≥0,∑i=1mαilabelist.\alpha\geq0,\sum^{m}_{i=1}\alpha_{i}label^{i}st.α≥0,i=1∑mαilabeli
引入松弛变量c，允许有些分类在错误的一侧
一旦求出了所有的α\alphaα，分隔超平面就很容易能表达出来，SVM的求解目的就在于求出满足条件的α\alphaα
书上没有详细证明，可以参考这个教程：
我是SVM算法推导
SMO算法
上面有个遗留问题，就是在这样的带约束条件的凸优化问题中，怎么求解α\alphaα？
这里的α\alphaα是一个向量，里面有m个α1,...,αm\alpha_{1},...,\alpha_{m}α1,...,αm
找到α\alphaα以后，w和b都能确定，整个超平面也就能够确定了。
在SMO算法中每次循环选择两个α\alphaα进行优化
SMO算法中几个杂七杂八的函数：
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):j=iwhile(j==i)j=int(random.uniform(0,m))return j##对太大太小的数值进行调整```pythondef clipAlpha(aj,H,L):if aj>H:aj=Hif L>aj:aj=Lreturn aj
简化版SMO算法：
代码中误差和alpha、eta计算公式的推导
def smoSimple(dataMatIn,classLabels,C,toler,maxIter):dataMatrix=mat(dataMatIn)labelMat=mat(classLabels).transpose()b=0m,n=shape(dataMatrix)alphas=mat(zeros(m,1))iter=0while(iter<maxIter):##记录alpha是否被优化过alphaParisChanged=0for i in range(m)：##对于数据集中的每个数据向量##fXi是分类的类别，Ei是误差fXi=float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T))+b##.T是转置操作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)fXj=float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T))+bEj=fXj-float(labelMat[j])##新旧值的比较,将alpha[j]调整到0与C之间alphaIold=alphas[i].copy()alphaJold=alphas[j].copy()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##eta是alpha[j]的最佳修改量eta=2.0*dataMatrix[i:]*dataMatrix[j:].T-dataMatrix[i,:]*dataMatrix[i,:].T-dataMatrix[j,:]*dataMatrix[j,:].Tif eta>=0:print"eta>=0"continue##如果eta==0计算新的alpha[j]alphas[j]-=labelMat[j]*(Ei-Ej)/eta##限定一下alphas[j]的范围alphas[j]=clipAlpha(alphas[j],H,L)if(abs(alpha[j]-alphaJold)<0.00001):print"j is not moving enough"continue##更新a[i]的值##a[j]减小a[i]就增大，数值相同方向相反alpha[i]+=labelMat[j]*labelMat[i]*(alphas[j]-alphaJold)*dataMatrix[j:]##为a[i]设置常数项b1=b-Ei-labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i:].T-labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i:]*dataMatrix[j:].T##为a[j]设置常数项b2=b-Ej-labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i:].T-labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i:]*dataMatrix[j:].T##统一常数项if(0<alphas[i])and(C>alphas[i]):b=b1elif(0<alphasp[j])and(C>alphas[j]):b=b2else:b=(b1+b2)/2.0alphaPairsChanged+=1print"iter:%d i:%d pairschanged:%d"%(iter,i,alphaPairsChanged)##检查alpha是否更新，更新了则重置iter为0##因为要在所有数据上遍历maxiter次，发现没有可更新的alpha后##整个程序才会完美结束if(alphaPairsChanged==0):iter+=1else:iter=0print"iteration number:%d"%iterreturn b,alphas
完整Platt SMO算法加速优化
优化点：
建立非边界α\alphaα值的缓存用于保存误差值，并且从中选择使步长(误差之差)最大的α\alphaα
完整版中杂七杂八的函数
class optStruct:##误差缓存def __init__(self,dataMatIn,classLabels,C,toler):self.X=dataMatInself.labelMat=classLabelsself.C=Cself.tol=tolerself.m=shape(dataMatIn)[0]self.alpha=mat(zeros(self,m,1))self.b=0self.eCache=mat(zeros(self.m,2))def calcEk(oS,k):##将算误差值的函数单独拎出来fXk=float(multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k:].T))+oS.bEk=fXk-float(oS.labelMat[k])return Ekdef selectJ(i,oS,Ei):maxK=-1maxDeltaE=0 Ej=0oS.eCache[i]=[i,Ei]##选择具有最大步长的jdef selectJ(i,oS,Ei):maxK=-1maxDeltaE=0Ej=0oS.eCache[i]=[1,Ei]##非0 E值所对应的alpha值列表validEcacheList=nonzero(oS.eCache[:,0].A)[0]if(len(validEcacheList))>1:for k in validEcacheList:if k==i:continueEk=calcEk(oS,k)deltaE=abs(Ei-Ek)##选取改变量最大的值if(deltaE>maxDeltaE):maxK=kmaxDeltaE=deltaEEj=Ekreturn maxK,Ejelse:j=selectJrand(i,oS.m)Ej=calcEk(oS,j)return j,Ej##计算误差值并存入缓存，对alpha进行优化的时候会用到def updateEk(oS,k):Ek=calcEk(oS,k)oS.eCache[k]=[1,Ek]
完整版中内循环(用缓存存储误差,和之前的内循环算法没有差别)
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)alphaIold=oS.alphas[i].copy()alphaJold=oS.alpha[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 0eta=2.0*oS.X[i,:]*oS.X[j,:].T-oS.X[i:]*oS.X[i:].T-oS.X[j,:]*oS.X[j,:].Tif eta>=0:print"eta>=0"return 0oS.alphas[j]-=oS.labelMat[j]*(Ei-Ej)/etaoS.alphas[j]=clipAlpha(oS.alphas[j],H,L)##更新误差缓存updateEK(oS,j)if(abs(oS.alpha[j]-alphaJold)<0.00001):print"j not moving enough"return 0oS.alphas[i]+=oS.labelMat[j]*oS.labelMat[i]*(alphaJold-oS.alphas[j])##更新误差缓存updateEK(os,i)##下面和简单SMO都一样，只不过用了optStruct类，懒得打了b1=。。。。
完整版SMO外循环代码
def smoP(dataMatIn,classLabels,C,toler,maxIter,kTup=('lin',0)):oS=optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler)iter=0entireSet=TruealphaPairsChange=0##循环条件：没到最大迭代次数##在所有数据集上扫描/alpha的值已经被改变过while(iter<maxIter)and((alphaPairsChanged>0)or(entireSet)):alphaPairsChanged=0if entireSet:##遍历所有值for i in range(oS,m):##用内循环选alpha[j]alphaPairsChanged+=innerL(i,oS)print"fullSet,iter:%d i:%d,pairs changed %d" % (iter,i,alphaPairsChanged)iter+=1else:##寻找不在边界的alphanonBoundIs=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+=1if entireSet:entireSet=Falseelif(alphaPairsChanged==0):entireSet=Trueprint:"iteration number:%d" % iterreturn oS.b,oS.alphas
计算w，就可以得到我们的分隔超平面
def calcWs(alphas,dataArr,classLabels):##将数组用mat转化为矩阵并转置X=mat(dataArr)labelMat=mat(classLabels).tanspose()m,n=shape(X)w=zeros((n,1))for i in range(m):##非零alpha对应的都是支持向量##所以最后乘进去的只有支持向量w+=multiply(alphas[i]*labelMat[i],X[i,:].T)return w
核函数
核函数完成从一个特征空间到另一个特征空间的映射
在SVM中我们可以把一个内积运算替换成核函数
SVM中常用径向基核函数
k(x,y)=exp(−∣∣x−y∣∣22σ2)k(x,y)=exp(\frac{-||x-y||^{2}}{2\sigma^{2}})k(x,y)=exp(2σ2−∣∣x−y∣∣2)
σ\sigmaσ是个函数值跌落到0的速度参数
核转换函数：
def kernelTrans(X,A,kTup):##元祖kTup给出核函数的信息##kTup[0]是描述核函数类型的一个字符串##kTup[1]是deltam,n=shape(X)K=mat(zeros(m,1))if KTup[0]=='lin':K=X*A.Telif kTup[0]=='rbf':for j in range(m):deltaRow=X[j,:]-AK[j]=deltaRow*deltaRow.TK=exp(K/(-1*kTup[1]**2))else:raise NameError('Houston We Have a Problem That Kernel is not recognized')return K
对原本的optStruct类进行修改
class optStruct:##误差缓存def __init__(self,dataMatIn,classLabels,C,toler,kTup):self.X=dataMatInself.labelMat=classLabelsself.C=Cself.tol=tolerself.m=shape(dataMatIn)[0]self.alpha=mat(zeros(self,m,1))self.b=0self.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)
在训练集中，用K值改变b的结果
在测试集中，确定sigma的大小，并用该核函数构建一个分类器
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()##找出alpha值非0的支持向量和他的类别标签svInd=nonzero(alphas.A>0)[0]sVs=datMat[svInd]labelSV=labelMar[svInd]print"there are %d Support Vectors"% shape(sVs)[0]m,n=shape(datMat)##对于训练数据集errorCount=0for i in range(m):##得到转换为核函数的数据kernelEval=kernelTrans(sVs,datMat[i,:],('rbf',k1))##核函数与alpha和类别标签求积predict=kernelEval.T*multiply(labelSV,alpha[svInd])+bif sign(predict)!=sign(labelArr[i]):errorCount+=1print"the training error rate id:%f"%(float(errorCount)/m)errorCount=0datMat=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])+bif sign(predict)!=sign(labelArr[i]):errorCount+=1print"the test error rate is:%f"%(float(errorCount)/m)
六、Adaboost
集成方法是对其他算法进行组合的一种方式
bagging：数据重新抽样构建分类器，如随机森林（多个决策树同时分类，看哪种结果多，详情见别人的随机森林博客）
boosting：boosting中的分类器权重并不相等，每个权重代表的是其对应分类器在上一轮迭代中的成功度
一开始，所有权重都初始化为相等的值
α\alphaα是根据每个弱分类器的错误率计算的
α=12ln(1−ϵϵ)\alpha=\frac{1}{2}ln(\frac{1-\epsilon}{\epsilon})α=21ln(ϵ1−ϵ)
算出α\alphaα后，对权重向量D进行更新。
如果某个样本被正确分类
D(i+1)i=Di(t)e−αSum(D)D^{(i+1)_{i}}=\frac{D^{(t)}_{i}e^{-\alpha}}{Sum(D)}D(i+1)i=Sum(D)Di(t)e−α
如果未被正确分类就是eαe^{\alpha}eα
不断迭代直到训练错误率到达指定值
单层决策树生成函数
构建一个弱分类器
def stumpClassify(dataMatrix,dimen,threshVal,threshIneq):retArray=ones((shape(dataMatrix[0],1)))##将所有不满足要求的设置为-1类if threshIneq=='lt':retArray[dataMatrix[:，dimen]<=threshVal]=-1.0else:retArray[dataMatrix[:,dimen]>threshVal]=-1.0return retArray##找到数据集上最佳的单层决策树def buildStump(dataArr,classLabels,D):dataMatrix=mat(dataArr)labelMat=mat(classLabels).Tm,n=shape(dataMatrix)numSteps=10.0bestStump={}bestClasEst=mat(zeros((m,1)))minError=inf##遍历数据集的每一个特征for i in range(n):rangeMin=dataMatrix[:,i].min()rangeMax=dataMatrix[:,i].max()##用最大最小值计算步长stepSize=(rangeMax-rangeMin)/numSteps##对于步长中的每个值for j in range(-1,int(numSteps)+1):##考虑大于等于和小于等于两种情况for inequal in ['lt','gt']:threshVal=(rangeMin+float(j)*stepSize)##计算分类预测结果predictedVals=stumpClassify(dataMatrix,i,threshVal,inequal)errArr=mat(ones((m,1)))##分类对了就为0，分类错了就为1errArr[predictedVals==labelMat]=0##计算加权错误率weightedError=D.T*errArr##找到加权错误率最小的if weightedError<minError:minError=weightedErrorbestClasEst=predictedVals.copy()bestStump['dim']=ibestStump['thresh']==threshValbestStump['ineq']=inequalreturn bestStump,minError,bestClasEst
使用多个弱分类器构建Adaboost代码
任何分类器都可以作为基础分类器，本书讲到的任何一个算法都行
def adaBoostTrainDS(dataArr,classLabels,numIt=40):weakClassArr=[]m=shape(dataArr)[0]D=mat(ones((m,1))/m)aggClassEst=mat(zeros((m,1)))for i in range(numIt):bestStump,error,classEst=buildStump(dataArr,classLabels,D)print"D:",D.T##计算alpha值alpha=float(0.5*log((1.0-error)/max(error,1e-16)))##将alpha值加入bestStump字典中bestStump['alpha']=alpha##字典添加到弱分类器列表中weakClassArr.append(bestStump)print"ClassEst:",ClassEst.T##计算新权重Dexpon=multiply(-1*alpha*mat(classLabels).T,classEst)D=multiply(D,exp(expon))D=D/D.sum()##aggClassEst是每个数据点的类别估计的累积值##可以通过aggClassEst的符号来了解总类别aggClassEst+=alpha*classEstprint "aggClassEst:",aggClassEst.T##计算累积错误率##sign()正数是1，负数是-1,0是0aggErrors=multiply(sign(aggClassEst)!=mat(classLabels).T,ones((m,1)))errorRate=aggErrors.sum()/mprint"totalerror:",errorRate,"\n"##分类错误率为0，结束循环if errorRate==0.0:breakreturn weakClassArr
Adaboost分类函数
def adaClassify(datToClass,classifierArr):dataMatrix=mat(datToClass)m=shape(dataMatrix)[0]aggClassEst=mat(zeros((m,1)))##遍历所有弱分类器for i in range(len(classifierArr)):##用stumpClassify对每个分类器估计一个值classEst=stumpClassify(dataMatrix,classifierArr[i]['dim'],classifierArr[i]['thresh'],classifierArr[i]['ineq'])##aggClassEst:每个数据点类别估计的累积值aggClassEst+=classifierArr[i]['alpha']*classEstprint aggClassEstreturn sign(aggClassEst)
非均衡分类问题
有时候每个分类的代价不同，比如合法邮件案例中，宁愿误判不合法的邮件为合法邮件，也不要漏掉合法邮件造成损失。
代价敏感学习/对训练数据进行改造（欠抽样、过抽样）
正确率，召回率，ROC曲线
被西瓜书第二章折磨过的懂的都懂
ROC曲线下面的面积就是AUC，给出分类器的平均性能值，完美分类器的AUC是1.0，随机分类器的AUC是0.5
用代码画出ROC曲线
分类样例按照预测强度排序，依次作为正反例分割点，计算ROC曲线中对应的值
def plotROC(preStrengths,classLabels):import matplotlib.pyplot as pltcur=(1.0,1.0)ySum=0.0numPosClas=sum(array(classLabels)==1.0)##计算x轴和y轴上的步长yStep=1/float(numPosClas)xStep=1/float(len(classLabels)-numPosClas)##将分类器的预测强度排序sortedIndicies=predStrengths.argsort()fig=plt.figure()fig.clf()ax=plt.subplot(111)for index in sortedIndicies.tolist()[0]if classLabels[index]==1.0:delX=0delY=yStepelse:delX=xStepdelY=0##ySum计算AUC的值,对小矩形累加ySum+=cur[1]ax.plot([cur[0],cur[0]-delX],[cur[1],cur[1]-delY)cur=(cur[0]-delX,cur[1]-delY)ax.plot([0,1],[0,1],'b--')plt.xlabel('False Positive Rate')plt.ylabel('True Positive Rate')plt.title('ROC curve for AdaBoost')plt.show()##xStep是小矩形的宽print"the Area Under the Curve is:",ySum*xStep
