朴素贝叶斯

发表于 2018-08-28 | 分类于机器学习专题

字数统计: 3.5k | 阅读时长 ≈ 15

注：这篇文章原文出处放在了下面的连接，原文已经写得非常清晰明了，这里有适当的修改和增加内容，写在这里是为了更好地查阅巩固。

戳我原文

朴素贝叶斯（Naive Bayes）是基于贝叶斯定理与特征条件假设的分类方法。

对于给定的训练数据集，首先基于特征条件独立假设学习输入、输出的联合分布；然后基于此模型，对给定的输入x，利用贝叶斯定理求出后验概率最大的输出y。

朴素贝叶斯实现简单，学习与预测的效率都很高，是一种常用的方法。

朴素贝叶斯的学习与分类

贝叶斯定理

先看什么是条件概率

P(A|B)表示事件B已经发生的前提下，事件A发生的概率，叫做事件B发生下事件A的条件概率。其基本求解公式为

$P\left(A|B\right)=\frac{P\left(AB\right)}{P\left(B\right)}$

贝叶斯定理便是基于条件概率，通过P(A|B)来求P(B|A)：

$P\left(B|A\right)=\frac{P\left(A|B\right)·P\left(B\right)}{P\left(A\right)}$

顺便提一下，上式中的分母，可以根据全概率公式分解为：

$P\left(A\right)=\sum_{i=1}^n{P\left(B_i\right)P\left(A|B_i\right)}$

特征条件独立假设

这一部分开始朴素贝叶斯的理论推导，从中你会深刻地理解什么是特征条件独立假设。

给定训练数据集(X,Y)，其中每个样本X都包括nn维特征，即$x=(x_1,x_2,···,x_n)$，类标记集合含有K种类别，即$y=(y_1,y_2,···,y_k)$

如果现在来了一个新样本x我们要怎么判断它的类别?从概率的角度来看，这个问题就是给定x，它属于哪个类别的概率更大。那么问题就转化为求解$P(y_1|x),P(y_2|x),P(y_k|x)$中最大的那个，即求后验概率最大的输出：$arg\underset{y_k}{\max}P\left(y_k|x\right)$ 。那$P(y_k|x)$怎么求解？答案就是贝叶斯定理

$P\left(y_k|x\right)=\frac{P\left(x|y_k\right)·P\left(y_k\right)}{P\left(x\right)}$

根据全概率公式，可以进一步分解上式中的分母：

$P\left(y_k|x\right)=\frac{P\left(x|y_k\right)·P\left(y_k\right)}{\sum_{i=1}^n{P\left(x|y_k\right)P\left(y_k\right)}} （公式1）$

先不管分母，分子中的$P(y_k)$是先验概率，根据训练集就可以简单地计算出来，而条件概率$P(x|y_k)=P(x_1,x_2,···,x_n|y_k)$它的参数规模是指数数量级别的，假设第$i$维特征$x_i$可取值的个数有$S_i$个，类别取值个数为k个，那么参数个数为$k\prod_{j=1}^n{S_j}$

这显然是不可行的。针对这个问题，朴素贝叶斯算法对条件概率分布做了独立性的假设，通俗地讲就是说假设各个维度的特征 $x_1,x_2,···,x_n$互相独立，由于这是一个较强的假设，朴素贝叶斯算法也因此得名。在这个假设的前提上，条件概率可以转化为：

$P\left(x|y_i\right)=P\left(x_1,x_2,···,x_n|y_i\right)=\prod_{i=1}^n{P\left(x_i|y_i\right)} （公式2）$

这样参数规模就降到了$\sum_{i=1}^n{S_ik}$

以上就是针对条件概率所作出的特征条件独立性假设，至此，先验概率$P(y_k)$和条件概率$P(x|y_k)$的求解问题就都解决了，那么我们是不是可以求解我们所需要的后验概率$P(y_k|x)$了？

答案是肯定的。我们继续上面关于$P(y_k|x)$的推导，将公式2代入公式1中得到：

$P\left(y_k|x\right)=\frac{P\left(y_k\right)\prod_{i=1}^n{P\left(x_i|y_k\right)}}{\sum_k{P\left(y_k\right)\prod_{i=1}^n{P\left(x_i|y_k\right)}}}$

于是朴素贝叶斯分类器可表示为：

$f\left(x\right)=arg\underset{y_k}{\max}P\left(y_k|x\right)=arg\underset{y_k}{\max}\frac{P\left(y_k\right)\prod_{i=1}^n{P\left(x_i|y_k\right)}}{\sum_k{P\left(y_k\right)\prod_{i=1}^n{P\left(x_i|y_k\right)}}}$

因为对于所有的$y_k$，上式中的分母的值都是一样的（为什么？注意到全加符号就容易理解了），所以可以忽略分母部分，朴素贝叶斯分裂期最终表示为：

$f\left(x\right)=arg\underset{y_k}{\max}P\left(y_k\right)\prod_{i=1}^n{P\left(x_i|y_k\right)}$

朴素贝叶斯法的参数估计

极大似然估计

根据上述，可知朴素贝叶斯要学习的东西就是$P(Y=c_k)$和$P(X^{j}=a_{jl}|Y=c_k)$可以应用极大似然估计法估计相应的概率（简单讲，就是用样本来推断模型的参数，或者说是使得似然函数最大的参数）

先验概率$P(Y=c_k)$的极大似然估计是

$P\left(Y=c_k\right)=\frac{\sum_{i=1}^N{I\left(y_i=c_k\right)}}{N},\,\,k=1,2,···,K$

也就是用样本中$c_k$的出现次数除以样本容量。

设第$j$个特征$x^{(j)}$可能取值的集合为${a_{j1},a_{j2},···,a_{jl}}$，条件概率 $P(X^{j}=a_{jl}|Y=c_k)$的极大似然估计是：

$P\left(X^{\left(j\right)}=a_{jl}|Y=c_k\right)=\frac{\sum_{i=1}^N{I\left(x_{i}^{\left(j\right)}=a_{jl},y_{i=}c_k\right)}}{\sum_{i=1}^N{I\left(y_i=c_k\right)}}$

式中，$x_i^j$是第$i$个样本的第$j$个特征。

一个例题

例题如下

贝叶斯估计

极大似然估计有一个隐患，假设训练数据中没有出现某种参数与类别的组合怎么办？比如上例中当Y=1对应的$X^{(1)}$的取值只有1和2。这样可能会出现所要估计的概率值为0的情况，但是这不代表真实数据中就没有这样的组合。这时会影响到后验概率的计算结果，使分类产生偏差。解决办法是贝叶斯估计。

条件概率的贝叶斯估计：

$P_{\lambda}\left(X^{\left(j\right)}=a_{jl}\parallel Y=c_k\right)=\frac{\sum_{i=1}^N{I\left(x_{i}^{\left(j\right)}=a_{jl},y_{i=}c_k\right)}+\lambda}{\sum_{i=1}^N{I\left(y_i=c_k\right)}+S_j\lambda}$

其中$\lambda≥0$，$S_j$表示第$x_j$个特征可能取值的个数。分子和分母分别比极大似然估计多了一点东西，其意义为在随机变量各个取值的频数上赋予一个正数$λ≥0$。当$λ=0$时就是极大似然估计。常取$λ=1$，这时称为拉普拉斯平滑。

先验概率的贝叶斯估计：（K为类别个数）

$P_{\lambda}\left(Y=c_k\right)=\frac{\sum_{i=1}^N{I\left(y_i=c_k\right)}+\lambda}{N+K\lambda}$

例题如下：套入公式计算即可

上诉说的是$X_j$为普通离散型的情况。如果$X_j$是稀疏的离散值，即各个特征的出现概率很低，那么可以假设$X_j$服从伯努利分布，即特征$X_j$出现记为1，不出现记为0。即我们只关注$X_j$是否出现，不关注$X_j$出现的次数，这样得到的$P\left(X^{\left(j\right)}=a_{jl}|Y=c_k\right)$是在是在样本类别$c_k$中$a_{jl}$出现的频率，公式如下所示。其中$a_{jl}$取值为0和1。

$P(X_j=a_{jl}|Y=C_k)=P(X_j|Y=C_k)a_{jl}+(1-P(X_j|Y=C_k))(1-a_{jl})$

如果$X_j$是连续值，那么假设$X_j$的先验概率为高斯分布(正态分布)，这样假设$P\left(X^{\left(j\right)}=a_{jl}|Y=c_k\right)$的概率分布公式如下所示。其中$\mu_k$和$\sigma_{k}^2$是正态分布的期望和方差，$\mu_k$为样本类别$C_k$中，所有$X_j$的平均值，$\sigma_{k}^2$为样本类别$C_k$中，所有$X_j$的方差，$\mu_k$和$\sigma_{k}^2$可以通过极大似然估计求得。

$P(X_j=a_{jl}|Y=C_k)=\frac{1}{\sqrt{2 \pi \sigma^2}}exp(-\frac{(a_{jl}- \mu_k)^2}{2\sigma_{k}^{2}})$

python代码实现

朴素贝叶斯文档分类

# -*- coding: utf-8 -*-
"""
Created on 下午5:28 22 03 2017
bayes algorithm: classify a words as good or bad   [text classify]
@author: plushunter
"""
from numpy import *
class Naive_Bayes:
    def __init__(self):
        self._creteria = "NB"
    #创建不重复词集
    def _creatVocabList(self,dataSet):
        vocabSet = set([])  # 创建一个空的SET
        for document in dataSet:
            vocabSet = vocabSet | set(document)  # 并集
        return list(vocabSet)  # 返回不重复词表（SET的特性）
    #文档词集向量模型
    def _setOfWordToVec(self,vocabList, inputSet):
        """
        功能:给定一行词向量inputSet，将其映射至词库向量vocabList，出现则标记为1，否则标记为0.
        """
        returnVec = [0] * len(vocabList)
        for word in inputSet:
            if word in vocabList:
                returnVec[vocabList.index(word)] = 1
        return returnVec
    #文档词袋模型
    def _bagOfsetOfWordToVec(self,vocabList, inputSet):
        """
        功能：对每行词使用第二种统计策略，统计单个词的个数，然后映射到此库中
        输出：一个n维向量，n为词库的长度，每个取值为单词出现的次数
        """
        returnVec = [0] * len(vocabList)
        for word in inputSet:
            if word in vocabList:
                returnVec[vocabList.index(word)] += 1 #更新此处代码
        return returnVec
    def _trainNB0(self,trainMatrix, trainCategory):
        """
        输入：训练词矩阵trainMatrix与类别标签trainCategory,格式为Numpy矩阵格式
        功能：计算条件概率p0Vect、p1Vect和类标签概率pAbusive
        """
        numTrainDocs = len(trainMatrix)#样本个数
        numWords = len(trainMatrix[0])#特征个数，此处为词库长度
        pAbusive = sum(trainCategory) / float(numTrainDocs)#计算负样本出现概率（先验概率）
        p0Num = ones(numWords)#初始词的出现次数为1，以防条件概率为0，影响结果
        p1Num = ones(numWords)#同上
        p0Denom = 2.0#类标记为2，使用拉普拉斯平滑法,
        p1Denom = 2.0
        #按类标记进行聚合各个词向量
        for i in range(numTrainDocs):
            if trainCategory[i] == 0:
                p0Num += trainMatrix[i]
                p0Denom += sum(trainMatrix[i])
            else:
                p1Num += trainMatrix[i]
                p1Denom += sum(trainMatrix[i])
        p1Vect = log(p1Num / p1Denom)#计算给定类标记下，词库中出现某个单词的概率
        p0Vect = log(p0Num / p0Denom)#取log对数，防止条件概率乘积过小而发生下溢
        return p0Vect, p1Vect, pAbusive
    def _classifyNB(self,vec2Classify, p0Vec, p1Vec, pClass1):
        """
        该算法包含四个输入:
        vec2Classify表示待分类的样本在词库中的映射集合，
        p0Vec表示条件概率P(wi|c=0)P(wi|c=0)，
        p1Vec表示条件概率P(wi|c=1)P(wi|c=1)，
        pClass1表示类标签为1时的概率P(c=1)P(c=1)。
        p1=ln[p(w1|c=1)p(w2|c=1)…p(wn|c=1)p(c=1)]
        p0=ln[p(w1|c=0)p(w2|c=0)…p(wn|c=0)p(c=0)]
        log取对数为防止向下溢出
        功能:使用朴素贝叶斯进行分类,返回结果为0/1
        """
        p1 = sum(vec2Classify * p1Vec) + log(pClass1)
        p0 = sum(vec2Classify * p0Vec) + log(1 - pClass1)
        if p1 > p0:
            return 1
        else:
            return 0
    #test
    def testingNB(self,testSample):
        "step1：加载数据集与类标号"
        listOPosts, listClasses = loadDataSet()
        "step2：创建词库"
        vocabList = self._creatVocabList(listOPosts)
        "step3：计算每个样本在词库中出现的情况"
        trainMat = []
        for postinDoc in listOPosts:
            trainMat.append(self._bagOfsetOfWordToVec(vocabList, postinDoc))
        p0V, p1V, pAb = self._trainNB0(trainMat, listClasses)
        "step4：测试"
        thisDoc = array(self._bagOfsetOfWordToVec(vocabList, testSample))
        result=self._classifyNB(thisDoc, p0V, p1V, pAb)
        print testSample, 'classified as:', result
        # return result
###
# 加载数据集
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]  # 1 is abusive, 0 not
    return postingList, classVec
#测试
if __name__=="__main__":
    clf = Naive_Bayes()
    testEntry = [['love', 'my', 'girl', 'friend'],
                 ['stupid', 'garbage'],
                 ['Haha', 'I', 'really', "Love", "You"],
                 ['This', 'is', "my", "dog"],
                 ['maybe','stupid','worthless']]
    for item in testEntry:
        clf.testingNB(item)

使用朴素贝叶斯过滤垃圾邮件

# -*- coding: utf-8 -*-
"""
Created on 下午8:47 22 03 2017
Email_Classify 
@author: plushunter 
"""
import re
import Bayes
from numpy import *
# mysent='This book is the best book on Python or M.L I have ever laid eyes upon.'
# regEx = re.compile('\\W*')
# listOfTokens=regEx.split(mysent)
# tok=[tok.upper() for tok in listOfTokens if len(tok)>0]
# print tok
#
# emailText=open('email/ham/6.txt').read()
# listOfTokens=regEx.split(emailText)
# print listOfTokens
def textParse(bigString):
    import re
    listOfTokens=re.split(r'\w*',bigString)
    return [tok.lower() for tok in listOfTokens if len(tok)>2]
def spamTest():
    clf = Bayes.Naive_Bayes()
    docList=[]
    classList=[]
    fullText=[]
    for i in range(1,26):
        wordList=textParse(open('email/spam/%d.txt'%i).read())
        docList.append(wordList)
        fullText.extend(wordList)
        classList.append(1)
        wordList=textParse(open('email/ham/%i.txt'%i).read())
        docList.append(wordList)
        fullText.extend(wordList)
        classList.append(0)
    vocabList=clf._creatVocabList(docList)
    trainingSet=range(50);testSet=[]
    for i in range(10):
        randIndex=int(random.uniform(0,len(trainingSet)))
        testSet.append(trainingSet[randIndex])
        del(trainingSet[randIndex])
    trainMatix=[];trainClasses=[]
    for docIndex in trainingSet:
        trainMatix.append(clf._bagOfsetOfWordToVec(vocabList,docList[docIndex]))
        trainClasses.append(classList[docIndex])
    p0V,p1V,pSpam=clf._trainNB0(array(trainMatix),array(trainClasses))
    errorCount = 0
    for docIndex in testSet:
        wordVector = clf._bagOfsetOfWordToVec(vocabList,docList[docIndex])
        if clf._classifyNB(array(wordVector), p0V, p1V, pSpam)!=classList[docIndex]:
            errorCount+=1
    print 'the error rate is :',float(errorCount)/len(testSet)

朴素贝叶斯优缺点

优点

具有稳定的分类效率。
对缺失数据不敏感，算法也比较简单。
对小规模数据表现良好，能处理多分类任务，适合增量式训练。尤其是数据量超出内存后，我们可以一批批的去增量训练。

缺点

对输入数据的表达形式比较敏感，需针对不同类型数据采用不同模型。
由于我们是使用数据加先验概率预测后验概率，所以分类决策存在一定的错误率。
假设各特征之间相互独立，但实际生活中往往不成立，因此对特征个数比较多或特征之间相关性比较大的数据来说，分类效果可能不是太好。

最大期望算法

发表于 2018-08-27 | 分类于机器学习专题

字数统计: 3.6k | 阅读时长 ≈ 15

注：这篇文章原文出处放在了下面的连接，有适当的修改和增加内容，写在这里是为了更好地查阅巩固。

戳我原文

EM算法简介

最大期望算法（Expectation Maximization，简称EM算法）是在概率模型中寻找参数最大似然估计或者最大后验估计的算法，其中概率模型依赖于无法观测的隐藏变量。其主要思想就是通过迭代来建立完整数据的对数似然函数的期望界限，然后最大化不完整数据的对数似然函数。最大期望算法是一种迭代优化算法，其计算方法是每次迭代分为期望(E)步和最大(M)步。

EM算法实例

假如现在我们有两枚硬币1和2，随机抛掷后面朝上概率分别为P1，P2。为了估计两硬币概率，我们做如下实验，每次取一枚硬币，连掷5次后得到如下结果。

我们可以很方便的估计出硬币1概率P1=0.4，硬币2概率P2=0.5。

$P1=\frac{3+1+2}{15}=0.4$ $P2=\frac{2+3}{10}=0.5$

下面我们增加问题难度。如果并不知道每次投掷时所使用的硬币标记，那么如何估计P1和P2呢?

此时我们加入隐含变量z，可以把它认为是一个5维的向量(z1,z2,z3,z4,z5)，代表每次投掷时所使用的硬币。比如z1就代表第一轮投掷时所使用的是硬币1还是2，我们必须先估计出z，然后才能进一步估计P1和P2

我们先随机初始化一个P1和P2，用它来估计z，然后基于z按照最大似然概率方法去估计新的P1和P2.例如假设P1=0.2和P2=0.7，然后我们看看第一轮投掷的最可能是哪个硬币。如果是硬币1，得出3正2反的概率为0.2×0.2×0.2×0.8×0.8=0.00512，如果是硬币2，得出3正2反的概率为0.7×0.7×0.7×0.3×0.3=0.03087。然后依次求出其他4轮中的相应概率，接下来便可根据最大似然方法得到每轮中最有可能的硬币。

我们把上面的值作为z的估计值(2,1,1,2,1)，然后按照最大似然概率方法来估计新的P1和P2。得到

$P1=\frac{2+1+2}{15}=0.33$ $P2=\frac{3+3}{10}=0.6$

可以预估，我们继续按照上面思路，用估计出的P1和P2再来估计z，再用z来估计新的P1和P2，反复迭代下去，可以最终得到P1=0.4，P2=0.5。然后无论怎样迭代，P1和P2的值都会保持0.4和0.5不变，于是我们就找到P1和P2的最大似然估计。

上面我们用最大似然方法估计出z值，然后再用z值按照最大似然概率方法估计新的P1和P2。也就是说，我们使用了最有可能的z值，而不是所有的z值。如果考虑所有可能的z值，对每一个z值都估计出新的P1和P2，将每一个z值概率大小作为权重，将所有新的P1和P2分别加权相加，这样估计出的P1和P2是否会更优呢?

但所有的z值共有2^5=32种，我们是否进行32次估计呢？当然不是，我们利用期望来简化运算。

利用上面表格，我们可以算出每轮投掷种使用硬币1或者使用硬币2的概率。比如第一轮使用硬币1的概率

$z_1=\frac{0.00512}{0.00512+0.03087}=0.14$

相应的算出其他4轮的概率。

上表中表示期望值，例如0.86表示此轮中使用硬币2的概率是0.86。前面方法我们按照最大似然概率直接将第一轮估计为硬币2，此时我们更加严谨，只说有0.14的概率是硬币1，有0.86的概率是硬币2。这样我们在估计P1或者P2时，就可以用上全部的数据，而不是部分的数据。此步我们实际上是估计出z的概率分布，称为E步。

按照期望最大似然概率的法则来估计出新的P1和P2。以P1估计为例，第一轮的3正2反相当于有0.14×3=0.42的概率为正，有0.14×2的概率为反。然后依次计算出其他四轮。那么我们可以得到P1概率，可以看到改变了z的估计方法后，新估计出的P1要更加接近0.4，原因是我们使用了所有抛掷的数据，而不是部分的数据。此步中我们根据E步中求出z的概率分布，依据最大似然概率法则去估计P1和P2，称为M步。

$P1=\frac{4.22}{4.22+7.98}=0.35$

上面我们是通过迭代来得到P1和P2，结果更接近真实的P1和P2。

EM算法推导

Jensen不等式

Jensen不等式表述如下：

如果f是凸函数，x是随机变量，那么：$E[f(x)]\geq f(E[x])$ ,特别地，如果f是严格凸函数，那么当且仅当$P(x=E[x])=1$(也就是说x是常量)，$E[f(x)]=f(E[x])$

通过下面这张图，我们可以加深理解：

上图中，函数f是凸函数，X是随机变量，有0.5的概率为a，有0.5的概率是b（就像抛硬币一样）。X的期望值就是a和b的中值了，图中可以看到$E[f(x)]\geq f(E[x])$成立.

极大似然估计

EM算法推导过程中，会使用到极大似然估计参数。

极大似然估计是一种概率论在统计学的应用。已知某个随机样本满足某种概率分布，但是其中具体的参数不清楚，参数估计就是通过若干次试验，观察结果，利用结果推出参数的大概值。极大似然估计建立在这样的思想上：已知某个参数能使这个样本出现的概率最大，我们当然不会再去选择其他小概率的样本，所以干脆就把这个参数作为估计的真实值。

这里再给出求极大似然估计值的一般步骤：

1）写出似然函数；
2）对似然函数取对数，并整理；
3）求导数，令导数为0，得到似然方程；
4）解似然方程，得到的参数即为所求；

EM算法推导

对于m个样本观察数据$x=(x^{(1)},x^{(2)},x^{(3)},…,x^{(m)})$，找出样本的模型参数θ，极大化模型分布的对数似然函数如下所示

$\theta =\arg \max_{\theta} \sum _{i=1}^{m}logP(x^{(i)};\theta)$

如果我们得到的观察数据有未观察到的隐含数据$z=(z^{(1)},z^{(2)},z^{(3)},…,z^{(m)})$，此时我们极大化模型分布的对数似然函数如下

$\theta =\arg \max_{\theta} \sum _{i=1}^{m}logP(x^{(i)};\theta) \\ \\ \\ =\arg \max_{\theta} \sum _{i=1}^{m}log\sum _{z^{(i)}}P(x^{(i)},z^{(i)};\theta)$

上面方程是无法直接求出θ的，因此需要一些特殊技巧，在此我们引入Jensen不等式。

我们再回到上述推导过程，得到如下方程式。

$\sum _{i=1}^{m}log\sum _{z^{(i)}}P(x^{(i)},z^{(i)};\theta) \\ \\ \\ =\sum _{i=1}^{m}log\sum _{z^{(i)}}Q_i(z^{(i)})\frac{P(x^{(i)},z^{(i)};\theta) }{Q_i(z^{(i)})}\ \ \ \ \ \ \ (1) \\ \\ \\ \ge \sum _{i=1}^{m}\sum _{z^{(i)}}Q_i(z^{(i)})log\frac{P(x^{(i)},z^{(i)};\theta) }{Q_i(z^{(i)})} \ \ \ \ \ \ \ (2)$

我们来解释下怎样得到的方程式(1)和方程式(2)，上面(1)式中引入一个未知的新的分布$Q_i(z^{(i)})$，第二式用到Jensen不等式。

首先log函数是凹函数，那么E[f(X)]≤f(E[X])，也就是f(E(X))≥E(f(X))。其中$\sum _{z^{(i)}}Q_i(z^{(i)})\frac{P(x^{(i)},z^{(i)};\theta) }{Q_i(z^{(i)})}$ 是$\frac{P(x^{(i)},z^{(i);\theta})}{Q_i(z^{(i)})}$的数学期望，那么$log\sum _{z^{(i)}}Q_i(z^{(i)})\frac{P(x^{(i)},z^{(i)};\theta) }{Q_i(z^{(i)})}$ 便相当于f(E(X))，同时$\sum _{z^{(i)}}Q_i(z^{(i)})log\frac{P(x^{(i)},z^{(i)};\theta) }{Q_i(z^{(i)})}$相当于E(f(X))，因此我们便得到上述方程式(1)(2)。

如果要满足Jensen不等式的等号，那么需要满足X为常量，即为

$\frac{P(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}=c,\ c为常量$

那么稍加改变能够得到

$c Q_i(z^{(i)})=P(x^{(i)},z^{(i)};\theta)\ ,c为常量$ $\sum_z c\ Q_i(z^{(i)})= \sum _z P(x^{(i)},z^{(i)};\theta)= c$

因此得到下列方程，其中方程(3)利用到条件概率公式。

$Q_i(z^{(i)})=\frac{P(x^{(i)},z^{(i)};\theta)}{c}=\frac{P(x^{(i)},z^{(i)};\theta)}{\sum _z P(x^{(i)},z^{(i)};\theta)}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$ $=\frac{P(x^{(i)},z^{(i)};\theta)}{ P(x^{(i)};\theta)}=P(z^{(i)}|x^{(i)};\theta) \ \ \ \ \ \ \ \ \ \ \ \ \ \ (3)$

如果$Q_i(z^{(i)})=P(z^{(i)}|x^{(i)};\theta)$ 那么第(2)式就是我们隐藏数据的对数似然的下界。如果我们能极大化方程式(2)的下界，则也在尝试极大化我们的方程式(1)。即我们需要最大化下式

$\arg \max _{\theta} \sum _{i=1}^{m}\sum _{z^{(i)}}Q_i(z^{(i)})log\frac{P(x^{(i)},z^{(i)};\theta) }{Q_i(z^{(i)})}$

去掉上式中常数部分，则我们需要极大化的对数似然下界为

$\arg \max _{\theta} \sum _{i=1}^{m}\sum _{z^{(i)}}Q_i(z^{(i)})[log P(x^{(i)},z^{(i)};\theta)-log {Q_i(z^{(i)})}]$ $=\arg \max _{\theta} \sum _{i=1}^{m}\sum _{z^{(i)}}Q_i(z^{(i)})log P(x^{(i)},z^{(i)};\theta)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (4)$

注意到上式中$Q_i(z^{(i)})$ 是一个分布，因此$\sum _{z^{(i)}}Q_i(z^{(i)})log P(x^{(i)},z^{(i)};\theta)$ 可以理解为$logP(x^{(i)},z^{(i)};\theta)$ 基于条件概率分布$Q_i(z^{(i)})$的期望，也就是我们EM算法中E步。极大化方程式(4)也就是我们EM算法中的M步。

到这里，我们推出了在固定参数θ后，使下界拉升的Q(z)的计算公式就是后验概率（条件概率），解决了Q(z)如何选择的问题。此步就是EM算法的E步，目的是建立L(θ)的下界。接下来的M步，目的是在给定Q(z)后，调整θ，从而极大化L(θ)的下界J（在固定Q(z)后，下界还可以调整的更大）。那么一般的EM算法的步骤如下：

第一步：初始化分布参数θ；
第二步：重复E步和M步直到收敛：
E步：根据参数的初始值或上一次迭代的模型参数来计算出的因变量的后验概率（条件概率），其实就是隐变量的期望值，来作为隐变量的当前估计值：
$Q_i\left(z^{\left(i\right)}\right)=p\left(z^{\left(i\right)}|x^{\left(i\right)};\theta\right)$
M步：最大化似然函数从而获得新的参数值：
$\theta :=arg\underset{\theta}{\max}\sum_i{\sum_{z^{\left(i\right)}}{Q_i\left(z^{\left(i\right)}\right)\log\frac{p\left(x^{\left(i\right)},z^{\left(i\right)};\theta\right)}{Q_i\left(z^{\left(i\right)}\right)}}}$

EM算法可以保证收敛到一个稳定点，但是却不能保证收敛到全局的极大值点，因此它是局部最优的算法。当然，如果我们的优化目标L(θ,θj)是凸的，则EM算法可以保证收敛到全局最大值，这点和梯度下降法中迭代算法相同。

Sklearn实现EM算法

高斯混合模型(GMM)使用高斯分布作为参数模型，利用期望最大(EM)算法进行训练，有兴趣了解高斯混合模型的同学可以去这儿。下列代码来自于Sklearn官网GMM模块，利用高斯混合模型确定iris聚类。

import matplotlib as mpl
import matplotlib.pyplot as plt

import numpy as np

from sklearn import datasets
from sklearn.mixture import GaussianMixture
from sklearn.model_selection import StratifiedKFold

colors = ['navy', 'turquoise', 'darkorange']

def make_ellipses(gmm, ax):
    for n, color in enumerate(colors):
        if gmm.covariance_type == 'full':
            covariances = gmm.covariances_[n][:2, :2]
        elif gmm.covariance_type == 'tied':
            covariances = gmm.covariances_[:2, :2]
        elif gmm.covariance_type == 'diag':
            covariances = np.diag(gmm.covariances_[n][:2])
        elif gmm.covariance_type == 'spherical':
            covariances = np.eye(gmm.means_.shape[1]) * gmm.covariances_[n]
        v, w = np.linalg.eigh(covariances)
        u = w[0] / np.linalg.norm(w[0])
        angle = np.arctan2(u[1], u[0])
        angle = 180 * angle / np.pi  # convert to degrees
        v = 2. * np.sqrt(2.) * np.sqrt(v)
        ell = mpl.patches.Ellipse(gmm.means_[n, :2], v[0], v[1],
                                  180 + angle, color=color)
        ell.set_clip_box(ax.bbox)
        ell.set_alpha(0.5)
        ax.add_artist(ell)

iris = datasets.load_iris()

# Break up the dataset into non-overlapping training (75%)
# and testing (25%) sets.
skf = StratifiedKFold(n_splits=4)
# Only take the first fold.
train_index, test_index = next(iter(skf.split(iris.data, iris.target)))


X_train = iris.data[train_index]
y_train = iris.target[train_index]
X_test = iris.data[test_index]
y_test = iris.target[test_index]

n_classes = len(np.unique(y_train))

# Try GMMs using different types of covariances.
estimators = dict((cov_type, GaussianMixture(n_components=n_classes,
                   covariance_type=cov_type, max_iter=20, random_state=0))
                  for cov_type in ['spherical', 'diag', 'tied', 'full'])

n_estimators = len(estimators)

plt.figure(figsize=(3 * n_estimators // 2, 6))
plt.subplots_adjust(bottom=.01, top=0.95, hspace=.15, wspace=.05,
                    left=.01, right=.99)


for index, (name, estimator) in enumerate(estimators.items()):
    # Since we have class labels for the training data, we can
    # initialize the GMM parameters in a supervised manner.
    estimator.means_init = np.array([X_train[y_train == i].mean(axis=0)
                                    for i in range(n_classes)])

    # Train the other parameters using the EM algorithm.
    estimator.fit(X_train)

    h = plt.subplot(2, n_estimators // 2, index + 1)
    make_ellipses(estimator, h)

    for n, color in enumerate(colors):
        data = iris.data[iris.target == n]
        plt.scatter(data[:, 0], data[:, 1], s=0.8, color=color,
                    label=iris.target_names[n])
    # Plot the test data with crosses
    for n, color in enumerate(colors):
        data = X_test[y_test == n]
        plt.scatter(data[:, 0], data[:, 1], marker='x', color=color)

    y_train_pred = estimator.predict(X_train)
    train_accuracy = np.mean(y_train_pred.ravel() == y_train.ravel()) * 100
    plt.text(0.05, 0.9, 'Train accuracy: %.1f' % train_accuracy,
             transform=h.transAxes)

    y_test_pred = estimator.predict(X_test)
    test_accuracy = np.mean(y_test_pred.ravel() == y_test.ravel()) * 100
    plt.text(0.05, 0.8, 'Test accuracy: %.1f' % test_accuracy,
             transform=h.transAxes)

    plt.xticks(())
    plt.yticks(())
    plt.title(name)

plt.legend(scatterpoints=1, loc='lower right', prop=dict(size=12))
plt.show()

EM算法优缺点

优点

聚类。
算法计算结果稳定、准确。
EM算法自收敛，既不需要事先设定类别，也不需要数据间的两两比较合并等操作。

缺点

对初始化数据敏感。
EM算法计算复杂，收敛较慢，不适于大规模数据集和高维数据。
当所要优化的函数不是凸函数时，EM算法容易给出局部最优解，而不是全局最优解。

浴室沉思（二）

发表于 2018-08-20 | 分类于浴室沉思

字数统计: 208 | 阅读时长 ≈ 1

地球那么大，为什么月亮的倒影能够准确地投映在我家院子里的一口井里？

你可以通过逃跑解决肥胖问题。

你看不见你的脖子，你也离不开它。

现金是电子支付的离线缓存。

给水果剥皮，就像是拆来自大地的礼物。

酒店房间里收费的饮品就像是手机应用的内购。

你不能想象一种不存在的颜色。

火山：睡上一万年，生起床气，到处乱丢东西，再睡上一万年。

如果你每天“带薪拉屎”十分钟，一年下来就能相当于休了个一周的“带薪拉屎假”。

10.

如果你的手臂可以拆卸，你也没办法自己把它们互换位置。。

浴室沉思（一）

发表于 2018-08-20 | 分类于浴室沉思

字数统计: 207 | 阅读时长 ≈ 1

总有一天，你会完全忘记你读到过这句话。

音乐是颤抖的风。

狗狗可能并不喜欢叼球回来，只是不想看到东西弄丢了。

根据美国的法律，在合法地饮用第一杯酒之前，你要先绕太阳飞行21圈。

人类殖民火星，就像是通过空气传播的病毒。

因为脑海里的声音不需要换气，脑海里可以一直在尖叫。

同义词没有同义词，但是反义词的反义词是同义词。

闪电是有人躲在云里偷拍我们，但忘了关闪光灯。

氧气有最严重的戒断症状。

10.

火车是一种在三维世界中建立在二维平面上只能一维移动的物体。

浴室沉思101

发表于 2018-08-19 | 分类于浴室沉思

字数统计: 1.1k | 阅读时长 ≈ 3

原帖地址：戳我

欢迎来到r/showerthoughts, 我们很高兴在这里见到你，并且希望你玩得开心。这个帖子是一个关于showerthoughts 是什么/不是什么的简介。

什么是Showerthoughts?

简单来说，Showerthoughts 就是那些让庸常生活充满趣味的想法。这些想法往往是“另一个角度看待日常生活”的结果，而”另一个角度看待日常生活“往往又会让大家看到新的细节。这些想法可能十分可笑，也可能有些苦涩，可能会给你带来思想上的启迪，或者仅仅就是傻乎乎的而已。但他总应该是让人感到“啊，我怎么就没发现呢。“

Showerthoughts, 这个词指的是我们在日常生活中做无脑无聊时候事情时候产生的想法，并不一定要发生在浴室，也可以是在通勤中，在开车，或者是在等客服接你的电话的时候。

这里有一些Showerthoughts的例子：

“你的胃觉得所有的土豆都是土豆泥。”
“当人们提到穿越回过去的时候，总是担心会导致改变今天的世界；但是他们却没有觉得，今天做的一切都会改变未来的世界。“
”如果你够蠢，所有的飞行物都是不明飞行物。”
“C3PO的确是卢克天行者的哥哥。”
“天鹅是一种聒噪，暴力，有进攻性的恐怖生物，同时也是爱情的象征。”

什么不是Showerthoughts？

就像上面说的那样，Showerthoughts 应该让人从一个新的角度看待已经存在的现象或者事物。它不应该是一个个人意见或者观点，或者是一种产品的改进建议，或者是某些能快速通过Google 解决的问题。Showerthoughts 应该是一个命题，而不是一个假设。 R/showerthoughts 里面的内容应该文字流畅，内容新颖。如果某个post 过于泛滥，或者说文字水平太低，抑或者太过粗鲁，管理员将会将其移出。
这些例子里的内容不是Showerthoughts：

“淦！我要买些肥皂了！”

（这是一条浴室观察，不是一条浴室沉思）

“我讨厌开车的时候被别人加塞，但是看着别人被加塞还蛮开心的。“

（这是个人观点）

“有人试过珠穆朗玛峰顶的味道嘛？“

（这已经很接近showerthoughts 了，但是这是一个疑问句，应该改写成一个陈述句。）

“地球每一天都再心的地方”

（文字不流畅，而且还有错别字）

“有没有可能有一群蜜蜂真的有蜂群思维呢？“

（这是个假设，而且还是个文字游戏）

很多时候，我们很难找到一个完全原创，文字流畅，充满洞见的想法，不过这才让showerthoughts 变得好玩：它让你发现在你无意识的时候，你也是在思考的；而且这些想法也让你的无聊日常变得有趣起来。

为什么要把Showerthoughts加入到自己的博客当中？

Showerthoughts是我在Reddit里面最喜欢的一个版块，也是个人很喜欢的一种文化。看似简单的句子里蕴藏了许多奇思妙想，看待事物的角度时常让人耳目一新。看过许多有趣的Showerthoughts之后，偶尔也会自己进行创作，虽然质量不高，但都会发在饭否上面记录。但由于没有统一地方记录这些Showerthoughts，许多之前写的Showerthoughts已经很难找到了。所以，打算在博客里专门创建一个Showerthoughts的Tag，记录收藏有趣的Showerthoughts。

这些Showerthoughts绝大多数会在Reddit上面的r/showerthoughts版块中找到完整的原句，自己只搬运翻译一些感兴趣的Showerthoughts；另外还会摘录一些关注的饭友写的showerthoughts；最后也会自己偶尔写一些较低质量的Showerthoughts。

希望大家在码完代码疲惫之后，可以看到一些有趣的Showerthoughts，消除倦意，继续战斗。

注：微博上也有相关的Showerthoughts搬运博主，偶尔会转发一些网友有趣的Showerthoughts：浴室沉思

浴室沉思101

发表于 2018-08-19 | 分类于浴室沉思

字数统计: 14 | 阅读时长 ≈ 1

思路

1、拆分数据集，验证集，测试集

Sklearn 与 TensorFlow 机器学习实用指南（七）：降维

发表于 2018-08-07 | 分类于 Sklearn 与 TensorFlow 机器学习实用指南

字数统计: 9.9k | 阅读时长 ≈ 36

先来了解一张降维汇总图，降维的算法比较多，这里就只简单说MDS，PCA以及流行学习的Isomap和LLE。

多维缩放MDS：Multiple Dimensional Scaling

低维嵌入：在很多时候，人们观测或收集到的数据样本虽是高维的，但与学习任务密切相关的也许仅是某个低维分布，即高维空间中的一个低维”嵌入” (embedding) . 图 10.2 给出了一个直观的例子. 原始高维空间中的样本点，在这个低维嵌入子空间中更容易进行学习。

MDS的目标是在降维的过程中将数据的dissimilarity(差异性)保持下来，也可以理解降维让高维空间中的距离关系与低维空间中距离关系保持不变。这里的距离用矩阵表示，N个样本的两两距离用矩阵D的每一项$dist_{ij}$ 表示，并且假设在低维空间中的距离是欧式距离。而降维后的数据表示为$Z_i$,那么就有

$dist_{ij}^2 = \left | \mathbf{z_i-z_j} \right |^2 = \left | \mathbf{z_i} \right |^2 + \left | \mathbf{z_j} \right |^2 - 2\mathbf{z_i}\mathbf{z_j}^T$

令$B=ZZ^T$ ,右边的三项统一用内积矩阵B来表示$b_{ij}=z_iz_j^T$ ,所以有

$dist_{ij}^2=b_{ii}+b_{jj}-2b_{ij} \quad(1)$

这时只要求出内积矩阵B即可求出降为后的矩阵Z（思路D-B-Z）。距离矩阵D去中心化之后(减去均值)，內积矩阵B的每一行每一列之和都是0，可以推导得出

其中 tr(.) 表示矩阵的迹(trace),$tr(E)=\sum_{i=1}^m\Vert z_i \Vert^2=\sum_{i=1}^mb_{ii}$,令

联立上(1)-(7),消除$b_{ii},b_{jj}$后，可以得到

$b_{ij} = -\frac{1}{2}(dist_{ij}^2-dist_{i\cdot} - dist_{\cdot j}-dist_{\cdot\cdot}^2)$

i⋅与⋅j是指某列或者某列总和，从而建立了距离矩阵D与内积矩阵B之间的关系.由此即可通过降维前后保持不变的距离矩阵 D 求取内积矩阵 B。对矩阵 B 做特征值分解(eigenvalue decomposition)，$B=VAV^T$，其中$A=diag(λ_1,λ_2,…λ_d)$为特征值构成的对角矩阵,$λ_1\geq λ_2\geq …\geq λ_d$，V 为特征向量矩阵.假定其中有d* 个非零特征值，它们构成对角矩阵$A=diag(λ_1,λ_2,…λ_{d})$,令 V*表示相应的特征向量矩阵，联立之前的$B=ZZ^T$,则最后Z可表达为

$Z=A_*^{1/2}V_*^T$

在现实应用中为了有效降维，往往仅需降维后的距离与原始空间中的距离尽可能接近?而不必严格相等.此时可取 d’<< d 个最大特征值构成对角矩阵。

主成分分析PCA: Principal Component Analysis

向量內积

由$A⋅B=|A||B|cos(a)$, A与B的内积等于A到B的投影长度乘以B的模。再进一步，如果我们假设B的模为1，则A与B的内积值等于A向B所在直线投影的矢量长度！

基

向量(x,y)实际上表示线性组合，$x(1,0)^T+y(0,1)^T$,不难证明所有二维向量都可以表示为这样的线性组合。此处（1，0）和（0，1）叫做二维空间中的一组基。但实际上任何两个线性无关的二维向量都可以成为一组基，所谓线性无关在二维平面内可以直观认为是两个不在一条直线上的向量.例如，(1,1)和(−1,1)也可以成为一组基。一般来说，我们希望基的模是1,实际上，对应任何一个向量我们总可以找到其同方向上模为1的向量，只要让两个分量分别除以模就好了。例如，上面的基可以变为$(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})$和$(-\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})$

基变换的矩阵表示：

下面我们找一种简便的方式来表示基变换。还是拿上面的例子，想一下，将(3,2)变换为新基上的坐标，就是用(3,2)与第一个基做内积运算，作为第一个新的坐标分量，然后用(3,2)与第二个基做内积运算，作为第二个新坐标的分量。实际上，我们可以用矩阵相乘的形式简洁的表示这个变换，其中矩阵的两行分别为两个基，乘以原向量，其结果刚好为新基的坐标：

协方差矩阵及优化目标

面我们讨论了选择不同的基可以对同样一组数据给出不同的表示，而且如果基的数量少于向量本身的维数，则可以达到降维的效果。但是我们还没有回答一个最最关键的问题：如何选择基才是最优的。或者说，如果我们有一组N维向量，现在要将其降到K维（K小于N），那么我们应该如何选择K个基才能最大程度保留原有的信息？

现在问题来了：如果我们必须使用一维来表示上面这些数据，又希望尽量保留原始的信息，你要如何选择？

通过上一节对基变换的讨论我们知道，这个问题实际上是要在二维平面中选择一个方向，将所有数据都投影到这个方向所在直线上，用投影值表示原始记录。这是一个实际的二维降到一维的问题。

那么如何选择这个方向（或者说基）才能尽量保留最多的原始信息呢？一种直观的看法是：希望投影后的投影值尽可能分散。

我们希望投影后投影值尽可能分散，而这种分散程度，可以用数学上的方差来表述。此处，一个字段的方差可以看做是每个元素与字段均值的差的平方和的均值，即：

$Var(a)=\frac{1}{m}\sum^m_{i=1}(a_i-\mu)^2$

通过去中心化（字段所有数值减去字段均值），将每个字段的均值都化为0了，因此方差可以直接用每个元素的平方和除以元素个数表示：

$Var(a)=\frac{1}{m}\sum^m_{i=1}a_i^2$

于是上面的问题被形式化表述为：寻找一个一维基，使得所有数据变换为这个基上的坐标表示后，方差值最大.

协方差

对于上面二维降成一维的问题来说，找到那个使得方差最大的方向就可以了。不过对于更高维，还有一个问题需要解决。考虑三维降到二维问题。与之前相同，首先我们希望找到一个方向使得投影后方差最大，这样就完成了第一个方向的选择，继而我们选择第二个投影方向。

如果我们还是单纯只选择方差最大的方向，很明显，这个方向与第一个方向应该是“几乎重合在一起”，显然这样的维度是没有用的，因此，应该有其他约束条件。从直观上说，让两个字段尽可能表示更多的原始信息，我们是不希望它们之间存在（线性）相关性的，因为相关性意味着两个字段不是完全独立，必然存在重复表示的信息。

数学上可以用两个字段的协方差表示其相关性，由于已经让每个字段均值为0，则：

$Cov(a,b) = \frac{1}{m}\sum_{i=1}^ma_ib_i$

可以看到，在字段均值为0的情况下，两个字段的协方差简洁的表示为其内积除以元素数m。

当协方差为0时，表示两个字段完全独立。为了让协方差为0，我们选择第二个基时只能在与第一个基正交的方向上选择。因此最终选择的两个方向一定是正交的。

至此，我们得到了降维问题的优化目标：将一组N维向量降为K维（K大于0，小于N），其目标是选择K个单位（模为1）正交基，使得原始数据变换到这组基上后，各字段两两间协方差为0，而字段的方差则尽可能大（在正交的约束下，取最大的K个方差）。

协方差矩阵

我们看到，最终要达到的目的与字段内方差及字段间协方差有密切关系。因此我们希望能将两者统一表示，仔细观察发现，两者均可以表示为内积的形式，而内积又与矩阵相乘密切相关。

假设我们只有a和b两个字段，那么我们将它们按行组成矩阵X：

然后我们用X乘以X的转置，并乘上系数1/m：

奇迹出现了！这个对角线上的两个元素分别是两个字段的方差，而其它元素是a和b的协方差。两者被统一到了一个矩阵的。

根据矩阵相乘的运算法则，这个结论很容易被推广到一般情况：设我们有m个n维数据记录，将其按列排成n乘m的矩阵X，设$C=\frac{1}{m}XX^𝖳$ ，则C是一个对称矩阵，其对角线分别个各个字段的方差，而第i行j列和j行i列元素相同，表示i和j两个字段的协方差。

协方差矩阵对角化

根据上述推导，我们发现要达到优化条件，等价于将协方差对角化：即除对角线外的其它元素化为0，并且在对角线上将元素按大小从上到下排列，这样我们就达到了优化目的。这样说可能还不是很明晰，我们进一步看下原矩阵与基变换后矩阵协方差矩阵的关系：

设原始数据矩阵X对应的协方差矩阵为C，而P是一组基按行组成的矩阵，设$Y=PX$，则Y为X对P做基变换后的数据。设Y的协方差矩阵为D，我们推导一下D与C的关系：

$D = \frac{1}{m}YY^T\\ \quad =\frac{1}{m}(PX)(PX)^T\\ =\frac{1}{m}PXX^TP\\=PCP^T$

现在事情很明白了，我们要找的P不是别的，而是能让原始协方差矩阵对角化的P。换句话说，优化目标变成了寻找一个矩阵P，满足$PCP^T$是一个对角矩阵，并且对角元素按从小到大依次排列，那么P的前K行就是要寻找的基，用P的前K行组成的矩阵乘以X就使得X从N维降到了K维并满足上述优化条件。

现在所有焦点都聚焦在了协方差矩阵对角化问题上，由上文知道，协方差矩阵C是一个是对称矩阵，在线性代数上，实对称矩阵有一系列非常好的性质：

1）实对称矩阵不同特征值对应的特征向量必然正交。
2）设特征向量λ重数为r，则必然存在r个线性无关的特征向量对应于λ，因此可以将这r个特征向量单位正交化。

由上面两条可知，一个n行n列的实对称矩阵一定可以找到n个单位正交特征向量，设这n个特征向量为$e_1,e_2,⋯,e_n$，我们将其按列组成矩阵：

$E = (e_1\ e_2 \ ··· \ e_n)$

则对协方差矩阵C有如下结论：

$E^TCE = \varLambda\ =\left[\begin{matrix} \lambda_1& & & \\ & \lambda_2& & \\ & & ···& \\ & & & \lambda_n\\ \end{matrix}\right]$

其中Λ为对角矩阵，其对角元素为各特征向量对应的特征值（可能有重复）。以上结论不再给出严格的数学证明，对证明感兴趣的朋友可以参考线性代数书籍关于“实对称矩阵对角化”的内容。

到这里，我们发现我们已经找到了需要的矩阵P：P是协方差矩阵的特征向量单位化后按行排列出的矩阵，其中每一行都是C的一个特征向量。如果设P按照Λ中特征值的从大到小，将特征向量从上到下排列，则用P的前K行组成的矩阵乘以原始数据矩阵X，就得到了我们需要的降维后的数据矩阵Y

于是，只需对协方差矩阵$XX^T$进行特征值分解，将求得的特征值排序:$λ_1\geq λ_2\geq …\geq λ_d$,再再取前 d’ 个特征值对应的特征向量构成 $W=(w_1,w_2,..,w_d’)$.这就是主成分分析的解.要注意降维后低维空间的维数 d’ 通常是由用户事先指定。

简单回顾一下：给定原始数据矩阵X，其协方差矩阵$C=\frac{1}{M}XX^T$对角元素代表了方差，其余元素代表相关性。假定降维后的数据矩阵为Y，其协方差D会满足对角元素最大，其余元素为0，这样才会代表降维后的要求。设Y=PX，P代表要与矩阵相乘的基，有$D=PCP^T$的关系。这个等式的形式与实对称矩阵对角化一样，这时候只要找出C的特征值对应的特征向量，组合起来即为我们需要求得P。最后Y=PX得到降维后的数据矩阵。

PCA 仅需保留 W 与样本的均值向量即可通过简单的向量减法和矩阵”向量乘法将新样本投影至低维空间中 . 显然，低维空间与原始高维空间必有不同，因为对应于最小的d-d’个特征值的特征向量被舍弃了，这是降维导致的结果.但舍弃这部分信息往往是必要的- 一方面舍弃这部分信息之后能使样本的采样密度增大，这正是降维的重要动机; 另一方面，当数据受到噪声影响时，最小的特征值所对应的特征向量往往与噪声有关?将它们舍弃能在一定程度上起到去噪的效果.

流行学习算法Isomap：等度量映射

等度量映射(Isometric Mapping，简称 Isomap) 的基本出发点，是认为低维流Î~嵌入到高维空间之后，直接在高维空间中计算直线距离具有误导性，因为高维空间中的直线距离在低维嵌入流形上是不可达的.如图所示，低维嵌入流形上两点间的距离是”测地线” (geodesic)距离。

那么，如何计算测地线距离呢?这时我们可利用流形在局部上与欧氏空间同胚这个性质，对每个点基于欧氏距离找出其近邻点，然后就能建立一个近邻连接图，图中近邻点之间存在连接，而非近邻点之间不存在连接，于是，计算两点之间测地线距离的问题就转变为计算近邻连接图上两点之间的最短路径问题.

在近邻连接图上计算两点间的最短路径?可采用著名的Dijkstra算法或Floyd算法，在得到任意两点的距离之后，就可通过MDS 方法来获得样本点在低维空间中的坐标。

流行学习算法LLE：局部线性嵌入

LLE首先假设数据在较小的局部是线性的，也就是说，某一个数据可以由它邻域中的几个样本来线性表示。比如我们有一个样本x1,我们在它的原始高维邻域里用K-近邻思想找到和它最近的三个样本x2,x3,x4. 然后我们假设x1可以由x2,x3,x4线性表示，即

$x_1 = w_{12}x_2 + w_{13}x_3 +w_{14}x_4$

其中，w12，w13，w14为权重系数。在我们通过LLE降维后，我们希望x1在低维空间对应的投影x′1和x2,x3,x4对应的投影x′2,x′3,x′4也尽量保持同样的线性关系，即

$x_1' \approx w_{12}x_2' + w_{13}x_3' +w_{14}x_4'$

也就是说，投影前后线性关系的权重系数w12，w13，w14是尽量不变或者最小改变的。

从上面可以看出，线性关系只在样本的附近起作用，离样本远的样本对局部的线性关系没有影响，因此降维的复杂度降低了很多。

对于LLE算法，我们首先要确定邻域大小的选择，即我们需要多少个邻域样本来线性表示某个样本。假设这个值为k。我们可以通过和KNN一样的思想通过距离度量比如欧式距离来选择某样本的k个最近邻。

在寻找到某个样本的xi的k个最近邻之后我们就需要找到找到xi和这k个最近邻之间的线性关系，也就是要找到线性关系的权重系数。找线性关系，这显然是一个回归问题。假设我们有m个n维样本{x1,x2,…,xm},我们可以用均方差作为回归问题的损失函数：即：

$J(w) = \sum\limits_{i=1}^{m}||x_i-\sum\limits_{j=1}^{k}w_{ij}x_j||_2^2$

一般我们也会对权重系数$w_{ij}$做归一化的限制，即权重系数需要满足

$\sum\limits_{j=1}^{k}w_{ij} = 1$

对于不在样本$x_i$邻域内的样本$x_j$，我们令对应的$w_{ij}=0$.也就是我们需要通过上面两个式子求出我们的权重系数。一般我们可以通过矩阵和拉格朗日子乘法来求解这个最优化问题。(这个推导就不写了）

最后得到

$W_i = \frac{Z_i^{-1}1_k}{1_k^TZ_i^{-1}1_k}$

其中$W_i=(w_{i1},w_{i2},…w_{ik})^T$ ,矩阵$Z_i=(x_i−x_j)^T(x_i−x_j)$,其中$1_k$ 为k维全1向量。

在我们得到了高维的权重系数，那么我们希望这些权重系数对应的线性关系在降维后的低维一样得到保持。假设我们的n维样本集{x1,x2,…,xm}在低维的d维度对应投影为{y1,y2,…,ym}, 则我们希望保持线性关系，也就是希望对应的均方差损失函数最小，即最小化损失函数J(Y)如下：

$J(y) = \sum\limits_{i=1}^{m}||y_i-\sum\limits_{j=1}^{k}w_{ij}y_j||_2^2$

这个优化目标与之前的同形，唯一的区别是之前需要确定权重系数$w_i$，而现在是知道权重系数，需要确定的是$x_i$对应的低维空间坐标$y_i$。

令$M=(I-W)^T(I-W)$ ,则优化函数转变为最小化下式：$J(Y) = tr(Y^TMY)$,tr为迹函数。约束函数矩阵化为：$Y^TY=mI$

上式可通过特征值分解求解:M 最小的 d’ 个特征值对应的特征向量组成的矩阵即为 $Z^T$.
算法流程：

从图中可以看出，LLE算法主要分为三步，第一步是求K近邻的过程，这个过程使用了和KNN算法一样的求最近邻的方法。第二步，就是对每个样本求它在邻域里的K个近邻的线性关系，得到线性关系权重系数W，第三步就是利用权重系数来在低维里重构样本数据。

应用

我们将会展示两种主要的降维方法：投影（projection）和流形学习（Manifold Learning），同时我们还会介绍三种流行的降维技术：主成分分析（PCA），核主成分分析（Kernel PCA）和局部线性嵌入（LLE）。

主成分分析（PCA）

主成分分析（Principal Component Analysis）是目前为止最流行的降维算法。首先它找到接近数据集分布的超平面，然后将所有的数据都投影到这个超平面上。

保留（最大）方差

在将训练集投影到较低维超平面之前，您首先需要选择正确的超平面。例如图左侧是一个简单的二维数据集，以及三个不同的轴（即一维超平面）。图右边是将数据集投影到每个轴上的结果。正如你所看到的，投影到实线上保留了最大方差，而在点线上的投影只保留了非常小的方差，投影到虚线上保留的方差则处于上述两者之间。

选择保持最大方差的轴看起来是合理的，因为它很可能比其他投影损失更少的信息。证明这种选择的另一种方法是，选择这个轴使得将原始数据集投影到该轴上的均方距离最小。这是就 PCA 背后的思想，相当简单。

主成分（Principle Componets）

PCA 寻找训练集中可获得最大方差的轴。在上图中，它是一条实线。它还发现了一个与第一个轴正交的第二个轴，选择它可以获得最大的残差。在这个 2D 例子中，没有选择：就只有这条点线。但如果在一个更高维的数据集中，PCA 也可以找到与前两个轴正交的第三个轴，以及与数据集中维数相同的第四个轴，第五个轴等。定义第i个轴的单位矢量被称为第i个主成分（PC）。在图中，第一个 PC 是c1，第二个 PC 是c2。在投影图中，前两个 PC 用平面中的正交箭头表示，第三个 PC 与上述 PC 形成的平面正交（指向上或下）

概述：主成分的方向不稳定：如果您稍微打乱一下训练集并再次运行 PCA，则某些新 PC 可能会指向与原始 PC 方向相反。但是，它们通常仍位于同一轴线上。在某些情况下，一对 PC 甚至可能会旋转或交换，但它们定义的平面通常保持不变。

那么如何找到训练集的主成分呢？幸运的是，有一种称为奇异值分解（SVD）的标准矩阵分解技术，可以将训练集矩阵X分解为三个矩阵$U·Σ·V^T$的点积，其中$V^T$$包含我们想要的所有主成分，如下所示。

下面的 Python 代码使用了 Numpy 提供的svd()函数获得训练集的所有主成分，然后提取前两个 PC:

X_centered=X-X.mean(axis=0)    # 中心化
U,s,V=np.linalg.svd(X_centered)
c1=V.T[:,0]
c2=V.T[:,1]

警告：PCA 假定数据集以原点为中心。正如我们将看到的，Scikit-Learn 的PCA类负责为您的数据集中心化处理。但是，如果您自己实现 PCA（如前面的示例所示），或者如果您使用其他库，不要忘记首先要先对数据做中心化处理。

投影到d维空间：
一旦确定了所有的主成分，你就可以通过将数据集投影到由前d个主成分构成的超平面上，从而将数据集的维数降至d维。选择这个超平面可以确保投影将保留尽可能多的方差。

为了将训练集投影到超平面上，可以简单地通过计算训练集矩阵X和Wd的点积，Wd定义为包含前d个主成分的矩阵（即由V^T的前d列组成的矩阵）

将训练集投影到d维空间的公式：

$X_{d-proj} = X \cdot W_d$

下面的 Python 代码将训练集投影到由前两个主成分定义的超平面上：

1 2	W2=V.T[:,:2] # 降为2维 X2D=X_centered.dot(W2)

使用 Scikit-Learn
Scikit-Learn 的 PCA 类使用 SVD 分解来实现，就像我们之前做的那样。以下代码应用 PCA 将数据集的维度降至两维（请注意，它会自动处理数据的中心化）：

from sklearn.decomposition import PCA

pca=PCA(n_components=2)
X2D=pca.fit_transform(X)

将 PCA 转化器应用于数据集后，可以使用components_访问每一个主成分（注意，它返回以 PC 作为水平向量的矩阵，因此，如果我们想要获得第一个主成分则可以写成pca.components_.T[:,0]）。

方差解释率（Explained Variance Ratio）

另一个非常有用的信息是每个主成分的方差解释率，可通过explained_variance_ratio_变量获得。它表示位于每个主成分轴上的数据集方差的比例。例如，让我们看下图中表示的三维数据集前两个分量的方差解释率：

1 2	>>> print(pca.explained_variance_ratio_) array([0.84248607, 0.14631839])

这表明，84.2% 的数据集方差位于第一轴，14.6% 的方差位于第二轴。第三轴的这一比例不到1.2％，因此可以认为它可能没有包含什么信息

选择正确的维度
通常我们倾向于选择加起来到方差解释率能够达到足够占比（例如 95%）的维度的数量，而不是任意选择要降低到的维度数量。当然，除非您正在为数据可视化而降低维度 — 在这种情况下，您通常希望将维度降低到 2 或 3。

下面的代码在不降维的情况下进行 PCA，然后计算出保留训练集方差 95% 所需的最小维数：

pca=PCA()
pac.fit(X)
cumsum=np.cumsum(pca.explained_variance_ratio_)
d=np.argmax(cumsum>=0.95)+1

你可以设置n_components = d并再次运行 PCA。但是，有一个更好的选择：不指定你想要保留的主成分个数，而是将n_components设置为 0.0 到 1.0 之间的浮点数，表明您希望保留的方差比率：

1 2	pca=PCA(n_components=0.95) X_reduced=pca.fit_transform(X)

另一种选择是画出方差解释率关于维数的函数（简单地绘制cumsum）。曲线中通常会有一个肘部，方差解释率停止快速增长。您可以将其视为数据集的真正的维度。在这种情况下，您可以看到将维度降低到大约100个维度不会失去太多的可解释方差。

PCA 压缩

显然，在降维之后，训练集占用的空间要少得多。例如，尝试将 PCA 应用于 MNIST 数据集，同时保留 95% 的方差。你应该发现每个实例只有 150 多个特征，而不是原来的 784 个特征。因此，尽管大部分方差都保留下来，但数据集现在还不到其原始大小的 20%！这是一个合理的压缩比率，您可以看到这可以如何极大地加快分类算法（如 SVM 分类器）的速度。

通过应用 PCA 投影的逆变换，也可以将缩小的数据集解压缩回 784 维。当然这并不会返回给你最原始的数据，因为投影丢失了一些信息（在5％的方差内），但它可能非常接近原始数据。原始数据和重构数据之间的均方距离（压缩然后解压缩）被称为重构误差（reconstruction error）。例如，下面的代码将 MNIST 数据集压缩到 154 维，然后使用inverse_transform()方法将其解压缩回 784 维。图 8-9 显示了原始训练集（左侧）的几位数字在压缩并解压缩后（右侧）的对应数字。您可以看到有轻微的图像质量降低，但数字仍然大部分完好无损。

1
2
3

pca=PCA(n_components=154)
X_mnist_reduced=pca.fit_transform(X_mnist)
X_mnist_recovered=pca.inverse_transform(X_mnist_reduced)

PCA逆变换公式，回退到原来的数据维度

$X_{recovered} = X_{d-proj} \cdot W_d^T$

增量 PCA（Incremental PCA）

先前 PCA 实现的一个问题是它需要在内存中处理整个训练集以便 SVD 算法运行。幸运的是，我们已经开发了增量 PCA（IPCA）算法：您可以将训练集分批，并一次只对一个批量使用 IPCA 算法。这对大型训练集非常有用，并且可以在线应用 PCA（即在新实例到达时即时运行）。

下面的代码将 MNIST 数据集分成 100 个小批量（使用 NumPy 的array_split()函数），并将它们提供给 Scikit-Learn 的IncrementalPCA类，以将 MNIST 数据集的维度降低到 154 维（就像以前一样）。请注意，您必须对每个最小批次调用partial_fit()方法，而不是对整个训练集使用fit()方法

from sklearn.decomposition import IncrementalPCA

n_batches=100
inc_pca=IncrementalPCA(n_components=154)
for X_batch in np.array_spplit(X_mnist,n_batches): #分100批次数据
    inc_pca.partial_fit(X_batch)    # 必须批次调用partial_fit()方法
X_mnist_reduced=inc_pca.transform(X_mnist)

或者，您可以使用 NumPy 的memmap类，它允许您操作存储在磁盘上二进制文件中的大型数组，就好像它完全在内存中；该类仅在需要时加载内存中所需的数据。由于增量 PCA 类在任何时间内仅使用数组的一小部分，因此内存使用量仍受到控制。这可以调用通常的fit()方法，如下面的代码所示：

X_mm=np.memmap(filename,dtype='float32',mode='readonly',shape=(m,n))
batch_size=m//n_batches
inc_pca=IncrementalPCA(n_components=154,batch_size=batch_size)
inc_pca.fit(X_mm)

随机 PCA（Randomized PCA）

Scikit-Learn 提供了另一种执行 PCA 的选择，称为随机 PCA。这是一种随机算法，可以快速找到前d个主成分的近似值。它的计算复杂度是O(m × d^2) + O(d^3)，而不是O(m × n^2) + O(n^3)，所以当d远小于n时，它比之前的算法快得多。

1 2	rnd_pca=PCA(n_components=154,svd_solver='randomized') X_reduced=rnd_pca.fit_transform(X_mnist)

核 PCA（Kernel PCA）

例如，下面的代码使用 Scikit-Learn 的KernelPCA类来执行带有 RBF 核的 kPCA

from sklearn.decomposition import KernelPCA

rbf_pca=KernelPCA(n_components=2,kernel='rbf',gamma=0.04)
X_reduced=rbf_pca.fit_transform(X)

选择一种核并调整超参数
由于 kPCA 是无监督学习算法，因此没有明显的性能指标可以帮助您选择最佳的核方法和超参数值。但是，降维通常是监督学习任务（例如分类）的准备步骤，因此您可以简单地使用网格搜索来选择可以让该任务达到最佳表现的核方法和超参数。例如，下面的代码创建了一个两步的流水线，首先使用 kPCA 将维度降至两维，然后应用 Logistic 回归进行分类。然后它使用Grid SearchCV为 kPCA 找到最佳的核和gamma值，以便在最后获得最佳的分类准确性：

from sklearn.model_selection import GridSearchCV from sklearn.linear_model import LogisticRegression from sklearn.pipeline import Pipeline
clf = Pipeline([
        ("kpca", KernelPCA(n_components=2)),
        ("log_reg", LogisticRegression())
])
param_grid = [{
        "kpca__gamma": np.linspace(0.03, 0.05, 10),
        "kpca__kernel": ["rbf", "sigmoid"]
    }]
grid_search = GridSearchCV(clf, param_grid, cv=3)
grid_search.fit(X, y)

你可以通过调用best_params_变量来查看使模型效果最好的核和超参数：

1 2	>>> print(grid_search.best_params_) {'kpca__gamma': 0.043333333333333335, 'kpca__kernel': 'rbf'}

另一种完全为非监督的方法，是选择产生最低重建误差的核和超参数。但是，重建并不像线性 PCA 那样容易。这里是原因：图 8-11 显示了原始瑞士卷 3D 数据集（左上角），并且使用 RBF 核应用 kPCA 后生成的二维数据集（右上角）。由于核技巧，这在数学上等同于使用特征映射φ将训练集映射到无限维特征空间（右下），然后使用线性 PCA 将变换的训练集投影到 2D。请注意，如果我们可以在缩减空间中对给定实例实现反向线性 PCA 步骤，则重构点将位于特征空间中，而不是位于原始空间中（例如，如图中由x表示的那样）。由于特征空间是无限维的，我们不能找出重建点，因此我们无法计算真实的重建误差。幸运的是，可以在原始空间中找到一个贴近重建点的点。这被称为重建前图像（reconstruction pre-image）。一旦你有这个前图像，你就可以测量其与原始实例的平方距离。然后，您可以选择最小化重建前图像错误的核和超参数。

您可能想知道如何进行这种重建。一种解决方案是训练一个监督回归模型，将预计实例作为训练集，并将原始实例作为训练目标。如果您设置了fit_inverse_transform = True，Scikit-Learn 将自动执行此操作，代码如下所示：

1
2
3

rbf_pca = KernelPCA(n_components = 2, kernel="rbf", gamma=0.0433,fit_inverse_transform=True)
X_reduced = rbf_pca.fit_transform(X)
X_preimage = rbf_pca.inverse_transform(X_reduced)

概述：默认条件下，fit_inverse_transform = False并且KernelPCA没有inverse_tranfrom()方法。这种方法仅仅当fit_inverse_transform = True的情况下才会创建。

你可以计算重建前图像误差：

1 2	>>> from sklearn.metrics import mean_squared_error >>> mean_squared_error(X, X_preimage) 32.786308795766132

现在你可以使用交叉验证的方格搜索来寻找可以最小化重建前图像误差的核方法和超参数。

局部线性嵌入LLE

局部线性嵌入（Locally Linear Embedding）是另一种非常有效的非线性降维（NLDR）方法。这是一种流形学习技术，不依赖于像以前算法那样的投影。简而言之，LLE 首先测量每个训练实例与其最近邻（c.n.）之间的线性关系，然后寻找能最好地保留这些局部关系的训练集的低维表示（稍后会详细介绍）。这使得它特别擅长展开扭曲的流形，尤其是在没有太多噪音的情况下。

例如，以下代码使用 Scikit-Learn 的LocallyLinearEmbedding类来展开瑞士卷。得到的二维数据集如图所示。正如您所看到的，瑞士卷被完全展开，实例之间的距离保存得很好。但是，距离不能在较大范围内保留的很好：展开的瑞士卷的左侧被挤压，而右侧的部分被拉长。尽管如此，LLE 在对流形建模方面做得非常好。

from sklearn.manifold import LocallyLinearEmbedding

lle=LocallyLinearEmbedding(n_components=2,n_neighbors=10)
X_reduced=lle.fit_transform(X)

其他降维方法

多维缩放（MDS）在尝试保持实例之间距离的同时降低了维度

Isomap 通过将每个实例连接到最近的邻居来创建图形，然后在尝试保持实例之间的测地距离时降低维度。

t-分布随机邻域嵌入（t-Distributed Stochastic Neighbor Embedding，t-SNE）可以用于降低维度，同时试图保持相似的实例临近并将不相似的实例分开。它主要用于可视化，尤其是用于可视化高维空间中的实例（例如，可以将MNIST图像降维到 2D 可视化）。

线性判别分析（Linear Discriminant Analysis，LDA）实际上是一种分类算法，但在训练过程中，它会学习类之间最有区别的轴，然后使用这些轴来定义用于投影数据的超平面。LDA 的好处是投影会尽可能地保持各个类之间距离，所以在运行另一种分类算法（如 SVM 分类器）之前，LDA 是很好的降维技术。

练习题

减少数据集维度的主要动机是什么？主要缺点是什么？
什么是维度爆炸？
一旦对某数据集降维，我们可能恢复它吗？如果可以，怎样做才能恢复？如果不可以，为什么？
PCA 可以用于降低一个高度非线性对数据集吗？
假设你对一个 1000 维的数据集应用 PCA，同时设置方差解释率为 95%，你的最终数据集将会有多少维？
在什么情况下你会使用普通的 PCA，增量 PCA，随机 PCA 和核 PCA？
你该如何评价你的降维算法在你数据集上的表现？
将两个不同的降维算法串联使用有意义吗？

1、动机：为了加速后续训练算法（在某些情况下，它甚至可以消除噪声和冗余特征，使训练算法表现更好）；通过可视化数据深入了解最重要的特征；节省空间（压缩）。缺点：某些信息丢失，可能会降低后续训练算法的性能；计算密集；为机器学习管道增加了一些复杂性；转换后的功能通常难以解释。

2、在机器学习中，一个常见的表现是随机采样的高维向量通常非常稀疏，增加了过拟合的风险，并且在没有足够的训练数据的情况下很难识别数据中的模式。

3、一旦使用我们讨论过的算法减少了数据集的维数，几乎总是不可能完全逆转操作，因为在降维期间某些信息会丢失。此外，虽然一些算法（例如PCA）具有可以重建与原始数据相对类似的数据集的简单反向变换过程，但是其他算法（例如T-SNE）则不然。

4、PCA可用于显着降低大多数数据集的维度，即使它们是高度非线性的，因为它至少可以消除无用的维度。但是，如果没有无用的维度 - 例如，瑞士卷 - 那么使用PCA降低维数将失去太多信息。你想要展开瑞士卷，而不是挤压它。

5、这是一个棘手的问题：它取决于数据集。让我们看看两个极端的例子。首先，假设数据集由几乎完全对齐的点组成。在这种情况下，PCA可以将数据集减少到一维，同时仍然保留95％的方差。现在想象一下，数据集由完全随机的点组成，分散在1000个维度周围。在这种情况下，需要所有1,000个维度来保持95％的方差。所以答案是，它取决于数据集，它可以是1到1,000之间的任何数字。将解释的方差绘制为维数的函数是一种粗略了解数据集内在维度的方法。

6、常规PCA是默认值，但仅当数据集有足够内存时才有效。当您需要在每次新实例到达时动态应用PCA，增量PCA对于在线任务很有用。当您想要显着降低维度并且有足够内存时，随机PCA非常有用; 在这种情况下，它比普通PCA快得多。最后，Kernel PCA对非线性数据集非常有用。

7、直观地说，如果从数据集中消除了大量维度而不会丢失太多信息，则降维算法表现良好。衡量这种情况的一种方法是应用反向变换并测量重建误差。但是，并非所有降维算法都提供逆向变换。

8、链接两个不同的降维算法绝对有意义。一个常见的例子是使用PCA快速摆脱大量无用的维度，然后应用另一个慢得多的降维算法，如LLE。这种两步法可能会产生与仅使用LLE相同的性能，但只需要很短的时间。

Sklearn 与 TensorFlow 机器学习实用指南（六）：集成学习

发表于 2018-08-04 | 分类于 Sklearn 与 TensorFlow 机器学习实用指南

字数统计: 12.2k | 阅读时长 ≈ 45

集成学习(ensemble learning)通过构建并结合多个学习器来完成学习任务，比单一学习器获得显著优越的泛化性能。想要获得好的集成，个体学习器应”好而不同“，要保证准确性和多样性。要产生好而不同的个体学习器，恰是集成学习研究的核心

目前集成学习可分为两大类，即个体学习器之间有依赖关系，必须串行生成的序列化方法；以及个体学习器不存在强依赖关系，可同时生成的并行化方法。前者的代表是Boosting，最著名的是代表有Adaboost, GBDT和XGBOOST;后者的代表是Bagging和随机森林。对于学习器的结合策略有三大类：投票法（分类），平均法（连续数值），学习法（Stacking）

Boosting

Adaboost

算法思想

Adaboost提升方法是改变训练数据的概率分布（训练数据的权值分布），针对不同的训练数据分布调用弱学习算法学习一系列弱分类器。AdaBoost的做法是，提高那些被前一轮弱分类器错误分类样本的权值，而降低那些被正确分类样本的权值。这样一来，那些没有得到正确分类的数据，由于其权值的加大而受到后一轮的弱分类器的更大关注，于是，分类问题就被一系列的弱分类器“分而治之”。另外，对于弱分类器的组合，AdaBoost采取加权多数表决的方法。具体地，加大分类误差率小的弱分类器的权值，使其在表决中起较大的作用，减小分类误差率较大的弱分类器的权值，使其在表决中起较小的作用。(两个权重，一个是样本权重，另外一个是分类器的权重)

AdaboostBoost的算法的框架如下图所示

算法流程

具体算法流程如下图所示：

step2. 预期迭代T轮
step4. 计算分类器$h_t$的误差率$\epsilon_t$

$\epsilon_t=P\left(h_t\left(x\right)\ne y_i\right)=\frac{\sum_{i=1}^m{w_{ki}I\left(h_t\left(x_i\right)\ne y_i\right)}}{\sum_{i=1}^m{w_{ki}}}=\sum_{i=1}^m{w_{ki}I\left(G_m\left(x_i\right)\ne y_i\right)}$

这里$w_{ki}$表示第k轮（第k个分类器）中第i个实例的权重，$\sum_{i=1}^m{w_{ki}=1}$，I表示指示函数，代表满足条件的样本。这表明，误差率是被$h_t$分类错误的样本的权重之和。（这些样本的权重会在后面归一化）

step5. 分类器比随机分类效果还差则停止
step6. 根据误差率计算分类器的权重，表示最终分类器的重要程度。由表达式可知，当误差率小于等于1/2时，$\alpha_k$大于等于0。并且$\alpha_k$随着误差率的减小而增大，意味着误差越小的分类器最后的重要程度越大。
step7. 更新样本权重，$Z_k$为归一化因子，把最后的全部样本权重求和即可。

$w_{k+1,i}=\frac{w_{ki}}{Z_k}exp(-\alpha_k y_i G_k(x_i))$

Adaboost算法优缺点

Adaboost优点

不容易发生过拟合。
Adaboost是一种有很高精度的分类器。
当使用简单分类器时，计算出的结果是可理解的。
可以使用各种方法构建子分类器，Adaboost算法提供的是框架。

Adaboost缺点

训练时间过长。
执行效果依赖于弱分类器的选择。
对样本敏感，异常样本在迭代中可能会获得较高的权重，影响最终的强学习器的预测准确性。

GBTD

GBDT（Gradient Boosting Decision Tree）是一种迭代的决策树算法，又叫 MART（Multiple Additive Regression Tree)，它通过构造一组弱的学习器（树），并把多颗决策树的结果累加起来作为最终的预测输出。该算法将决策树与集成思想进行了有效的结合。并且GBDT每一棵树都是回归树CART.。由于GBDT的核心在与累加所有树的结果作为最终结果，而分类树得到的离散分类结果对于预测分类并不是这么的容易叠加。这是区别于分类树的一个显著特征（毕竟男+女=是男是女?，这样的运算是毫无道理的），GBDT在运行时就使用到了回归树的这个性质，它将累加所有树的结果作为最终结果。所以GBDT中的树都是回归树，而不是分类树，它用来做回归预测，当然回归树经过调整之后也能用来做分类。

这里要先介绍GBDT简单版本的提升树Boosting Decision Tree，后面再介绍GBDT。

提升树Boosting Decision Tree

提升树(Boosting Decision Tree)由于输出样本是连续值，因此我们通过迭代多棵回归树来共同决策（之前CART只是拟合一颗完整的回归树）。回归树的构造在上一节已经介绍过了，不再赘述。

我们利用平方误差来表示损失函数，其中每一棵回归树学习的是之前所有树的结论和残差，拟合得到一个当前的残差回归树。其中残差=真实值-预测值，提升树即是整个迭代过程生成的回归树的累加。提升树模型可以表示为决策树的加法模型：

$f_M\left(x\right)=\sum_{m=1}^M{T\left(x;\varTheta_m\right)}$

其中$T\left(x;\varTheta_m\right)$表示决策树；$\varTheta_m$为决策树的参数；M为树的个数。

提升树的过程如下，节点下所有点的均值作为该节点的预测值，例如左图的15与25。（这里是将特征分开处理并缩小了树的规模，若用CART可能会出现深度为3的回归树）

回归问题的提升树算法叙述如下：

对比初始的CART回归树与GBDT所生成的回归树，可以发现，最终的结果可能是相同的，那我们为什么还要使用GBDT呢？

答案就是对模型过拟合的考虑。过拟合是指为了让训练集精度更高，学到了很多“仅在训练集上成立的规律”，导致换一个数据集后，当前规律的预测精度就不足以使人满意了。毕竟，在训练精度和实际精度（或测试精度）之间，后者才是我们想要真正得到的。
在上面这个例子中，初始的回归树为达到100%精度使用了3个特征（上网时长、时段、网购金额），但观察发现，分枝“上网时长>1.1h”很显然过拟合了，不排除恰好A上网1.5h, B上网1小时，所以用上网时间是不是>1.1小时来判断所有人的年龄很显然是有悖常识的。
而在GBDT中，两棵回归树仅使用了两个特征（购物金额与对百度知道的使用方式）就实现了100%的预测精度，其分枝依据更合乎逻辑（当然这里是相比较于上网时长特征而言），算法在运行中也体现了“如无必要，勿增实体”的奥卡姆剃刀原理

梯度提升决策树Gradient Boosting Decision Tree

GBDT回归算法

提升树利用加法模型与向前分布算法实现学习的优化过程，即是通过迭代得到一系列的弱分类器，进而通过不同的组合策略得到相应的强学习器。在GBDT的迭代中，假设前一轮得到的强学习器为$f_{t−1}(x)$ ，对应的损失函数则为$L(y,f_{t−1}(x))$ 。因此新一轮迭代的目的就是找到一个弱学习器$h_t(x)$ ，使得损失函$L(y,f_{t−1}(x)+h_t(x))$ 达到最小。因此问题的关键就在于对损失函数的度量，这也正是难点所在。当损失函数是平方损失和指数损失时，每一步优化是很简单的。但对一般损失函数而言，往往每一步优化没那么容易，如绝对值损失函数和Huber损失函数。常见的损失函数及其梯度如下表所示：

那我们怎么样才能找到一种通用的拟合方法呢？针对这一问题，Freidman提出了梯度提升算法：利用最速下降的近似方法，即利用损失函数的负梯度在当前模型的值，进而拟合一个CART回归树。其中第t轮的第i个样本的损失函数的扶梯度表示为，右下角的等式是求偏导后带入计算的

$r_{ti}=-\left[\frac{\partial L\left(y,f\left(x_i\right)\right)}{\partial f\left(x_i\right)}\right]_{f\left(x\right)=f_{t-1}\;\,\left(x\right)}$

负梯度作为回归问题中提升树算法的残差的近似值（与其说负梯度作为残差的近似值，不如说残差是负梯度的一种特例，拟合一个回归树），这就是梯度提升决策树。假设样本数据为m，最大的迭代次数为T。其算法过程如下：

step1. 初始化弱分类器，计算出使损失函数极小化的一个常数值c，此时树仅有一个根结点。c的均值可取样本y的均值。这个初始化得到的c将用于第一次计算负梯度$f(x_i)=f(t-1)=f(0)$的代入计算得到残差近似值$r_{ti}$。(第t代第i个样本)
step2(a) .计算每个样本的$r_{ti}$
step2(b). 利用$(x_i,r_{ti})i=1,2,3,…,m$，我们可以拟合一颗CART回归树，得到第t棵回归树（注意回归树节点内的均值求法不再是我们想要的了），其对应的叶节点区域为$R_{tj},j=1,2,3,…,J$，其中J为叶子节点的个数。
step2(c). 接下来，针对每一个叶子节点中的样本，要拟合叶子结点最好的输出值$c_{tj}$(不再是简单的求节点均值)，使得求出的损失函数最小。回顾之前写的新一轮迭代的目的，这时候的输出值$c_{tj}$组成就是我们想要的第t棵弱学习器$h_t(x)$。其实就是在上一棵强学习器树稍加改变决策树中叶节点值，希望拟合的误差越来越小。

$c_{tj}=\underset{c}{\underbrace{\arg\min}} \sum_{x_i \in R_{tj}}L(y_i,f_{t-1}(x_i)+c)$

这样我们便得到本轮的弱学习器决策树拟合函数

$h_t(x)=\sum _{j=1} ^{J} c_{tj},I(x \in R_{tj})$

step2(d). 更新强学习器，上一个强学习器+弱学习器

$f_t(x)=f_{t-1}(x)+\sum_{j=1}^{J}c_{tj},I(x\in R_{tj})$

step3. 得到输出的最后一轮的最终强学习器模型（最大迭代T轮）

$f(x)=f_T(x)=f_0(x)+\sum_{t=1}^{T}\sum_{j=1}^{J}c_{tj},I(x\in R_{tj})$

GBDT分类算法

GBDT分类算法在思想上和回归算法没有区别，但是由于样本输出不是连续的值，而是离散的类别，导致我们无法直接从输出类别去拟合类别输出的误差。为解决此问题，我们尝试用类似于逻辑回归的对数似然损失函数的方法,也就是说我们用的是类别的预测概率值和真实概率值来拟合损失函数。对于对数似然损失函数，我们有二元分类和多元分类的区别。

二元GBDT分类算法

对于二元GBDT，如果用类似于逻辑回归的对数似然损失函数，则损失函数表示为

$L(y,f(x))=log(1+exp(-yf(x)))$

其中y∈{−1,1}。则此时的负梯度误差为

$r_{ti}=-\left[\frac{\partial L(y_i,f(x_i))}{\partial f(x_i)} \right]_{f(x)=f_{t-1}(x)}=\frac{y_i}{1+exp(y_if(x_i))}$

对于生成的决策树，我们各个叶子节点的最佳残差拟合值为

$c_{tj}=\underset{c}{\underbrace{\arg\min}} \sum _{x_i\in R_{tj}}log(1+exp(-y_i(f_{t-1}(x_i)+c)))$

由于上式比较难优化，我们一般使用近似值代替

$c_{tj}=\frac{\sum _{x_i\in R_{tj}}r_{ti}}{\sum _{x_i \in R_{tj}}|r_{ti}|(1-|r_{ti}|)}$

除了负梯度计算和叶子节点的最佳残差拟合的线性搜索外，二元GBDT分类和GBDT回归算法过程相同。

多元GBDT分类算法

多元GBDT要比二元GBDT复杂一些，对应的是多元逻辑回归和二元逻辑回归的复杂度差别。假如类别数为K，则我们的对数似然函数为

$L(y,f(x))=-\sum_{k=1}^{K}y_k log(p_k(x))$

其中如果样本输出类别为k，则$y_k=1$。第k类的概率$p_k(x)$的表达式为

$p_k(x)=\frac {exp(f_k(x))}{\sum _{l=1}^{K}exp(f_l(x))}$

集合上两式，我们可以计算出第t轮的第i个样本对应类别l的负梯度误差为

$r_{til}=-\left[\frac{\partial L(y_i,f(x_i))}{\partial f(x_i)} \right]_{f_k(x)=f_{t-1,l}(x)}=y_{il}-p_{t-1,l}(x_i)$

其实这里的误差就是样本i对应类别l的真实概率和t-1轮预测概率的差值。对于生成的决策树，我们各个叶子节点的最佳残差拟合值为

$c_{tjl}=\underset{cjl}{\underbrace{\arg\min}} \sum_{i=1}^{m} \sum_{k=1}^{K}L(y_k,f_{t-1,l}(x))+\sum _{j=1}^{J}c_{jl},I(x_i\in R_{tj})$

由于上式比较难优化，我们用近似值代替

$c_{tjl}=\frac{K-1}{K}=\frac{\sum_{x_i\in R_{tjl}}r_{til}}{\sum _{x_i\in R_{til}}|r_{til}|(1-|r_{til}|)}$

除了负梯度计算和叶子节点的最佳残差拟合的线性搜索，多元GBDT分类和二元GBDT分类以及GBDT回归算法过程相同。

XGBoost

数据建模中，经常采用Boosting方法通过将成百上千个分类准确率较低的树模型组合起来，成为一个准确率很高的预测模型。这个模型会不断地迭代，每次迭代就生成一颗新的树。但在数据集较复杂的时候，可能需要几千次迭代运算，这将造成巨大的计算瓶颈。

针对这个问题。华盛顿大学的陈天奇博士开发的XGBoost（eXtreme Gradient Boosting）基于C++通过多线程实现了回归树的并行构建，并在原有Gradient Boosting算法基础上加以改进，从而极大地提升了模型训练速度和预测精度.

梯度下降与牛顿法

在机器学习任务中，需要最小化损失函数L(θ)，其中θ是要求解的模型参数。梯度下降法常用来求解这种无约束最优化问题，它是一种迭代方法：选取初值$θ^0$ ,不断迭代，更新θ的值，进行损失函数的极小化。

迭代公式：$θ^t = θ_{t-1}+△θ$

梯度下降：将$L(θ^t)$ 在$θ^{t-1}$ 处进行一阶泰勒展开:

$L(θ^t)=L(θ^{t-1}+△θ)\approx L(θ^{t-1})+L'(θ^{t-1})△θ$

要使得$L(θ^t)<L(θ^{t-1})$ ,可取$△θ=-\alpha L’(θ^{t-1})$ ，则$θ^t = θ^{t-1}-\alpha L’(θ^{t-1})$
这里$\alpha$ 是步长，可通过line search确定，但一般直接赋一个小的数

牛顿法：将$L(θ^t)$ 在$θ^{t-1}$ 处进行二阶泰勒展开:

$L(θ^t)\approx L(θ^{t-1})+L'(θ^{t-1})△θ+L''(θ^{t-1})\frac{(△θ)^2}{2}$

可将一阶和二阶导数分别记为g 和 h,则:

$L(θ^t)\approx L(θ^{t-1})+g△θ+h\frac{(△θ)^2}{2}$

要使得$L(θ^t)$极小，即让$g△θ+h\frac{(△θ)^2}{2}$ 极小，可令其对△θ求偏导值为0，求得$△θ=-\frac{g}{h}$ ,故$θ^t = θ^{t-1}+△θ=θ^{t-1}-\frac{g}{h}$ ,将其推广到向量形式，有$θ^t = θ^{t-1}-H^{-1}g$

GBDT 在函数空间中利用梯度下降法进行优化
XGBoost 在函数空间中用牛顿法进行优化

XGBoost的推导过程

定义目标函数

相比原始的GBDT，XGBoost的目标函数多了正则项，使得学习出来的模型更加不容易过拟合。有哪些指标可以衡量树的复杂度？树的深度，内部节点个数，叶子节点个数(T)，叶节点分数(w)…
XGBoost采用的是：T代表叶子数量，w代表叶子预测权值

而XGBoost第t次迭代后，模型的预测等于前t-1次的模型预测加上第t棵树的预测

即每次迭代生成一棵新的回归树，从而使预测值不断逼近真实值（即进一步最小化目标函数）。每一次保留原来的模型不变，加入一个新的函数f到模型里面。其中$\hat{y}_i\left(t-1\right)$就是t-1轮的模型预测，$f_t{(x_i)}$为新t轮加入的预测函数。选取的$f_t{(x_i)}$必须使我们的目标函数尽量最大地降低（这里应用到了Boosting的基本思想，即当前的基学习器重点关注以前所有学习器犯错误的那些数据样本，以此来达到提升的效果）

此时目标损失函数可写作：

如果我们考虑平方误差作为损失函数，公式可改写为：

这里的化简没有看懂，感觉少了一项残差的平方。不过原等式$+f_t{(x_i)}$看做泰勒展开的$+\bigtriangleup x$ .

公式中$y_i,\widehat y_i$ 都是已知的，模型要学习的只有第t棵树$f_t$ .另外对于损失函数不是平方误差的情况，我们可以采用如下的泰勒展开近似来定义一个近似的目标函数，方便我们进行下一步的计算。其中$g_i$和$h_i$为损失函数对$\widehat y_i^{t-1}$的一阶和二阶倒数。

这时候，所有东西都准备好了，最后我们怎么定义$f_t$呢？它可以代表一颗具有预测结果的树，即叶子节点有预测权重。我们定义w为树叶的权重序列，q为树的结构，那么q(x)代表样本x落在树叶的位置。

目标函数的最小化

得到了目标函数，接下来是最关键的一步，在这种新的定义下，我们可以把目标函数进行如下改写，其中$I$被定义为每个叶子上面样本集合$I_j=\{i| q(x_i)=j\}$ 。$L(y_i,\widehat y_i^{t-1})$为真实值与前一个函数计算所得残差是已知的(我们都是在已知前一个树的情况下计算下一颗树的)，同时，在同一个叶子节点上的数的函数值是相同的，可以做合并，于是：

对上诉目标函数求导等于0，可以得到最后的结果：

得到这个目标函数，乍一看目标函数的计算与回归树的结构q函数没有什么关系，但是如果我们仔细回看目标函数的构成，就会发现其中$G_j$和$H_j$的取值是由第j个树叶上数据样本所决定的。而第jj个树上所具有的数据样本则是由树结构q函数决定的。也就是说，一旦回归树的结构q确定，那么相应的目标函数就能够根据上式计算出来。那么回归树的生成问题也就转换为找到一个最优的树结构q，使得它具有最小的目标函数。

计算求得的Obj目标函数代表了当指定一个树的结构的时候，目标函数上面最多减少多少。我们可以把它叫做结构分数（structure score）。可以把它认为是类似于基尼系数一样更加一般的对于树结构进行打分的函数。

当回归树的结构确定时，我们前面已经推导出其最优的叶节点分数以及对应的最小损失值，问题是怎么确定树的结构？才能让得到的结构分数最好，目标函数损失降低最大。主要有以下两种方法

暴力枚举所有可能的树结构，选择损失值最小的 - NP难问题（树的结构有无穷种）
贪心法，每次尝试分裂一个叶节点，计算分裂前后的增益，选择增益最大的（主要）

确定树的结构（贪心法）

分裂前后的增益怎么计算？

ID3算法采用信息增益
C4.5算法采用信息增益比
CART采用Gini系数
XGBoost采用上诉优化函数的打分

即每一次尝试区队已有的叶子加入一个分割。对于一个剧透的分割方案，我们可以获得的增益可以由如下公式计算得到：

这个公式形式上跟ID3算法（采用信息熵计算增益）或者CART算法（采用基尼指数计算增益）是一致的，都是用分裂后的某种值减去分裂前的某种值，从而得到增益。为了限制树的生长，我们可以加入阈值，当增益大于阈值时才让节点分裂，上式中的$\gamma$即阈值，它是正则项里叶子节点数T的系数，所以xgboost在优化目标函数的同时相当于做了预剪枝。另外，上式中还有一个系数$\lambda$，是正则项里leaf score的L2模平方的系数，对leaf score做了平滑，也起到了防止过拟合的作用，这个是传统GBDT里不具备的特性。

但需要注意是：引入的分割不一定会使得情况变好，因为在引入分割的同时也引入新叶子的惩罚项。所以通常需要设定一个阈值，如果引入的分割带来的增益小于一个阀值的时候，我们可以剪掉这个分割。此外在XGBoost的具体实践中，通常会设置树的深度来控制树的复杂度，避免单个树过于复杂带来的过拟合问题。

如何使用及参数

一些常见的问题

1、机器学习算法中GBDT和XGBOOST的区别有哪些？

2、为什么在实际的 kaggle 比赛中 gbdt 和 random forest 效果非常好？

3、为什么xgboost/gbdt在调参时为什么树的深度很少就能达到很高的精度？

这里就不花时间写了，可以参考知乎和一些博客文章。

Bagging和随机森林

如果采样出的每个子集都完全不同，则每个基学习器只用到一小部分训练数据，甚至不足以进行有效学习，这显然无法确保产生出比较好的基学习器。为了解决这个问题，我们可考虑使用相互有交叠的采样子集。

Bagging

随机取出一个样本放入采样集中，再把该样本放回初始数据集。这样的自助采样过程还给 Bagging带来了另一个优点：由于每个基学习器只使用了初始训练集中约 63.2%的样本，剩下约 36.8%的样本可用作验证集来对泛化性能进行包外估计(out-of-bag estimate)。Bagging通常对分类任务使用简单的投票法，对回归任务使用简单平均法。
Bagging的算法描述如下所诉($D_{bs}$是自助采样产生的样本分布,输出采用投票法)：

随机森林

随机森林(Random Forest)是Bagging的一个扩展变体。随机森林在以决策树为基学习器构建 Bagging集成的基础上，进一步在决策树的训练过程中引入了随机属性选择。具体来说，传统决策树在选择划分属性时是在当前结点的属性集合（假定有d个属性）中选择一个最优属性；而在随机森林中，对基决策树的每个结点，先从该结点的属性集合中随机选择一个包含k个属性的子集，然后再从这个子集中选择一个最优属性用于划分。这里的参数k控制了随机性的引入程度：若令 k=d，则基决策树的构建与传统决策树相同；若令k=1,则是随机选择一个属性用于划分；一般情况下，推荐值 k=log2d.随机森林简单、容易实现、计算开销小.

随机森林的训练效率常优于 Bagging，因为在个体决策树的构建过程中，Bagging使用的是“确定型”决策树，在选择划分属性时要对结点的所有属性进行考察，而随机森林使用的“随机型”决策树则只需考察一个属性子集。

由于这些树是随机生成的，大部分的树对解决分类或回归问题是没有意义的，那么生成上万的树有什么好处呢？好处便是生成的决策树中有少数非常好的决策树。当你要做预测的时候，新的观察值随着决策树自上而下的预测并被赋予一个预测值或标签。一旦森林中的每棵树都有了预测值或标签，所有的预测结果将被归总到一起，所有树的投票做为最终的预测结果。简单来说，会像大数原理一样，抛硬币的次数接近无穷，其正反概率会越接近真实概率1/2。大部分的树会相互抵消，最后得到一个泛化较好的结果，从而得到一个好的预测结果。

学习器结合策略

假定集成包含 T 个基学习器h1,h2,…ht,其中hi在示例x上的输出为hi(x)下面介绍几种对hi进行结合的常见策略。

投票法 1）硬投票 2）相对多数投票 3）加权投票
平均法 1）简单平均法 2）加权平均法（较好）
学习法 1）Stacking

Stacking
当训练数据很多时，一种更为强大的结合策略是使用“学习法”，即通过另一个学习器来进行结合。Stacking是学习法的典型代表。这里我们把个体学习器称为初级学习器，用于结合的学习器称为次级学习器或元学习(metalearner)。

Stacking先从初始数据集训练出初级学习器，然后“生成”一个新数据集用于训练次级学习器。在这个新数据集中，初级学习器的输出被当做样例输入特征，而初始样本的标记仍被当作样例标记。Stacking的算法描述如下，这里假定初级学习器使用不同的学习算法产生，即初级集成是异质的（初级学习器也可是同质的）。

train数据是初级训练器5折交叉得到的输出，test是初级训练器预测test后的均值

应用

投票分类

像我们之前讨论的一样，我们会在一个项目快结束的时候使用集成算法，一旦你建立了一些好的分类器，就把他们合并为一个更好的分类器。事实上，在机器学习竞赛中获得胜利的算法经常会包含一些集成方法。

接下来的代码创建和训练了在 sklearn 中的投票分类器。这个分类器由三个不同的分类器组成（训练集是第五章中的 moons 数据集）：

>>> from sklearn.ensemble import RandomForestClassifier 
>>> from sklearn.ensemble import VotingClassifier 
>>> from sklearn.linear_model import LogisticRegression 
>>> from sklearn.svm import SVC
>>> log_clf = LogisticRegression() 
>>> rnd_clf = RandomForestClassifier() 
>>> svm_clf = SVC()
>>> voting_clf = VotingClassifier(estimators=[('lr', log_clf), ('rf', rnd_clf), >>> ('svc', svm_clf)],voting='hard') 
>>> voting_clf.fit(X_train, y_train)

让我们看一下在测试集上的准确率：

>>> from sklearn.metrics import accuracy_score 
>>> for clf in (log_clf, rnd_clf, svm_clf, voting_clf): 
>>>     clf.fit(X_train, y_train) 
>>>     y_pred = clf.predict(X_test) 
>>>     print(clf.__class__.__name__, accuracy_score(y_test, y_pred)) 
LogisticRegression 0.864 
RandomForestClassifier 0.872 
SVC 0.888 
VotingClassifier 0.896

你看！投票分类器比其他单独的分类器表现的都要好。

如果所有的分类器都能够预测类别的概率（例如他们有一个predict_proba()方法），那么你就可以让 sklearn 以最高的类概率来预测这个类，平均在所有的分类器上。这种方式叫做软投票。他经常比硬投票表现的更好，因为它给予高自信的投票更大的权重。你可以通过把voting=”hard”设置为voting=”soft”来保证分类器可以预测类别概率。然而这不是 SVC 类的分类器默认的选项，所以你需要把它的probability hyperparameter设置为True（这会使 SVC 使用交叉验证去预测类别概率，其降低了训练速度，但会添加predict_proba()方法）。如果你修改了之前的代码去使用软投票，你会发现投票分类器正确率高达 91%

Bagging 和 Pasting

就像之前讲到的，可以通过使用不同的训练算法去得到一些不同的分类器。另一种方法就是对每一个分类器都使用相同的训练算法，但是在不同的训练集上去训练它们。有放回采样被称为装袋（Bagging，是 bootstrap aggregating 的缩写）。无放回采样称为粘贴（pasting）

换句话说，Bagging 和 Pasting 都允许在多个分类器上对训练集进行多次采样，但只有 Bagging 允许对同一种分类器上对训练集进行进行多次采样。

当所有的分类器被训练后，集成可以通过对所有分类器结果的简单聚合来对新的实例进行预测。聚合函数通常对分类是统计模式（例如硬投票分类器）或者对回归是平均。每一个单独的分类器在如果在原始训练集上都是高偏差，但是聚合降低了偏差和方差。通常情况下，集成的结果是有一个相似的偏差，但是对比与在原始训练集上的单一分类器来讲有更小的方差。

sklearn 中的 Bagging 和 Pasting

sklearn 为 Bagging 和 Pasting 提供了一个简单的API：BaggingClassifier类（或者对于回归可以是BaggingRegressor。接下来的代码训练了一个 500 个决策树分类器的集成，每一个都是在数据集上有放回采样 100 个训练实例下进行训练（这是 Bagging 的例子，如果你想尝试 Pasting，就设置bootstrap=False）。n_jobs参数告诉 sklearn 用于训练和预测所需要 CPU 核的数量。（-1 代表着 sklearn 会使用所有空闲核）：

>>>from sklearn.ensemble import BaggingClassifier 
>>>from sklearn.tree import DecisionTreeClassifier
>>>bag_clf = BaggingClassifier(DecisionTreeClassifier(), n_estimators=500, 	max_samples=100, bootstrap=True, n_jobs=-1) 
>>>bag_clf.fit(X_train, y_train) 
>>>y_pred = bag_clf.predict(X_test)

如果基分类器可以预测类别概率（例如它拥有predict_proba()方法），那么BaggingClassifier会自动的运行软投票，这是决策树分类器的情况。

下图对比了单一决策树的决策边界和 Bagging 集成 500 个树的决策边界，两者都在 moons 数据集上训练。正如你所看到的，集成的分类比起单一决策树的分类产生情况更好：集成有一个可比较的偏差但是有一个较小的方差（它在训练集上的错误数目大致相同，但决策边界较不规则）。

Bootstrap 在每个预测器被训练的子集中引入了更多的分集，所以 Bagging 结束时的偏差比 Pasting 更高，但这也意味着预测因子最终变得不相关，从而减少了集合的方差。总体而言，Bagging 通常会导致更好的模型，这就解释了为什么它通常是首选的。然而，如果你有空闲时间和 CPU 功率，可以使用交叉验证来评估 Bagging 和 Pasting 哪一个更好。

Out-of-Bag 评价

对于 Bagging 来说，一些实例可能被一些分类器重复采样，但其他的有可能不会被采样。BaggingClassifier默认采样。BaggingClassifier默认是有放回的采样m个实例（bootstrap=True），其中m是训练集的大小，这意味着平均下来只有63%的训练实例被每个分类器采样，剩下的37%个没有被采样的训练实例就叫做 Out-of-Bag 实例。注意对于每一个的分类器它们的 37% 不是相同的。

因为在训练中分类器从来没有看到过 oob 实例，所以它可以在这些实例上进行评估，而不需要单独的验证集或交叉验证。你可以拿出每一个分类器的 oob 来评估集成本身。

在 sklearn 中，你可以在训练后需要创建一个BaggingClassifier来自动评估时设置oob_score=True来自动评估。接下来的代码展示了这个操作。评估结果通过变量oob_score_来显示：

>>> bag_clf = BaggingClassifier(DecisionTreeClassifier(), n_estimators=500,bootstrap=True, n_jobs=-1, oob_score=True)
>>> bag_clf.fit(X_train, y_train) 
>>> bag_clf.oob_score_ 
0.93066666666666664

根据这个 obb 评估，BaggingClassifier可以再测试集上达到93.1%的准确率，让我们修改一下：

>>> from sklearn.metrics import accuracy_score 
>>> y_pred = bag_clf.predict(X_test) 
>>> accuracy_score(y_test, y_pred) 
0.93600000000000005

我们在测试集上得到了 93.6% 的准确率，足够接近了！

对于每个训练实例 oob 决策函数也可通过oob_decision_function_变量来展示。在这种情况下（当基决策器有predict_proba()时）决策函数会对每个训练实例返回类别概率。例如，oob 评估预测第二个训练实例有 60.6% 的概率属于正类（39.4% 属于负类）：

1
2
3

>>> bag_clf.oob_decision_function_ 
array([[ 0.,  1.], [ 0.60588235,  0.39411765],[ 1., 0. ], 
...  [ 1. ,  0. ],[ 0.,  1.],[ 0.48958333,  0.51041667]])

随机森林

正如我们所讨论的，随机森林是决策树的一种集成，通常是通过 bagging 方法（有时是 pasting 方法）进行训练，通常用max_samples设置为训练集的大小。与建立一个BaggingClassifier然后把它放入 DecisionTreeClassifier 相反，你可以使用更方便的也是对决策树优化够的RandomForestClassifier（对于回归是RandomForestRegressor）。接下来的代码训练了带有 500 个树（每个被限制为 16 叶子结点）的决策森林，使用所有空闲的 CPU 核：

>>>from sklearn.ensemble import RandomForestClassifier
>>>rnd_clf = RandomForestClassifier(n_estimators=500, max_leaf_nodes=16, n_jobs=-1) 
>>>rnd_clf.fit(X_train, y_train)
>>>y_pred_rf = rnd_clf.predict(X_test)

除了一些例外，RandomForestClassifier使用DecisionTreeClassifier的所有超参数（决定数怎么生长），把BaggingClassifier的超参数加起来来控制集成本身

随机森林算法在树生长时引入了额外的随机；与在节点分裂时需要找到最好分裂特征相反，它在一个随机的特征集中找最好的特征。它导致了树的差异性，并且再一次用高偏差换低方差，总的来说是一个更好的模型。以下是BaggingClassifier大致相当于之前的randomforestclassifier：

1	>>>bag_clf = BaggingClassifier(DecisionTreeClassifier(splitter="random", max_leaf_nodes=16),n_estimators=500, max_samples=1.0, bootstrap=True, n_jobs=-1)

极端随机树

当你在随机森林上生长树时，在每个结点分裂时只考虑随机特征集上的特征（正如之前讨论过的一样）。相比于找到更好的特征我们可以通过使用对特征使用随机阈值使树更加随机（像规则决策树一样）。

这种极端随机的树被简称为 Extremely Randomized Trees（极端随机树），或者更简单的称为 Extra-Tree。再一次用高偏差换低方差。它还使得 Extra-Tree 比规则的随机森林更快地训练，因为在每个节点上找到每个特征的最佳阈值是生长树最耗时的任务之一。

你可以使用 sklearn 的ExtraTreesClassifier来创建一个 Extra-Tree 分类器。他的 API 跟RandomForestClassifier是相同的，相似的， ExtraTreesRegressor 跟RandomForestRegressor也是相同的 API。

我们很难去分辨ExtraTreesClassifier和RandomForestClassifier到底哪个更好。通常情况下是通过交叉验证来比较它们（使用网格搜索调整超参数）

特征重要度

最后，如果你观察一个单一决策树，重要的特征会出现在更靠近根部的位置，而不重要的特征会经常出现在靠近叶子的位置。因此我们可以通过计算一个特征在森林的全部树中出现的平均深度来预测特征的重要性。sklearn 在训练后会自动计算每个特征的重要度。你可以通过feature_importances_变量来查看结果。例如如下代码在 iris 数据集（第四章介绍）上训练了一个RandomForestClassifier模型，然后输出了每个特征的重要性。看来，最重要的特征是花瓣长度（44%）和宽度（42%），而萼片长度和宽度相对比较是不重要的（分别为 11% 和 2%）：

>>> from sklearn.datasets import load_iris 
>>> iris = load_iris() 
>>> rnd_clf = RandomForestClassifier(n_estimators=500, n_jobs=-1) 
>>> rnd_clf.fit(iris["data"], iris["target"]) 
>>> for name, score in zip(iris["feature_names"], rnd_clf.feature_importances_): 
>>>     print(name, score) 
sepal length (cm) 0.112492250999
sepal width (cm) 0.0231192882825 
petal length (cm) 0.441030464364 
petal width (cm) 0.423357996355

相似的，如果你在 MNIST 数据及上训练随机森林分类器（在第三章上介绍），然后画出每个像素的重要性，你可以得到下图

随机森林可以非常方便快速得了解哪些特征实际上是重要的，特别是你需要进行特征选择的时候

提升

下图显示连续五次预测的 moons 数据集的决策边界（在本例中，每一个分类器都是高度正则化带有 RBF 核的 SVM）。第一个分类器误分类了很多实例，所以它们的权重被提升了。第二个分类器因此对这些误分类的实例分类效果更好，以此类推。右边的图代表了除了学习率减半外（误分类实例权重每次迭代上升一半）相同的预测序列（误分类的样本权重提升速率即学习率）。你可以看出，序列学习技术与梯度下降很相似，除了调整单个预测因子的参数以最小化代价函数之外，AdaBoost 增加了集合的预测器，逐渐使其更好。

klearn 通常使用 Adaboost 的多分类版本 SAMME（这就代表了分段加建模使用多类指数损失函数）。如果只有两类别，那么 SAMME 是与 Adaboost 相同的。如果分类器可以预测类别概率（例如如果它们有predict_proba()），如果 sklearn 可以使用 SAMME 叫做SAMME.R的变量（R 代表“REAL”），这种依赖于类别概率的通常比依赖于分类器的更好。

接下来的代码训练了使用 sklearn 的AdaBoostClassifier基于 200 个决策树桩 Adaboost 分类器（正如你说期待的，对于回归也有AdaBoostRegressor）。一个决策树桩是max_depth=1的决策树-换句话说，是一个单一的决策节点加上两个叶子结点。这就是AdaBoostClassifier的默认基分类器：

1
2
3

>>>from sklearn.ensemble import AdaBoostClassifier
>>>ada_clf = AdaBoostClassifier(DecisionTreeClassifier(max_depth=1), n_estimators=200,algorithm="SAMME.R", learning_rate=0.5) 
>>>ada_clf.fit(X_train, y_train)

如果你的 Adaboost 集成过拟合了训练集，你可以尝试减少基分类器的数量或者对基分类器使用更强的正则化。

梯度提升

另一个非常著名的提升算法是梯度提升。与 Adaboost 一样，梯度提升也是通过向集成中逐步增加分类器运行的，每一个分类器都修正之前的分类结果。然而，它并不像 Adaboost 那样每一次迭代都更改实例的权重，这个方法是去使用新的分类器去拟合前面分类器预测的残差。

让我们通过一个使用决策树当做基分类器的简单的回归例子（回归当然也可以使用梯度提升）。这被叫做梯度提升回归树（GBRT，Gradient Tree Boosting 或者 Gradient Boosted Regression Trees）。首先我们用DecisionTreeRegressor去拟合训练集（例如一个有噪二次训练集）：

1
2
3

>>>from sklearn.tree import DecisionTreeRegressor 
>>>tree_reg1 = DecisionTreeRegressor(max_depth=2) 
>>>tree_reg1.fit(X, y)

现在在第一个分类器的残差上训练第二个分类器：

1
2
3

>>>y2 = y - tree_reg1.predict(X) 
>>>tree_reg2 = DecisionTreeRegressor(max_depth=2) 
>>>tree_reg2.fit(X, y2)

随后在第二个分类器的残差上训练第三个分类器：

1
2
3

>>>y3 = y2 - tree_reg1.predict(X) 
>>>tree_reg3 = DecisionTreeRegressor(max_depth=2) 
>>>tree_reg3.fit(X, y3)

现在我们有了一个包含三个回归器的集成。它可以通过集成所有树的预测来在一个新的实例上进行预测。

1	>>>y_pred = sum(tree.predict(X_new) for tree in (tree_reg1, tree_reg2, tree_reg3))

下图左栏展示了这三个树的预测，在右栏展示了集成的预测。在第一行，集成只有一个树，所以它与第一个树的预测相似。在第二行，一个新的树在第一个树的残差上进行训练。在右边栏可以看出集成的预测等于前两个树预测的和。相同的，在第三行另一个树在第二个数的残差上训练。你可以看到集成的预测会变的更好。

我们可以使用 sklean 中的GradientBoostingRegressor来训练 GBRT 集成。与RandomForestClassifier相似，它也有超参数去控制决策树的生长（例如max_depth，min_samples_leaf等等），也有超参数去控制集成训练，例如基分类器的数量（n_estimators）。接下来的代码创建了与之前相同的集成：

1
2
3

>>>from sklearn.ensemble import GradientBoostingRegressor
>>>gbrt = GradientBoostingRegressor(max_depth=2, n_estimators=3, learning_rate=1.0) 
>>>gbrt.fit(X, y)

超参数learning_rate 确立了每个树的贡献。如果你把它设置为一个很小的树，例如 0.1，在集成中就需要更多的树去拟合训练集，但预测通常会更好。这个正则化技术叫做 shrinkage。下图展示了两个在低学习率上训练的 GBRT 集成：其中左面是一个没有足够树去拟合训练集的树，右面是有过多的树过拟合训练集的树。

为了找到树的最优数量，你可以使用早停技术。最简单使用这个技术的方法就是使用staged_predict()：它在训练的每个阶段（用一棵树，两棵树等）返回一个迭代器。加下来的代码用 120 个树训练了一个 GBRT 集成，然后在训练的每个阶段验证错误以找到树的最佳数量，最后使用 GBRT 树的最优数量训练另一个集成

>>>import numpy as np 
>>>from sklearn.model_selection import train_test_split
>>>from sklearn.metrics import mean_squared_error

>>>X_train, X_val, y_train, y_val = train_test_split(X, y)
>>>gbrt = GradientBoostingRegressor(max_depth=2, n_estimators=120) 
>>>gbrt.fit(X_train, y_train)
>>>errors = [mean_squared_error(y_val, y_pred)         
     for y_pred in gbrt.staged_predict(X_val)] 
>>>bst_n_estimators = np.argmin(errors)
>>>gbrt_best = GradientBoostingRegressor(max_depth=2,n_estimators=bst_n_estimators) 
>>>gbrt_best.fit(X_train, y_train)

验证错误在图的左面展示，最优模型预测被展示在右面

你也可以早早的停止训练来实现早停（与先在一大堆树中训练，然后再回头去找最优数目相反）。你可以通过设置warm_start=True来实现，这使得当fit()方法被调用时 sklearn 保留现有树，并允许增量训练。接下来的代码在当一行中的五次迭代验证错误没有改善时会停止训练：

>>>gbrt = GradientBoostingRegressor(max_depth=2, warm_start=True)
min_val_error = float("inf") 
error_going_up = 0 
for n_estimators in range(1, 120):    
    gbrt.n_estimators = n_estimators    
    gbrt.fit(X_train, y_train)    
    y_pred = gbrt.predict(X_val)    
    val_error = mean_squared_error(y_val, y_pred)    
    if val_error < min_val_error:        
        min_val_error = val_error        
        error_going_up = 0    
    else:        
        error_going_up += 1        
        if error_going_up == 5:            
            break  # early stopping

GradientBoostingRegressor也支持指定用于训练每棵树的训练实例比例的超参数subsample。例如如果subsample=0.25，那么每个树都会在 25% 随机选择的训练实例上训练。你现在也能猜出来，这也是个高偏差换低方差的作用。它同样也加速了训练。这个技术叫做随机梯度提升。

也可能对其他损失函数使用梯度提升。这是由损失超参数控制（见 sklearn 文档）。

Stacking

本章讨论的最后一个集成方法叫做 Stacking（stacked generalization 的缩写）。这个算法基于一个简单的想法：不使用琐碎的函数（如硬投票）来聚合集合中所有分类器的预测，我们为什么不训练一个模型来执行这个聚合？图展示了这样一个在新的回归实例上预测的集成。底部三个分类器每一个都有不同的值（3.1，2.7 和 2.9），然后最后一个分类器（叫做 blender 或者 meta learner ）把这三个分类器的结果当做输入然后做出最终决策（3.0）

练习题

如果你在相同训练集上训练 5 个不同的模型，它们都有 95% 的准确率，那么你是否可以通过组合这个模型来得到更好的结果？如果可以那怎么做呢？如果不可以请给出理由。
软投票和硬投票分类器之间有什么区别？
是否有可能通过分配多个服务器来加速 bagging 集成系统的训练？pasting 集成，boosting 集成，随机森林，或 stacking 集成怎么样？
out-of-bag 评价的好处是什么？
是什么使极端随机树Extra-Tree 比规则随机森林更随机呢？这个额外的随机有什么帮助呢？那这个 Extra-Tree 比规则随机森林谁更快呢？
如果你的 Adaboost 模型欠拟合，那么你需要怎么调整超参数？
如果你的梯度提升过拟合，那么你应该调高还是调低学习率呢？

1、只要模型间多样性较大，组合成一个集合模型，会起到一定的效果。

2、硬投票分类器只计算集成中每个分类器的投票数，并选择得票最多的类。软投票分类器计算每个类的平均估计概率，并选择具有最高概率的类。软投票使概率大的类别权重更高，通常表现更好，但只有在可以估计类概率的分类器时才有效（例如，对于Scikit-Learn中的SVM分类器，您必须设置probability = True）。

3、 bagging，pasting，随机森林是可以的。boosting 集成因为学习器需要基于先前的学习器构建，因此训练是连续的，故不适合。至于stacking 集成，给定层中的所有预测变量彼此独立，因此可以在多个服务器上并行训练它们。但是，一层中的预测变量只能在前一层中的预测变量都经过训练后才能进行训练。

4、未被训练过得实例可以当做验证集

5、在随机森林中，只考虑特征的随机子集并在每个节点处进行分割。对于Extra-Trees也是如此，但它们更进一步：不是像常规决策树一样搜索最佳阈值，而是为每个特征使用随机阈值。这种额外的随机性就像一种正则化的形式：如果随机森林过度拟合训练数据，极端随机树可能表现更好。此外，由于Extra-Trees不会搜索最佳阈值，因此它们比随机森林训练要快得多。然而，在做出预测时，它们既不比随机森林更快也不慢。

6、如果您的AdaBoost集合欠拟合，可以尝试增加学习器的数量或减少基本学习器的正则化超参数。也可以尝试稍微提高学习率。

7、如果您的Gradient Boosting过拟合，应该尝试降低学习率。你也可以使用早期停止法。

Sklearn 与 TensorFlow 机器学习实用指南（五）：决策树

发表于 2018-07-23 | 分类于 Sklearn 与 TensorFlow 机器学习实用指南

字数统计: 3.4k | 阅读时长 ≈ 13

一些基本概念

决策树简述

决策树（Decision Tree）是数据挖掘中一种基本的分类和回归方法，它呈树形结构，在分类问题中，表示基于特征对实例进行分类的过程，可以认为是if−thenif−then规则的集合。决策树模型的主要优点是模型具有可读性，分类速度快。在学习时，利用训练数据，根据损失函数最小化原则建立决策树模型；而在预测时，对新的数据，利用决策树模型进行分类。主要的决策树算法有ID3算法、C4.5算法和CART算法。一个决策树的学习过程包括三个步骤：特征选择、决策树的生成以及决策树的修剪。

信息熵

信息熵是衡量样本纯度的一种指标，嘉定当前样本集合D中第k类样本所占的比例为$p_k(k=1,2,…,|y|)$,则D的信息熵定义为

$Ent(D)=-\sum_{k=1}^{\vert y\vert}p_klog_2{p_k}$

Ent(D)的值越小，则D的纯度越高。

以周志华西瓜书P76的西瓜数据集中的17个样本数据为例子。

显然，标签类别|y|=2.其中正例占$p_1$=8/17，反例占$p_2$=9/17。于是根节点的信息熵为

$Ent(D)=-\sum_{k=1}^2p_klog_2^{p_k}=-(\frac{8}{17}log_2\frac{8}{17}+\frac{9}{17}log_2\frac{9}{17})=0.998$

条件熵

$Ent(Y|X)$表示在已知随机变量X的条件下随机变量Y的不确定性，定义为X给定条件下Y的条件概率分布的熵对X的数学期望

$Ent\left( Y|X \right) =\sum_{i=1}^n{p_iEnt\left( Y|X=x_i \right)}$

特征选择

信息增益

信息增益表示得知特征X的信息而使得类Y的信息的不确定性减少的程度。特征A对训练数据集D的信息增益Gain(D,a)定义为集合D的经验熵$Ent(D)$与特征A给定条件下D的经验条件熵$Ent(Y | X)$之差，即假定离散属性a有n个可能的取值${a^1,a^2,…a^n}$，若使用属性a的取值来对样本集合划分，则会产生n个分支节点子集$D_1,D_2,..,Dn$，$|D_i|$ 为$D_i$的样本个数。(考虑不同的分支节点所包含的样本数不同，用$\frac{|D^i|}{|D|}$代替期望概率$p_i$，即样本越多的分支节点的影响越大)

$Gain(D,a)=Ent(D)-\sum_{i=1}^n\frac{|D^i|}{|D|}End(D^i)$

以西瓜数据集的“色泽”属性为例，它有3个可能的取值{青绿，乌黑，浅白}。若使用该属性对D角线划分，则可以得到三个子集，分布记为$D^1$(色泽=青绿)={1,4,6,10,13,17}，正反例率分别为$p_1=\frac{3}{6}$，$p_2=\frac{3}{6}$；$D^2$(色泽=乌黑)={2,3,7,8,9,15}，正反例率分别为$p_1=\frac{4}{6}$，$p_2=\frac{2}{6}$ ；$D^3$(色泽=乌黑)={5,11,12,14,16}，正反例率分别为$p_1=\frac{1}{5}$，$p_2=\frac{4}{5}$。

首先，根据信息熵公式可计算出“色泽”划分之后所获得的3个分支节点的信息熵为

$Ent(D^1)=-(\frac{3}{6}log_2\frac{3}{6}+\frac{3}{6}log_2\frac{3}{6})=1.000 \\ \\ \\ Ent(D^2)=-(\frac{4}{6}log_2\frac{4}{6}+\frac{2}{6}log_2\frac{2}{6})=0.918 \\ \\ \\ Ent(D^3)=-(\frac{1}{5}log_2\frac{1}{5}+\frac{4}{5}log_2\frac{4}{5})=0.722$

之后根据上面公式可计算出属性“色泽”的信息增益为

$Gain(D,色泽)=Ent(D)-\sum_{i=1}^3\frac{|D^i|}{|D|}Ent(D^i) \\ \\ \\ =0.998-((\frac{6}{17}\times1.000+\frac{6}{17}\times0.918+\frac{5}{17}\times0.722)=0.109$

一般而言，信息增益越大，则意味着使用属性a来进行划分所获得的“纯度提升”越大。其中ID3决策树就是以信息增益为准则来选择划分属性。根据得到的属性最大信息增益来划分结果示例如下所示。之后再根据子集样本进一步划分，直到只有一个样本个体或者样本个体都为同一类标签为止。

若我们把样本编号1-17也作为一个划分属性，则根据信息增益公式可计算出它的信息增益为0.998，远大于其他候选划分属性。这很容易理解，“编号”将产生17个分支，每个分支仅包含一个样本，这些分支节点纯度已打最大。然而这些侧介绍显然不具有泛化能力。

所以实际上，信息增益准则对可取数目较多的属性有所偏好；为了减少这种偏好，著名的C4.5决策树算法不直接采用信息增益，而是使用增益率来做最优划分属性。

增益率

增益率定义为其信息增益Gain(D，a)与训练数据集D关于特征a样本数的信息熵之比，增益率的数学定义为

$Gain_{ratio}(D,a)=\frac{Gain(D,a)}{IV(a)}$ $IV(a)=-\sum_{v=1}^n\frac{|D^i|}{|D|}log_2\frac{|D^i|}{|D|}$

（注意子集的划分是属性a样本的可取类别数目n，而不是标签类别|y|了啊），IV(a)可以成为属性a的内在信息，若属性a的可取数目越大，则IV(a)的值通常越大。这样就可以一定程度平衡信息增益对可取数目较多的属性的偏好。以“色泽”为例

$IV(色泽)=-\frac{6}{17}log_2\frac{6}{17}-\frac{6}{17}log_2\frac{6}{17}-\frac{5}{17}log_2\frac{5}{17}=1.580$

另外，需要注意的是，增益率准则对可取数值数目较少的属性有所偏好，因此，C4.5算法并不是直接选择增益率最大的候选划分属性，而是先从候选划分属性中找出信息增益高于平均水平的属性，再从中选择增益率的最高的。

基尼指数

数据集D的纯度可以用基尼指数来度量，Gini（D）越小，则数据集D的纯度越高。假设有K个类，样本点属于第k类的概率为$p_k$，则概率分布的基尼系数定义为

$Gini\left( p \right) =\sum_{k=1}^K{p_k\left( 1-p_k \right) =1-\sum_{k=1}^K{p_k^2}}$

根据基尼指数定义，可以得到样本集合D的基尼指数，其中$D_k$表示数据集D中属于第k类的样本子集

$Gini(D)=1-\sum_{k=1}^{K}\left(\frac{|D_k|}{|D|} \right)^2$

若样本集合D根据特征A是否取某一可能值a被分割成$D_1$和$D_2$两部分，即

$D_1=\left\{ \left( x,y \right) \in D|A\left( x \right) =0 \right\} \mathrm{，}D_2=D-D_1$

则在特征A的条件下，集合D的基尼指数定义为

$Gini\left( D,A \right) =\frac{|D_1|}{|D|}Gini\left( D_1 \right) +\frac{|D_2|}{|D|}Gini\left( D_2 \right)$

基尼系数Gini(D)表示集合D的不确定性，基尼系数Gini(D,A)表示A=a分割后集合D的不确定性。基尼指数越大，样本集合的不确定性越大。对于属性A，分别计算任意属性值将数据集划分为两部分之后的Gain_Gini，选取其中的最小值，作为属性A得到的最优二分方案。然后对于训练集S，计算所有属性的最优二分方案，选取其中的最小值，作为样本及S的最优二分方案。

$\min_{i\epsilon A}(Gain\_Gini(D,A))$ $\min_{A\epsilon Attribute}(\min_{i\epsilon A}(Gain\_Gini(D,A)))$

剪枝

剪枝是决策树学习算法对付“过拟合”的主要手段。在决策树学习中，为了尽可能正确分类样本，节点划分过程不断重复，有时会造成决策树分支过多，这时有时候把自身特点当做所有数据都具有的一般性质而导致过拟合，因此，可以通过主动去掉一些分支来降低过拟合的风险

基本策略有预剪枝和后剪枝，预剪枝是指在决策树生成过程中，对每个节点在划分前先进行估计，若当前节点的划分不能带来决策树泛化性能的提升，则停止划分当前节点标记为叶节点；后剪枝则是先从训练集生成一颗完整的决策树，然后自底向上地对非叶节点进行考察，若将该节点对应的字数替换为叶节点能带来决策树泛化性能提升，则将该子树替换为叶节点。

CART决策树的生成

ID3算法和C4.5按照信息增益和增益率最大的特征作为节点特征，然后递归构建。比较简单，还要注意这两种方法都容易过拟合。这里主要讲一下CART分类与回归树。

分类树与回归树（classification and regression tree，CART）模型（Breiman）由特征选择、树生成及剪枝组成，既可用于分类也可用于回归。CART算法采用二分递归分割的技术将当前样本集分为两个子样本集，使得生成的每个非叶子节点都有两个分支。因此CART算法生成的决策树是结构简洁的二叉树。CART可以处理连续型变量和离散型变量，利用训练数据递归的划分特征空间进行建树，用验证数据进行剪枝。利用训练数据递归的划分特征空间进行建树，用验证数据进行剪枝。

如果待预测分类是离散型数据，则CART生成分类决策树。
如果待预测分类是连续性数据，则CART生成回归决策树。

CART分类树

对分类树用基尼系数（Gini index）最小化准则，进行特征选择，生成二叉树。

具体算法步骤如下：

1）设结点的训练数据集为D，计算现有特征对该数据集的基尼指数。此时，对每一个特征A，对其可能取的每个值a，根据样本点对A=a的测试为”是”或者“否”将D分割为D1和D2两部分，计算其基尼系数。
2）在所有可能的特征A以及他们所有可能的切分点a中，选择基尼系数最小的特征及其对应的切分点作为最优特征与最优切分点。依最优特征与最优切分点，从现结点生成两个子结点，将训练数据集依特征分配到两个子结点中去。
3）对两个子结点递归地调用上述两个步骤，直至满足停止条件。
4）生成CART决策树

以生物特征分类数据为例

针对上述离散型数据，按照体温为恒温和非恒温(多类别属性也按’非’分成两类别)进行划分。其中恒温时包括哺乳类5个、鸟类2个，非恒温时包括爬行类3个、鱼类3个、两栖类2个，如下所示我们计算以体温为划分特征D1,D2的基尼指数。

$Gini(D_1)=1-[ (\frac{5}{7})^2+(\frac{2}{7})^2]=\frac{20}{49}$ $Gini(D_2)=1-[ (\frac{3}{8})^2+(\frac{3}{8})^2+(\frac{2}{8})^2]=\frac{42}{64}$

然后计算得到特征体温下数据集的Gini指数，最后我们选择Gain_Gini最小的特征和相应的划分

$Gain\_Gini(D,体温)=\frac{7}{15}*\frac{20}{49}+\frac{8}{15}*\frac{42}{64}$

在所有可能的特征A以及他们所有可能的切分点a中的二分划分中，选择基尼系数最小的特征及其对应的切分点作为最优特征与最优切分点，依次递归。

CART回归树

回归树衡量最好的标准不再是最大熵，而是最小化均方差。而且在每个节点（不一定是叶子节点）都会得一个预测值，这个预测值可为所有样本的平均值。我们利用最小二乘回归树生成算法来生成回归树f(x)，即在训练数据集所在的输入空间中，递归地将每个区域分为两个子区域并决定每个子区域上的输出值，构建二叉决策树。

回归树算法流程：j为选定的某个特征属性，s为在这个特征上的切分数值（需要遍历所有特征和切分点来选到最小化均方差），R1R2为切分的样本，C1C2为切分区域样本中的特征均值。

实例详解：

考虑如上所示的连续性变量，根据给定的数据点，考虑1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5共9个自定均分的切分点。对各切分点依次求出R1,R2,c1,c2及m(s)，例如当切分点s=1.5时，得到R1={1},R2={2,3,4,5,6,7,8,9,10}，其中c1,c2,m(s)如下所示

$c_1=\frac{1}{N_m}\sum_{x_i\epsilon R_m(j,s)}y_i=\frac{1}{1}\sum_{x_i\epsilon R_1(1,1.5)}5.56=5.56$ $c_2=\frac{1}{N_m}\sum_{x_i\epsilon R_m(j,s)}y_i=\frac{1}{9}\sum_{x_i\epsilon R_2(1,1.5)}(5.70+5.91+...+9.05)=7.50$ $m(s)=\min_{j,s}[\min_{c_1}\sum _{x_i\epsilon R_i(j,s)}(y_i-c_1)^2+\min_{c_2}\sum _{x_i\epsilon R_i(j,s)}(y_i-c_1)^2]=0+15.72=15.72$

依次改变(j,s)对，可以得到s及m(s)的计算结果，如下表所示。

当x=6.5时，此时R1={1,2,3,4,5,6},R2={7,8,9,10},c1=6.24,c2=8.9。回归树T1(x)为

$T_1(x)=\begin{cases} & 6.24,x<6.5 \\ & 8.91,x\ge 6.5 \end{cases}$ $f_1(x)=T_1(x)$

然后我们利用f1(x)拟合训练数据的残差(各目标值减去对应的c1,c2)，如下表所示

用f1(x)拟合训练数据得到平方误差

$L(y,f_1(x))=\sum_{i=1}^{10}(y_i-f_1(x_i))^2=1.93$

第二步求T2(x)与求T1(x)方法相同，只是拟合的数据是上表的残差。可以得到

$T_2(x)=\begin{cases} & -0.52,x<3.5 \\ & 0.22,x\ge 3.5 \end{cases}$ $f_2(x)=f_1(x)+T_2(x)= \begin{cases} & 5.72,x<3.5 \\ & 6.46,3.5\le x \le 6.5 \\ & 9.13,x\ge 6.5 \end{cases}$

用f2(x)拟合训练数据的平方误差

$L(y,f_2(x))=\sum_{i=1}^{10}(y_i-f_2(x_i))^2=0.79$

继续求得T3(x)、T4(x)、T5(x)、T6(x)，如下所示

$T_3(x)=\begin{cases} & 0.15,x<6.5 \\ & -0.22,x\ge 6.5 \end{cases} L(y,f_3(x))=0.47$ $T_4(x)=\begin{cases} & -0.16,x<4.5 \\ & 0.11,x\ge 4.5 \end{cases} L(y,f_4(x))=0.30$ $T_5(x)=\begin{cases} & 0.07,x<6.5 \\ & -0.11,x\ge 6.5 \end{cases} L(y,f_5(x))=0.23$ $T_6(x)=\begin{cases} & -0.15,x<2.5 \\ & 0.04,x\ge 2.5 \end{cases}$ $f_6(x)=f_5(x)+T_6(x)=T_1(x)+...+T_6(x)= \begin{cases} & 5.63,x<2.5 \\ & 5.82,2.5\le x \le 3.5 \\ & 6.56,3.5\le x \le 4.5 \\ & 6.83,4.5\le x \le 6.5 \\ & 8.95,x\ge 6.5 \end{cases}$

用f6(x)拟合训练数据的平方损失误差如下所示，假设此时已经满足误差要求，那么f(x)=f6(x)便是所求的回归树。

$L(y,f_6(x))=\sum_{i=1}^{10}(y_i-f_6(x_i))^2=0.71$

CART剪枝

们将一颗充分生长的树称为T0 ，希望减少树的大小来防止过拟化。但同时去掉一些节点后预测的误差可能会增大，即完整树T0拟合的损失是最小的。那么如何达到这两个变量之间的平衡则是问题的关键。因此我们用一个变量α 来平衡，定义损失函数如下

$C_\alpha(T)=C(T)+\alpha|T|$

T为任意子树，|T|为子树T的叶子节点个数。
α是参数，权衡拟合程度与树的复杂度。
C(T)为预测误差，可以是平方误差也可以是基尼指数，C(T)衡量训练数据的拟合程度。

那么我们如何找到这个合适的α来使拟合程度与复杂度之间达到最好的平衡呢？准确的方法就是将α从0取到正无穷，对于每一个固定的α，我们都可以找到使得$C_\alpha(T)$最小的最优子树T(α)。

当α很小的时候(相当于弱化正则)，T0 完整树是这样的最优子树.
当α很大的时候(相当于强化正则)，单独一个根节点就是最优子树。

Breiman等人证明：可以用递归地方法对树进行剪枝。将a从小增大，$0=a_0<a_1<…..a_n<+∞$产生一系列的区间$[a_i,a_i+1),i=0,1,…,n$剪枝得到的子树序列对应着区间$a∈[a_i,a_i+1)，i=0,1,2,…,n$的最优子树序列为${T_0,T_1,T_2,…,T_n}$,序列的子树是嵌套的。

具体地，从整体树$T_0$开始剪枝，对$T_0$的人以内部结点t，以t为单结点树的损失函数是

$C_a\left( t \right) =C\left( t \right) +a$

以t为根结点的子树$T_t$的损失函数是

$C_a(T_t)=C(T_t)+a|T_t|$

到这里，我们的目的就变得很明确了，当t为单节点树的损失比以t为根节点的子树$T_t$损失相等的时候，损失相同，而t的节点少，就进行剪枝操作。其过程会随a由0逐渐增大，单节点树t损失一开始大于子树$T_t$，随后不断减少，直到大于子树损失。具体的，当$\alpha =\frac{C\left( t \right) -C\left( T_t \right)}{|T_t|-1}$ ，$T_t$和t有相同的损失。

为此，对$T_0$中的每一个内部结点t，计算

$g\left( t \right) =\frac{C\left( t \right) -C\left( T_t \right)}{|T_t|-1}$

它表示剪枝后整体损失函数减少的程度。在$T_0$中减去g(t)最小的$T_t$，将得到的子树作为$T_1$，同时将最小的g(t)设为$a_1$，$T_1$为区间$[a_1,a_2)$的最优子树。如此剪枝下去，直至得到根结点。在这一过程中，不断得增加a的值，产生新的区间。

然后，在剪枝得到的子树序列$T_0,T_1,…,T_n$中通过交叉验证选取最优子树$T_a$.

具体地，利用独立的验证数据集，测试子树序列$T_0,T_1,…,T_n$中各棵子树的平方误差或基尼指数。平方误差或基尼指数最小的决策树被认为是最优的决策树。在子树序列中，每棵子树$T_0,T_1,…,T_n$都对应一个参数$a_1,a_2,…,a_n$。所以当最优子树$T_k$确定时，对应的$a_k$也就确定了，即得到最由决策树$T_a$.

应用

决策树分类边界

另外，决策树所形成的分类边界有一个明显的特点：轴平行，即它的分类边界有若干个与坐标轴平行的分段组成。

决策树的训练和可视化

下面的代码就是在我们熟知的鸢尾花数据集上进行一个决策树分类器的训练

决策树的众多特性之一就是，它不需要太多的数据预处理，尤其是不需要进行特征的缩放或者归一化。

from sklearn.datasets import load_iris
from sklearn.tree import DecisionTreeClassifier
iris = load_iris()
X = iris.data[:, 2:] # petal length and width y = iris.target
    tree_clf = DecisionTreeClassifier(max_depth=2)
    tree_clf.fit(X, y)

你可以通过使用export_graphviz()方法，通过生成一个叫做iris_tree.dot的图形定义文件将一个训练好的决策树模型可视化

from sklearn.tree import export_graphviz
    export_graphviz(
            tree_clf,
            out_file=image_path("iris_tree.dot"),
            feature_names=iris.feature_names[2:],
            class_names=iris.target_names,
            rounded=True,
            filled=True
        )

然后，我们可以利用graphviz package 中的dot命令行，将.dot文件转换成 PDF 或 PNG 等多种数据格式。例如，使用命令行将.dot文件转换成.png文件的命令如下：

Graphviz是一款开源图形可视化软件包，http://www.graphviz.org/

1	$ dot -Tpng iris_tree.dot -o iris_tree.png

我们的第一个决策树如图

节点的samples属性统计出它应用于多少个训练样本实例，例如，我们有一百个训练实例是花瓣长度大于 2.45 里面的（深度为 1，右侧），在这 100 个样例中又有 54 个花瓣宽度小于 1.75cm（深度为 2，左侧）；节点的value属性告诉你这个节点对于每一个类别的样例有多少个，例如：右下角的节点中包含 0 个 Iris-Setosa，1 个 Iris-Versicolor 和 45 个 Iris-Virginica；最后，节点的Gini属性用于测量它的纯度：如果一个节点包含的所有训练样例全都是同一类别的，我们就说这个节点是纯的（Gini=0）

下面公式显示了训练算法如何计算第i个节点的 gini 分数。例如，深度为 2 的左侧节点基尼指数为：

$G_i=\sum_{k=1}^nP_{i,k}^2$

$p_{i,k}$ 是第i个节点中训练实例为的k类实例的比例

Scikit-Learn 用的是 CART 算法， CART 算法仅产生二叉树：每一个非叶节点总是只有两个子节点（只有是或否两个结果）。然而，像 ID3 这样的算法可以产生超过两个子节点的决策树模型

下图显示了决策树的决策边界。粗的垂直线代表根节点（深度为 0）的决定边界：花瓣长度为 2.45 厘米。由于左侧区域是纯的（只有 Iris-Setosa），所以不能再进一步分裂。然而，右边的区域是不纯的，所以深度为 1 的右边节点在花瓣宽度为 1.75 厘米处分裂（用虚线表示）。又由于max_depth设置为 2，决策树在那里停了下来。但是，如果将max_depth设置为 3，两个深度为 2 的节点，每个都将会添加另一个决策边界（用虚线表示）

估计分类概率

决策树还可以估计某个实例属于特定类k的概率：首先遍历树来查找此实例的叶节点，然后它返回此节点中类k的训练实例的比例。

例如，假设你发现了一个花瓣长 5 厘米，宽 1.5 厘米的花朵。相应的叶节点是深度为 2 的左节点，因此决策树应该输出以下概率：Iris-Setosa 为 0%（0/54），Iris-Versicolor 为 90.7%（49/54），Iris-Virginica 为 9.3%（5/54）。当然，如果你要求它预测具体的类，它应该输出 Iris-Versicolor（类别 1），因为它具有最高的概率。我们了测试一下

>>> tree_clf.predict_proba([[5, 1.5]])
array([[ 0. , 0.90740741, 0.09259259]])
>>> tree_clf.predict([[5, 1.5]])
array([1])

CART 训练算法

Scikit-Learn 用分裂回归树（Classification And Regression Tree，简称 CART）算法训练决策树（也叫“增长树”）。这种算法思想真的非常简单：

首先使用单个特征k和阈值$t_k$ (例如，“花瓣长度≤2.45cm”）将训练集分成两个子集。它如何选择k和$t_k$ 呢？它寻找到能够产生最纯粹的子集一对(k,$t_k$) ，然后通过子集大小加权计算.

算法会尝试最小化成本函数。方法如公式

当它成功的将训练集分成两部分之后，它将会继续使用相同的递归式逻辑继续的分割子集，然后是子集的子集。当达到预定的最大深度之后将会停止分裂（由max_depth超参数决定），或者是它找不到可以继续降低不纯度的分裂方法的时候。几个其他超参数（之后介绍）控制了其他的停止生长条件（min_samples_split，min_samples_leaf，min_weight_fraction_leaf，max_leaf_nodes）。

正如您所看到的，CART 算法是一种贪婪算法：它贪婪地搜索最高级别的最佳分割方式，然后在每个深度重复该过程。它不检查分割是否能够在几个级别中的全部分割可能中找到最佳方法。贪婪算法通常会产生一个相当好的解决方法，但它不保证这是全局中的最佳解决方案

计算复杂度

在建立好决策树模型后，做出预测需要遍历决策树，从根节点一直到叶节点。决策树通常近似左右平衡，因此遍历决策树需要经历大致 $O(log_2^m)$ 个节点。由于每个节点只需要检查一个特征的值，因此总体预测复杂度仅为$O(log_2^m)$ ，与特征的数量无关。所以即使在处理大型训练集时，预测速度也非常快.

然而，训练算法的时候（训练和预测不同）需要比较所有特征（如果设置了max_features会更少一些）
在每个节点的所有样本上。就有了$O(nmlog_2^m)$的训练复杂度。对于小型训练集（少于几千例），Scikit-Learn 可以通过预先设置数据（presort = True）来加速训练，但是这对于较大训练集来说会显着减慢训练速度。

基尼不纯度或是信息熵

通常，算法使用 Gini 不纯度来进行检测，但是你也可以通过将标准超参数设置为”entropy”来使用熵不纯度进行检测。这里熵的概念是源于热力学中分子混乱程度的概念，当分子井然有序的时候，熵值接近于 0

那么我们到底应该使用 Gini 指数还是熵(ID3,C4.5决策树)呢？事实上大部分情况都没有多大的差别：他们会生成类似的决策树。基尼指数计算稍微快一点，所以这是一个很好的默认值。但是，也有的时候它们会产生不同的树，基尼指数会趋于在树的分支中将最多的类隔离出来，而熵指数趋向于产生略微平衡一些的决策树模型。

正则化超参数

决策树几乎不对训练数据做任何假设（于此相反的是线性回归等模型，这类模型通常会假设数据是符合线性关系的）。

如果不添加约束，树结构模型通常将根据训练数据调整自己，使自身能够很好的拟合数据，而这种情况下大多数会导致模型过拟合。

这一类的模型通常会被称为非参数模型，这不是因为它没有任何参数（通常也有很多），而是因为在训练之前没有确定参数的具体数量，所以模型结构可以根据数据的特性自由生长。

DecisionTreeClassifier类还有一些其他的参数用于限制树模型的形状:

min_samples_split（节点在被分裂之前必须具有的最小样本数），min_samples_leaf（叶节点必须具有的最小样本数），min_weight_fraction_leaf（和min_samples_leaf相同，但表示为加权总数的一小部分实例），max_leaf_nodes（叶节点的最大数量）和max_features（在每个节点被评估是否分裂的时候，具有的最大特征数量）。增加min_hyperparameters或者减少max_hyperparameters会使模型正则化
预剪枝与后剪枝

下图显示了对moons数据集进行训练生成的两个决策树模型，左侧的图形对应的决策树使用默认超参数生成（没有限制生长条件），右边的决策树模型设置为min_samples_leaf=4。很明显，左边的模型过拟合了，而右边的模型泛用性更好。

回归

决策树也能够执行回归任务，让我们使用 Scikit-Learn 的DecisionTreeRegressor类构建一个回归树，让我们用max_depth = 2在具有噪声的二次项数据集上进行训练。

1
2
3

from sklearn.tree import DecisionTreeRegressor
tree_reg = DecisionTreeRegressor(max_depth=2)
tree_reg.fit(X, y)

结果如图：

这棵树看起来非常类似于你之前建立的分类树，它的主要区别在于，它不是预测每个节点中的样本所属的分类，而是预测一个具体的数值。例如，假设您想对的新实例进行预测。从根开始遍历树，最终到达预测值等于 0.1106 的叶节点。该预测仅仅是与该叶节点相关的 110 个训练实例的平均目标值。而这个预测结果在对应的 110 个实例上的均方误差（MSE）等于 0.0151

下图的左侧显示的是模型的预测结果，如果你将max_depth=3设置为 3，模型就会如 6-5 图右侧显示的那样.注意每个区域的预测值总是该区域中实例的平均目标值。算法以一种使大多数训练实例尽可能接近该预测值的方式分割每个区域。

译者注：图里面的红线就是训练实例的平均目标值，对应上图中的value

CART 算法的工作方式与之前处理分类模型基本一样，不同之处在于，现在不再以最小化不纯度的方式分割训练集，而是试图以最小化 MSE 的方式分割训练集

下面公式显示了成本函数，该算法试图最小化这个成本函数。

和处理分类任务时一样，决策树在处理回归问题的时候也容易过拟合。如果不添加任何正则化（默认的超参数），你就会得到图 6-6 左侧的预测结果，显然，过度拟合的程度非常严重。而当我们设置了min_samples_leaf = 10，相对就会产生一个更加合适的模型了，就如图 6-6 所示的那样

不稳定性

它很容易理解和解释，易于使用且功能丰富而强大。然而，它也有一些限制，首先，你可能已经注意到了，决策树很喜欢设定正交化的决策边界，（所有边界都是和某一个轴相垂直的），这使得它对训练数据集的旋转很敏感，例如图 6-7 显示了一个简单的线性可分数据集。在左图中，决策树可以轻易的将数据分隔开，但是在右图中，当我们把数据旋转了 45° 之后，决策树的边界看起来变的格外复杂。尽管两个决策树都完美的拟合了训练数据，右边模型的泛化能力很可能非常差。

解决这个难题的一种方式是使用 PCA 主成分分析，这样通常能使训练结果变得更好一些。

更加通俗的讲，决策时的主要问题是它对训练数据的微小变化非常敏感，举例来说，我们仅仅从鸢尾花训练数据中将最宽的 Iris-Versicolor 拿掉（花瓣长 4.8 厘米，宽 1.8 厘米），然后重新训练决策树模型，你可能就会得到图 6-8 中的模型。正如我们看到的那样，决策树有了非常大的变化（原来的如图 6-2），事实上，由于 Scikit-Learn 的训练算法是非常随机的，即使是相同的训练数据你也可能得到差别很大的模型（除非你设置了随机数种子）

练习题

在 100 万例训练集上训练（没有限制）的决策树的近似深度是多少？
节点的基尼指数比起它的父节点是更高还是更低？它是通常情况下更高/更低，还是永远更高/更低？
如果决策树过拟合了，减少最大深度是一个好的方法吗？
如果决策树对训练集欠拟合了，尝试缩放输入特征是否是一个好主意？
如果对包含 100 万个实例的数据集训练决策树模型需要一个小时，在包含 1000 万个实例的培训集上训练另一个决策树大概需要多少时间呢？
如果你的训练集包含 100,000 个实例，设置presort=True会加快训练的速度吗？

1、$log_210^6\approx20$

2、子节点的基尼指数通常比父节点更低，这是由CART训练树损失函数决定的，它只会每次减少子节点权重求和后的Gini指数。然而子节点的Gini并不会永远小于父节点，比如现在根节点包括{A B A A A}五个样本，其Gini为$1-\frac{1}{5}^2-\frac{4}{5}^2=0.32$，划分样本后为{A B}和 {A A A} ，Gini分别为0.5和0。权重求和后的Gini为$\frac{2}{5}\times0.5+\frac{3}{5}\times0=0.2$，其Gini指数划分后只会比原来的更小。

3、如果决策树过度拟合训练集，那么减少max_depth是个不错的选择。

4、决策树不关心训练数据是否缩放或居中，这是他们的优势之一。因此，如果决策树如果欠拟合，缩放输入功能只会浪费时间。

5、训练决策树的时间复杂度为$O(nmlog_2m)$，m为样本数，n为特征数。所以训练样本增加十倍，其训练时间会增加$K=(n\cdot 10m\cdot log_2(10m))=10\times log(10m)/log(m)$.如果m=10^6,则K$\approx$11.7

6、仅当数据集小于几千个实例时，预先训练训练集能加速训练。如果它包含100,000个实例，则设置presort = True将大大减慢训练速度。

Sklearn 与 TensorFlow 机器学习实用指南（四）：支持向量机

发表于 2018-07-18 | 分类于 Sklearn 与 TensorFlow 机器学习实用指南

字数统计: 10.6k | 阅读时长 ≈ 41

推导

支持向量机条件描述

样本分类

Y∈{+1， -1}是样本的标签，分别代表两个不同的类。这里我们需要用这些样本去训练学习一个线性分类器（超平面）：$f(x)=sgn(w^Tx + b)$，也就是$w^Tx + b$大于0的时候，输出+1，小于0的时候，输出-1。sgn()表示取符号。而$g(x) =w^Tx + b=0$就是我们要寻找的分类超平面

也就是，对于任何一个正样本$y_i$=+1，它都要处于超平面的一边，也就是要保证：$y= w^Tx + b>0$。对于任何一个负样本$y_i=-1$，它都要处于超平面的另一边，也就是要保证：$y = w^Tx + b<0$这两个约束，其实可以合并成同一个式子：$y_i (w^Tx_i+b)\geq0$

最大间隔

点到直线距离距离引理：设直线 L 的方程为Ax+By+C=0，点 P 的坐标为（Xo，Yo），则点 P 到直线 L 的距离为

$\frac{\vert Ax_0+By_0+c\vert}{\sqrt{A^2+B^2}}$

间隔:间隔表示距离划分超平面最近的样本到划分超平面距离的两倍,即($x_i$为离直线最近的样本i)

$\gamma=2\min_i \frac{1}{\Vert w \Vert}\vert w^Tx_i+b \vert$

约束条件

所以，线性支持向量机的目标是找到一组合适的参数(w, b), 在满足分割样本的条件下，使得上诉间隔最大

$\max_{w,b} \min_i \frac{2}{\Vert w \Vert}\vert w^Tx_i+b \vert \\ \\ \\s.t. \quad y_i(w^T+b)>0, i=1,2,...,m$

上面的优化问题十分复杂, 难以处理. 为了能在现实中应用, 我们希望能对其做一些简化, 使其变为可以求解的, 经典的凸二次规划 (QP) 问题

凸二次规划的优化问题是指目标函数是凸二次函数, 约束是线性约束的一类优化问题

支持向量机的缩放引理：若（w* ,b*）是上面优化问题的解，那么对任何的r>0,(rw*,rb*)仍是该问题的解。

由于对 (w, b) 的放缩不影响解, 为了简化优化问题,我们约束 (w, b) 使得最短距离固定为1.

$\min_i \quad \vert w^Tx_i+b \vert = 1$

所以最后我们的优化目标可以等价为下面优化目标：

$\min_{w,b} \quad \frac{1}{2}w^Tw \\ \\ \\s.t. \quad y_i(w^T+b)\geq1, i=1,2,...,m$

这是一个凸二次规划问题，除了用解决QP问题的常规方法之外，还可以应用拉格朗日对偶性，通过求解对偶问题得到最优解，这就是线性可分条件下支持向量机的对偶算法，这样做的优点在于：一是对偶问题往往更容易求解；二者可以自然的引入核函数，进而推广到非线性分类问题。

拉格朗日函数与对偶形式

拉格朗日函数

对于优化问题

$\min_{u} \quad f(u) \\ \\ \\s.t. \quad g_i(u)\leq 0, i=1,2,...,m \\ \\ \\ \qquad \ \ \ h_j(u)=0, j=1,2,...,n \qquad（6）$

定义其拉格朗日函数为

$L(u,\alpha,\beta):=f(u)+\sum_{i=1}^m\alpha_ig_i(u)+\sum_{j=1}^n\beta_jh_j(u)\\ \\ \\s.t.\quad \alpha_i\geq0$

公式6的优化问题等价于

$\min \limits_{u} \max \limits_{\alpha,\beta} \quad L(u,\alpha,\beta) \\ \\ \\ s.t. \alpha_i\geq0, \ \ i=1,2,...,m \qquad(8)$

其证明过程如下

其中，当$g_i$不满足约束时，即$g_i(u)>0$,我们可以取$\alpha_i=\infty $;当$h_j$不满足约束时, 即$h_j(u)\neq 0$, 我们可以取$\beta_j=sign(h_j(u))\infty$，使得$\beta_jh_j(u)=\infty$. 当u满足约束时，由于$\alpha_i\geq 0,g_i(u)\leq 0$,则$\alpha_ig_i(u) \leq 0$.因此$\alpha_ig_i(u)$的最大值为0.

KKT条件和对偶形式

公式8描述的优化问题在最优值必须满足如下条件(KKT完整条件放后面)

主问题可行：$g_i(u)\leq 0,h_i(u)=0$;
对偶问题可行：$\alpha_i \geq 0$;
互补松弛：$\alpha_i g_i(u)=0$.

主问题可行和对偶问题可行是公式6和公式8描述的优化问题的约束项。松弛互补是在主问题和对偶问题都可行的条件下的最大值

在对偶问题里，公式6描述的优化问题，其等价形式公式8称为主问题，其对偶问题为(最大最小条件互换)

$\max \limits_{\alpha,\beta} \min \limits_{u} \quad L(u,\alpha,\beta) \\ \\ \\ s.t. \alpha_i\geq0, \ \ i=1,2,...,m$

引理一：对偶问题是主 (primal) 问题的下界（最大值里取最小值肯定会大于最小值里取最大值）

$\max \limits_{\alpha,\beta} \min \limits_{u} \quad L(u,\alpha,\beta)\leq \min \limits_{u} \max \limits_{\alpha,\beta} \quad L(u,\alpha,\beta)$

引理二(Slater条件)：当主问题为凸优化问题，即$f$和$g_i$ 为凸函数，$h_j$为仿射函数(类似于kx+b的一类函数)，且可行域中至少有一点使不等式约束严格成立时，对偶问题等价于原问题。(支持向量机的优化函数和约束函数都为凸函数)

线性支持向量机

线性支持向量机的拉格朗日函数为

$L(w,b,\alpha):=\frac{1}{2}w^tw+\sum_{i=1}^m\alpha_i(1-y_i(w^Tx_i+b))$

其对偶问题为（并非主问题）

$\max \limits_{\alpha}\min\limits_{w,b}\frac{1}{2}w^Tw+\sum_{i=1}^m\alpha_i(1-y_i(w^Tx_i+b)) \\ \\ \\ s.t.\quad \alpha_i\geq0,\quad i=1,2,...,m. \qquad(12)$

因为其内层的(w,b)属于无约束优化问题，我们可以通过令偏导等于0的方法得到(w,b)的最优值

$\frac{\partial L}{\partial w}=0\Rightarrow w=\sum_{i=1}^m{a_iy_ix_i}$ $\frac{\partial L}{\partial b}=0\Rightarrow \sum_{i=1}^m{a_iy_i=0}$

将上面公式带入(12)消去(w,b)即可得(变为min是因为化简过程有个负号，转为求最小)

$L\left( w,b,a \right) =\frac{1}{2}\sum_{i=1}^m{\sum_{j=1}^m{a_ia_jy_iy_j\left( x_i·x_j \right) -\sum_{i=1}^m{a_iy_i\left( \left( \sum_{j=1}^m{a_jy_jx_j} \right) ·x_i+b \right) +\sum_{i=1}^m{a_i}}}} \\ \\ \\ =-\frac{1}{2}\sum_{i=1}^m{\sum_{j=1}^m{a_ia_jy_iy_j\left( x_i·x_j \right) +\sum_{i=1}^m{a_i}}}$

故最终化简得到线性支持向量机对偶型

$\min \limits_{\alpha} \frac{1}{2}\sum_{i=1}^m\sum_{j=1}^m\alpha_i\alpha_jy_iy_jx_i^Tx_j-\sum_{i=1}^m\alpha_i \\ \\ \\ s.t.\sum_{i=1}^m\alpha_iy_i=0, \\ \\ \\ \alpha_i\geq 0,\quad i=1,2,...,m.$

上诉为线性支持向量机的最终化简对偶问题，其优化模型还要满足以下的KKT条件

主问题可行：$1-y_i(w^Tx_i+b)\leq0$
对偶问题可行：$\alpha_i\geq0$
互补松弛：$\alpha_i(1-y_i(w^Tx_i+b))=0$

根据KKT条件，我们可以得到两个非常重要的结论（这么多证明得出两个结论，不容易啊）

结论一：落在间隔边界上的支持向量，其对应的样本的对偶变量$\alpha_i>0$

证明：由KKT可知，根据$\alpha_i(1-y_i(w^Tx_i+b))=0$ 当;$\alpha_i>0$时，有$1-y_i(w^Tx_i+b)=0$ ，即$y_i(w^Tx_i+b)=1$

结论二：支持向量机的参数(w,b)仅仅由支持向量决定，与其他样本无关

证明：由于对偶变量$\alpha_i>0$对应的样本是支持向量，

$w=\sum_{i=1}^m\alpha_iy_ix_i=\sum_{i:\alpha_i=0}^m0\cdot y_ix_i+\sum_{i:\alpha_i>0}^m\alpha_iy_ix_i=\sum_{i\in SV}\alpha_iy_ix_i$

其中SV代表所有支持向量的集合.b可以由互补松弛算出，对于某一支持向量$x_s$ 及其标记$y_s$，由于$y_s(w^Tx_s+b)=1$,

则有

$b=y_s-w^Tx_s=y_s-\sum_{i\in SV}\alpha_iy_ix_i^Tx_s$

根据对偶公式(16)先计算最优解$\alpha_1,\alpha_2,…,\alpha_m$，然后可以得到w和b，这样我们就可以写出分类超平面$w^\ast\cdot x+b^\ast=0$ 和分类决策函数$f(x)=sign(w^\ast·x+b^\ast)$

在实践中，为了得到对b更稳健的估计,通常使用对所有支持向量求解得到b的平均值

软间隔支持向量机

对每个样本点引进一个松弛变量$\xi \geqslant 0 $，使函数间隔加上松弛变量大于等于1.这样，约束条件变成

$y_i\left( w·x_i+b \right) \geqslant 1-\xi _i$

当然，如果我们允许$\xi \geqslant 0 $任意大的话，那任意的超平面都是符合条件的了。这里我们引入松弛变量$\xi$，松弛变量值越大，当样本违背约束的程度越大，当松弛变量为0时，样本分类正确。所以，我们在原来的目标函数后面加上一项，使得这些 $\xi \geqslant 0 $的总和也要最小，优化最大间隔和违背约束程度越小，目标函数由原来的$\frac{1}{2}||w||^2$变成

$\frac{1}{2}||w||^2+C\sum_{i=1}^N{\xi _i}$

这里，$C>0$称为惩罚参数，一般事先由应用问题决定，用于权衡优化间隔和少量样本违背大间隔约束这两个目标，C越大时(看重这一项)对误分类的惩罚增大，，有更多的样本满足大间隔约束。C值小时对误分类的惩罚减小，允许有一些样本不满足大间隔约束.。最小化目标函数包含两层含义：使$\frac{1}{2}||w||^2$尽量最小间隔尽量大，同时使误分类点的个数尽量小，C是调和二者的系数。

则有软间隔支持向量机基本型：软间隔支持向量机旨在找到一组合适的参数 (w, b), 使得

$\underset{w,b,\xi}{\min}\ \frac{1}{2}||w||^2+C\sum_{i=1}^m{\xi _i} \\ \\ \\ s.t.\ \ y_i\left( w^T·x_i+b \right) \geqslant 1-\xi _i\ ,\ i=1,2,···,m \\ \\ \\ \xi _i\geqslant 0,\ i=1,2,···\mathrm{，}m \qquad(21)$

可证明w的解是唯一的，但b的解不唯一，b的解存在于一个区间。

用之前的方法将限制加入到目标函数中，得到如下原始最优化问题的拉格朗日函数(不等式为大于转为负号形式)：

$L\left( w,b,\xi ,\alpha,\beta \right) =\frac{1}{2}||w||^2+C\sum_{i=1}^m{\xi _i-\sum_{i=1}^m{\alpha_i\left( y_i\left( w^T·x_i+b \right) -1+\xi _i \right) -\sum_{i=1}^m{\beta_i\xi _i}}} \\ \\ \\ s.t.\alpha\geq0 ,\beta\geq0.$

其对偶问题为

$\max_\limits{\alpha,\beta}\min_\limits{w,b,\xi} L\left( \alpha,\beta,w,b,\xi \right) \qquad(23)$

首先求拉格朗日函数针对$w,b,\xi$求极小

$\frac{\partial L}{\partial w}=0\Rightarrow w=\sum_{i=1}^m{a_iy_ix_i} \\ \\ \\ \frac{\partial L}{\partial b}=0\Rightarrow \sum_{i=1}^m{\alpha_iy_i=0} \\ \\ \\ \frac{\partial L}{\partial \xi _i}=0\Rightarrow C-\alpha_i-\beta_i=0，i=1,2,3···,m \qquad(24)$

其中，有$\alpha\geq0,\beta\geq0$。因为$\beta_i=C-\alpha_i\geq0$，不失一般性，我们可以约束$0\leq\alpha_i\leq C$，从而去掉变量$\beta_i$.将这些求解变量带入拉格朗日对偶型即公式(23)，得到和原来一样的目标函数，唯一的区别就是现在拉格朗日乘子$\alpha$多了一个上限C。

最终化简的软间隔支持向量机对偶型问题(软间隔支持向量机的对偶问题等价于找到一组合适的$\alpha$使得)

$\underset{\alpha}{\min}\ \ \frac{1}{2}\sum_{i=1}^m{\sum_{j=1}^m{\alpha_i\alpha_jy_iy_j⟨x_i·x_j ⟩-\sum_{i=1}^m{\alpha_i}}} \\ \\ \\ s.t.\ \sum_{i=1}^m{a_iy_i=0} \\ \\ \\ 0\le a_i\le C\ ,\ \ i=1,2,···,m$

上诉为最终化简的软间隔支持向量机对偶型问题，其优化模型还要满足以下的KKT条件

主问题可行：$1-\xi_i-y_i(w^Tx_i+b)\leq0$,$-\xi\leq0$;
对偶问题可行：$\alpha_i\geq0$, $\beta_i\geq0$;
互补松弛：$\alpha_i(1-\xi_i-y_i(w^Tx_i+b))=0$,$\beta_i\xi_i=0$.

根据KKT条件，我们也可以得到两个结论：

结论一：软间隔支持向量机中, 支持向量落在最大间隔边界, 内部, 或被错误分类的样本

证明：由软间隔支持向量机的 KKT 条件可知，$\alpha_i(1-\xi_i-y_i(w^Tx_i+b))=0$且$\beta_i\xi_i=0$。当$\alpha_i>0$时，$1-\xi_i-y_i(w^Tx_i+b)=0$，进一步可以分为两种情况：

第一：$0<\alpha_i<’C$，此时$\beta_i$=C-$\alpha_i$>0(公式24的偏导约束)。因此$\xi_i$=0，即样本恰好落在最大间隔边界上。

第二：$\alpha_i=C$，此时$\beta_i=C-\alpha_i=0$. 若$\xi$ <1’则分类正确，该样本落在最大间隔内部(间隔边界与分离超平面之间)；若$\xi$=1，该样本在分隔超面上。若$\xi$>0样本位于分离超平面误分的一侧。

结论二：支持向量机的参数 (w, b) 仅由支持向量决定,与其他样本无关。

证明：和线性支持向量机证明方式相同.

求解w，b的方式也与线性支持向量机相同。注意：对任一适合条件0<a<C都可求得一个$b^∗$，即这些点正好是位于分隔边界上的点来求出b点的。但是由于原始问题对b的求解并不唯一，所以实际计算时可以取在所有符合条件的样本点上的平均值。

$w^*=\sum_{i=1}^N{a_{i}^{*}y_ix_i}$

选择a的一个分量的一个分量$a_i$适合约束条件适合约束条件0<$a_i$<C,计算

$b^*=y_j-\sum_{i=1}^N{a_{i}^{*}y_i⟨x_i·x_j ⟩}$

Hinge损失函数

线性支持向量机学习除了原始最优化问题，还有另外一种解释，就是最优化以下目标函数.下标”+”表示以下取正值的函数z=max(0，z)：

$\sum_i^{m}[1-y_i(w·x_i+b)]_++\lambda||w||^2$

目标函数的第一项是经验损失或经验风险，函数$L(y·(w·x+b))=[1-y(w·x+b)]_+$ ,称为合页损失函数（hinge loss function）.这就是说，当样本点$(x_i,y_i)$被分类正确且函数间隔（确信度）$y_i(w\cdot x_i+b)$大于1时，损失是0，否则损失是$1-y_i(w\cdot x_i+b)$。目标函数目标函数的第二项是系数的L2范数，是正则项。下面要证明上面的公式等价于线性支持向量机原始最优化问题，即(公式21)

证明：先令$[1-y_i(w·x_i+b)]_+=\xi_i$，则可以写成$\xi_i=\max(0, 1-y_i(w·x_i+b) )$ 。当样本满足约束分类正确时，$y_i(w\cdot x_i+b)>1$，有$1-y_i(w\cdot x_i+b)\leq0$，得$\xi_i=0$; 当当样本不满足约束时分类错误时，有$\xi_i=1-y_i(w\cdot x_i+b)$

故公式21的两个约束条件满足，其最优问题可以写作

$\underset{w,b}{min}\sum_{i=1}^m\xi_i+\lambda||w||^2$

若取$\lambda =\frac{1}{2C}$，则$\underset{w,b}{min} \frac{1}{C}(\frac{1}{2} ||w||^2+C\sum_{i=1}^N \xi_i)$ 与原始最优问题等价

图中还画出了0-1损失函数，可以认为它是一个二类分类问题的真正的损失函数，而合页损失函数是0-1损失函数的上界。由于0-1损失函数不是连续可导的，直接优化其构成的目标函数比较困难，可以认为线性支持向量机是优化由0-1损失函数的上界（合页损失函数）构成的目标函数。这时的上界损失函数又称为代理损失函数（surrogate function）。

图中虚线显示的是感知机的损失函数，相比之下，合页损失函数不仅要分类正确，而且确信度足够高时损失才是0，也就是说，合页损失函数对学习有更高的要求。

核函数

线性可分问题：既然在原始的特征空间$R^d$不是线性可分的, 支持向量机希望通过一个映射 $\phi$ :$R^d$ →
$R^{\widetilde{d}}$,使得数据在新的空间 $R^{\widetilde{d}}$ 是线性可分的.

令 $\phi$(x)代表将样本 x 映射到 $R^{\widetilde{d}}$中的特征向量,参数 w 的维数也要相应变为 $\widetilde{d}$维.

在上面的支持向量机对偶公式中，注意到，被映射到高维的特征向量总是以成对内积的形式存在，即 $\phi (x_i)^T\phi(x_j)$ ,如果先计算特征在$R^{\widetilde{d}}$ 空间的映射, 再计算内积, 复杂度是O($\widetilde{d}$) ,当特征被映射到非常高维的空间, 甚至是无穷维空间时, 这将会是沉重的存储和计算负担.

核技巧旨在将特征映射和内积这两步运算压缩为一步, 并且使复杂度由O($\widetilde{d}$)降为O($d$),即, 核技巧希望构造一个核函数 $\kappa (x_i,x_j) $ ,使得

$\kappa(x_i,x_j)=\phi (x_i)^T\phi(x_j)$

并且，$\kappa (x_i,x_j) $ 的计算复杂度是 O($d$) .

应用

软间隔分类

如果我们严格地规定所有的数据都不在“街道”上，都在正确地两边，称为硬间隔分类，硬间隔分类有两个问题，第一，只对线性可分的数据起作用，第二，对异常点敏感。

为了避免上述的问题，我们更倾向于使用更加软性的模型。目的在保持“街道”尽可能大和避免间隔违规（例如：数据点出现在“街道”中央或者甚至在错误的一边）之间找到一个良好的平衡。这就是软间隔分类。

在 Scikit-Learn 库的 SVM 类，你可以用C超参数（惩罚系数）来控制这种平衡：较小的C会导致更宽的“街道”，但更多的间隔违规。下图显示了在非线性可分隔的数据集上，两个软间隔SVM分类器的判定边界。左边图中，使用了较大的C值，导致更少的间隔违规，但是间隔较小。右边的图，使用了较小的C值，间隔变大了，但是许多数据点出现在了“街道”上。然而，第二个分类器似乎泛化地更好：事实上，在这个训练数据集上减少了预测错误，因为实际上大部分的间隔违规点出现在了判定边界正确的一侧

如果你的 SVM 模型过拟合，你可以尝试通过减小超参数C去调整(街道更宽)
SVM 特别适合复杂的分类，而中小型的数据集分类中很少用到。另外，SVM 对特征缩放比较敏感，要先做数据缩放处理。

以下的 Scikit-Learn 代码加载了内置的鸢尾花（Iris）数据集，缩放特征，并训练一个线性 SVM 模型（使用LinearSVC类，超参数C=1，hinge 损失函数）来检测 Virginica 鸢尾花。

import numpy as np
from sklearn import datasets
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import LinearSVC

iris = datasets.load_iris()
X = iris["data"][:, (2, 3)] # petal length, petal width
y = (iris["target"] == 2).astype(np.float64) # Iris-Virginica

svm_clf = Pipeline((
        ("scaler", StandardScaler()),
        ("linear_svc", LinearSVC(C=1, loss="hinge")),
    ))

svm_clf.fit(X_scaled, y)

Then, as usual, you can use the model to make predictions:

>>> svm_clf.predict([[5.5, 1.7]])
array([ 1.])

不同于 Logistic 回归分类器，SVM 分类器不会输出每个类别的概率。
SVM的分类器为SVC类，回归则为SVR类。

作为一种选择，你可以在 SVC 类，使用SVC(kernel="linear", C=1)，但是它比较慢，尤其在较大的训练集上，所以一般不被推荐。另一个选择是使用SGDClassifier类，即SGDClassifier(loss="hinge", alpha=1/(m*C))。它应用了随机梯度下降（SGD 见第四章）来训练一个线性 SVM 分类器。尽管它不会和LinearSVC一样快速收敛，但是对于处理那些不适合放在内存的大数据集是非常有用的，或者处理在线分类任务同样有用。

LinearSVC要使偏置项规范化，首先你应该集中训练集减去它的平均数，完成数据中心化。如果你使用了StandardScaler，那么它会自动处理。此外，确保你设置loss参数为hinge，因为它不是默认值。最后，为了得到更好的效果，你需要将dual参数设置为False，除非特征数比样本量多。

常用损失函数

铰链损失（Hinge Loss）：主要用于支持向量机（SVM）中，二分类SVM等于Hinge损失+ L2正则化。
互熵损失（Cross Entropy Loss，Softmax Loss ）：用于Logistic 回归与Softmax 分类中；
平方损失（Square Loss）：主要是最小二乘法（OLS）中；
指数损失（Exponential Loss） ：主要用于Adaboost 集成学习算法中；

Hinge loss 的叫法来源于其损失函数的图形，为一个折线，通用的函数表达式为：

$L(m_i) = max(0,1-m_i(w)) =\max(0, 1-y\tilde{y}), 其中y=\pm 1$

表示如果被正确分类，损失是0，否则损失就是 $1-m_i(w)$

在机器学习中，Hing 可以用来解 间距最大化 的问题，最有代表性的就是SVM 问题，最初的SVM 优化函数如下：

$\underset{w,\zeta}{argmin} \frac{1}{2}||w||^2+ C\sum_i \zeta_i \\ st.\quad \forall y_iw^Tx_i \geq 1- \zeta_i \\ \zeta_i \geq 0$

将约束项进行变形，则为：

$\zeta_i \geq 1-y_iw^Tx_i$

则损失函数可以进一步写为：

$\begin{equation}\begin{split}J(w)&=\frac{1}{2}||w||^2 + C\sum_i max(0,1-y_iw^Tx_i) \\ &= \frac{1}{2}||w||^2 + C\sum_i max(0,1-m_i(w)) \\ &= \frac{1}{2}||w||^2 + C\sum_i L_{Hinge}(m_i) \end{split}\end{equation}$

因此， SVM 的损失函数可以看作是 L2-norm 和 Hinge loss 之和。

非线性支持向量机分类

创建多项式特征

from sklearn.datasets import make_moons
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures

polynomial_svm_clf = Pipeline((
        ("poly_features", PolynomialFeatures(degree=3)),
        ("scaler", StandardScaler()),
        ("svm_clf", LinearSVC(C=10, loss="hinge"))
    ))

polynomial_svm_clf.fit(X, y)

多项式核

添加多项式特征很容易实现，不仅仅在 SVM，在各种机器学习算法都有不错的表现，但是低次数的多项式不能处理非常复杂的数据集，而高次数的多项式却产生了大量的特征，会使模型变得慢.

幸运的是，当你使用 SVM 时，你可以运用一个被称为“核技巧”（kernel trick）的神奇数学技巧。它可以取得就像你添加了许多多项式，甚至有高次数的多项式，一样好的结果。所以不会大量特征导致的组合爆炸，因为你并没有增加任何特征。超参数coef0控制了高阶多项式与低阶多项式对模型的影响。

from sklearn.svm import SVC
poly_kernel_svm_clf = Pipeline((
        ("scaler", StandardScaler()),
        ("svm_clf", SVC(kernel="poly", degree=3, coef0=1, C=5))
    ))
poly_kernel_svm_clf.fit(X, y)

增加相似特征

另一种解决非线性问题的方法是使用相似函数（similarity funtion）计算每个样本与特定地标（landmark）的相似度。例如，让我们来看看前面讨论过的一维数据集，并在x1=-2和x1=1之间增加两个地标。接下来，我们定义一个相似函数，即高斯径向基函数（Gaussian Radial Basis Function，RBF），设置γ = 0.3。

$RBF:\quad \phi_{\gamma}(x, \ell) = exp(-\gamma \|x - \ell \|^2)$

它是个从 0 到 1 的钟型函数，值为 0 的离地标很远，值为 1 的在地标上。现在我们准备计算新特征。例如，我们看一下样本x0=-1：它距离第一个地标距离(x1=-2)是 1，距离第二个地标(x1=1)是 2。因此它的新特征为x2=exp((-0.3 × 1)^2)≈0.74和x3=exp((-0.3 × 2)^2)≈0.30。右边的图显示了特征转换后的数据集（删除了原始特征），正如你看到的，它现在是线性可分了。

你可能想知道如何选择地标。最简单的方法是在数据集中的每一个样本的位置创建地标。这将产生更多的维度从而增加了转换后数据集是线性可分的可能性。但缺点是，m个样本，n个特征的训练集被转换成了m个实例，m个特征的训练集（假设你删除了原始特征）。这样一来，如果你的训练集非常大，你最终会得到同样大的特征

高斯 RBF 核

就像多项式特征法一样，相似特征法对各种机器学习算法同样也有不错的表现。但是在所有额外特征上的计算成本可能很高，特别是在大规模的训练集上。然而，“核” 技巧再一次显现了它在 SVM 上的神奇之处：高斯核让你可以获得同样好的结果成为可能，就像你在相似特征法添加了许多相似特征一样，但事实上，你并不需要在RBF添加它们。我们使用 SVC 类的高斯 RBF 核来检验一下。

rbf_kernel_svm_clf = Pipeline((
        ("scaler", StandardScaler()),
        ("svm_clf", SVC(kernel="rbf", gamma=5, C=0.001))
    ))
rbf_kernel_svm_clf.fit(X, y)

下图显示了用不同的超参数gamma (γ)和C训练的模型。增大γ使钟型曲线更窄，导致每个样本的影响范围变得更小：即判定边界最终变得更不规则，在单个样本周围环绕。相反的，较小的γ值使钟型曲线更宽，样本有更大的影响范围，判定边界最终则更加平滑。所以γ是可调整的超参数：如果你的模型过拟合，你应该减小γ值，若欠拟合，则增大γ（与超参数C相似）。

还有其他的核函数，但很少使用。例如，一些核函数是专门用于特定的数据结构。在对文本文档或者 DNA 序列进行分类时，有时会使用字符串核（String kernels）（例如，使用 SSK 核（string subsequence kernel）或者基于编辑距离（Levenshtein distance）的核函数）

这么多可供选择的核函数，你如何决定使用哪一个？一般来说，你应该先尝试线性核函数（记住LinearSVC比SVC(kernel=”linear”)要快得多），尤其是当训练集很大或者有大量的特征的情况下。如果训练集不太大，你也可以尝试高斯径向基核（Gaussian RBF Kernel），它在大多数情况下都很有效。如果你有空闲的时间和计算能力，你还可以使用交叉验证和网格搜索来试验其他的核函数，特别是有专门用于你的训练集数据结构的核函数

计算复杂性

LinearSVC类基于liblinear库，它实现了线性 SVM 的优化算法。它并不支持核技巧，但是它样本和特征的数量几乎是线性的：训练时间复杂度大约为O(m × n)。

如果你要非常高的精度，这个算法需要花费更多时间。这是由容差值超参数ϵ（在 Scikit-learn 称为tol）控制的。大多数分类任务中，使用默认容差值的效果是已经可以满足一般要求。

SVC 类基于libsvm库，它实现了支持核技巧的算法。训练时间复杂度通常介于O(m^2 × n)和O(m^3 × n)之间。不幸的是，这意味着当训练样本变大时，它将变得极其慢（例如，成千上万个样本）。这个算法对于复杂但小型或中等数量的数据集表现是完美的。然而，它能对特征数量很好的缩放，尤其对稀疏特征来说（sparse features）（即每个样本都有一些非零特征）。在这个情况下，算法对每个样本的非零特征的平均数量进行大概的缩放。下表对 Scikit-learn 的 SVM 分类模型进行比较。

SVM 回归

正如我们之前提到的，SVM 算法应用广泛：不仅仅支持线性和非线性的分类任务，还支持线性和非线性的回归任务。技巧在于逆转我们的目标：限制间隔违规的情况下，不是试图在两个类别之间找到尽可能大的“街道”（即间隔）。SVM 回归任务是限制间隔违规情况下，尽量放置更多的样本在“街道”上。“街道”的宽度由超参数ϵ控制。下图显示了在一些随机生成的线性数据上，两个线性 SVM 回归模型的训练情况。一个有较大的间隔（ϵ=1.5），另一个间隔较小（ϵ=0.5）

添加更多的数据样本在间隔之内并不会影响模型的预测，因此，这个模型认为是不敏感的（ϵ-insensitive）。

你可以使用 Scikit-Learn 的LinearSVR类去实现线性 SVM 回归。

1
2
3

from sklearn.svm import LinearSVR
svm_reg = LinearSVR(epsilon=1.5)
svm_reg.fit(X, y)

处理非线性回归任务，你可以使用核化的 SVM 模型。比如，图显示了在随机二次方的训练集，使用二次方多项式核函数的 SVM 回归。左图是较小的正则化（即更大的C值），右图则是更大的正则化（即小的C值）。

下面的代码使用了 Scikit-Learn 的SVR类（支持核技巧）。在回归任务上，SVR类和SVC类是一样的，并且LinearSVR是和LinearSVC等价。LinearSVR类和训练集的大小成线性（就像LinearSVC类），当训练集变大，SVR会变的很慢（就像SVC类）

from sklearn.svm import SVR

svm_poly_reg = SVR(kernel="poly", degree=2, C=100, epsilon=0.1)
svm_poly_reg.fit(X, y)

SVM 也可以用来做异常值检测，详情见 Scikit-Learn 文档

补充

拉格朗日子乘与KKT

等式约束问题

可以得到一个重要的结论：▽f(x)一定与▽h(x)平行，故其关系可以写成▽f(x) =λ ▽h(x).

扩展到高维度等式：f(x)为目标优化函数，hi(x)为约束等式，其表达式恒等于0。F(x)为等价的拉格朗日子乘函数。

计算 F 对x与$\lambda$ 的偏导数并令其为零，可得最优解的必要条件：

其中第一式为定常方程式(stationary equation)，第二式为约束条件。也就是之前已经证明过得结果

故等式约束要满足的条件为

不等式约束问题

这两种情况的最佳解具有两种不同的必要条件

内部解：h(x)<0，不满足h(x)的等式约束，故不在边界上，在可行域的内部。在约束条件无效的情形下， h(x)不起作用，约束优化问题退化为无约束优化问题，因此驻点(最优解)满足▽f(x)=0且$\lambda$=0 ($\lambda$为0才能消去F(x)的约束项)。
边界解：在约束条件有效的情形下，约束不等式变成等式约束，有h(x)=0。和等式约束做相同的处理。另外，存在$\lambda$ 使得$\bigtriangledown f=-\lambda \bigtriangledown g$，但这里$\lambda$ 的正负号是有其意义的。因为我们希望最小化 f，梯度$\bigtriangledown f$ (函数 f在点x的最陡上升方向)应该指向可行域的内部(因为你的最优解最小值是在边界取得的)，但 $\bigtriangledown g$指向可行域的的外部(即 g(x)>0的区域，因为你的约束是小于等于0)，因此$\lambda \geq0$，称为对偶可行性(dual feasibility)。

所以，不论是内部解或边界解，$\lambda h(x)=0$恒成立(不同情况总有一个为0)，称为互补松弛性(complementary slackness)。整合上述两种情况，最佳解的必要条件包括拉格朗日函数常定方程式、主问题可行，对偶可行，和互补松弛。

对偶性

一个优化问题，通过求出它的 dual problem ，在只有 weak duality 成立的情况下，我们至少可以得到原始问题的一个下界。而如果 strong duality 成立，则可以直接求解 dual problem 来解决原始问题，就如同经典的 SVM 的求解过程一样。有可能 dual problem 比 primal problem 更容易求解，或者 dual problem 有一些优良的结构（例如 SVM 中通过 dual problem 我们可以将问题表示成数据的内积形式从而使得 kernel trick 的应用成为可能）。此外，还有一些情况会同时求解 dual 和 primal problem ，比如在迭代求解的过程中，通过判断 duality gap 的大小，可以得出一个有效的迭代停止条件.

支持向量机的其他变体

Prob SVM.

Prob SVM. 对数几率回归可以估计出样本属于正类的概率, 而支持向量机只能判断样本属于正类或负类,无法得到概率.Prob SVM 先训练一个支持向量机, 得到参数(w,b)。再令$s_i:=y_iw^Tx_i+b$，将${(s_1,y_1),(s_2,y_2),…,(s_m,y_m)}$当做新的训练数据训练一个对数几率回归模型，得到参数$(\theta_1,\theta_0)$.因此，ProbSVM 的假设函数为

$h(x):=sigm(\theta_1(w^Tx_i+b)+\theta_0)$

对数几率回归模型可以认为是对训练得到的支持向量机的微调, 包括尺度 (对应$\theta_1)$和平移(对应$\theta_0$).

多分类支持向量机

支持向量机也可以扩展到多分类问题中. 对于 K 分类问题, 多分类支持向量机有 K组参数 ${(w_1,b_1),(w_2,b_2),…,(w_K,b_K)}$并希望模型对于属于正确标记的结果以 1 的间隔高于其他类的结果, 形式化如下

$\min \limits_{W,b} \frac{1}{m}\sum_{i=1}^m\sum_{k=1}^Kmax(0,(w_{yi}^Tx_i+b_{yi})-(w_{k}^Tx_i+b_{k}+1)+\frac{2}{\lambda}\sum_{k=1}^Kw_k^Tw_k$

支持向量回归 (SVR)

支持向量回归 (SVR). 经典回归模型的损失函数度量了模型的预测 $h(x_i)$与$y_i$的差别, 支持向量回归能够容忍$h(x_i)$与$y_i$之间小于$\epsilon$ 的偏差.令$s:=y-w^Tx+b$，我们定义$\epsilon$不敏感损失为

支持向量机和LR的异同

SVM与LR的相同点

LR和SVM都是分类算法，都是监督学习算法。
如果不考虑核函数，LR和SVM都是线性分类算法，也就是说他们的分类决策面都是线性的。LR也是可以用核函数的，至于为什么通常在SVM中运用核函数而不在LR中运用，后面讲到他们之间区别的时候会重点分析。总之，原始的LR和SVM都是线性分类器，这也是为什么通常没人问你决策树和LR什么区别，决策树和SVM什么区别，你说一个非线性分类器和一个线性分类器有什么区别？
LR和SVM都是判别模型。判别模型会生成一个表示P(Y|X)的判别函数（或预测模型），而生成模型先计算联合概率p(Y,X)然后通过贝叶斯公式转化为条件概率。简单来说，在计算判别模型时，不会计算联合概率，而在计算生成模型时，必须先计算联合概率。或者这样理解：生成算法尝试去找到底这个数据是怎么生成的（产生的），然后再对一个信号进行分类。基于你的生成假设，那么那个类别最有可能产生这个信号，这个信号就属于那个类别。判别模型不关心数据是怎么生成的，它只关心信号之间的差别，然后用差别来简单对给定的一个信号进行分类。常见的判别模型有：KNN、SVM、LR，常见的生成模型有：朴素贝叶斯，隐马尔可夫模型。当然，这也是为什么很少有人问你朴素贝叶斯和LR以及朴素贝叶斯和SVM有什么区别。
LR和SVM在学术界和工业界都广为人知并且应用广泛。

SVM与LR的不同点

1.损失函数

SVM的处理方法是只考虑support vectors，也就是和分类最相关的少数点，去学习分类器。而逻辑回归通过非线性映射，大大减小了离分类平面较远的点的权重，相对提升了与分类最相关的数据点的权重,两者的根本目的都是一样的。即支持向量机只考虑局部的边界线附近的点，而逻辑回归考虑全局（远离的点对边界线的确定也起作用）。

影响SVM决策面的样本点只有少数的支持向量，当在支持向量外添加或减少任何样本点对分类决策面没有任何影响；而在LR中，每个样本点都会影响决策面的结果

2.核技巧

在解决非线性问题时，支持向量机采用核函数的机制，而LR通常不采用核函数的方法。

这个问题理解起来非常简单。分类模型的结果就是计算决策面，模型训练的过程就是决策面的计算过程。通过上面的第二点不同点可以了解，在计算决策面时，SVM转化为对偶问题后，只有少数几个代表支持向量的样本参与了计算，也就是只有少数几个样本需要参与核计算（即kernal machine解的系数是稀疏的），这个在进行复杂核函数计算时优势很明显，能够大大简化模型和计算量。。然而，LR算法里，每个样本点都必须参与决策面的计算过程，也就是说，假设我们在LR里也运用核函数的原理，那么每个样本点都必须参与核计算，这带来的计算复杂度是相当高的。所以，在具体应用时，LR很少运用核函数机制。

3.正则项

根据需要，两个方法都可以增加不同的正则化项，如l1,l2等等。所以在很多实验中，两种算法的结果是很接近的。但是逻辑回归相对来说模型更简单，好理解，实现起来，特别是大规模线性分类时比较方便。而SVM的理解和优化相对来说复杂一些。但是SVM的理论基础更加牢固，有一套结构化风险最小化的理论基础，虽然一般使用的人不太会去关注。

4.异常值

两者对异常的敏感度也不一样。同样的线性分类情况下，如果异常点较多的话，无法剔除，首先LR，LR中每个样本都是有贡献的，最大似然后会自动压制异常的贡献，SVM+软间隔对异常还是比较敏感，因为其训练只需要支持向量，有效样本本来就不高，一旦被干扰，预测结果难以预料。

5.normalization

两个模型对数据和参数的敏感程度不同，Linear SVM比较依赖penalty的系数和数据表达空间的测度，而（带正则项的）LR比较依赖对参数做L1 regularization的系数。但是由于他们或多或少都是线性分类器，所以实际上对低维度数据overfitting的能力都比较有限，相比之下对高维度数据，LR的表现会更加稳定，为什么呢？

因为Linear SVM在计算margin有多“宽”的时候是依赖数据表达上的距离测度的，换句话说如果这个测度不好（badly scaled，这种情况在高维数据尤为显著），所求得的所谓Large margin就没有意义了，这个问题即使换用kernel trick（比如用Gaussian kernel）也无法完全避免。所以使用Linear SVM之前一般都需要先对数据做normalization，而求解LR（without regularization）时则不需要或者结果不敏感。

练习题

支持向量机背后的基本思想是什么
什么是支持向量
当使用 SVM 时，为什么标准化输入很重要？
分类一个样本时，SVM 分类器能够输出一个置信值吗？概率呢？
在一个有数百万训练样本和数百特征的训练集上，你是否应该使用 SVM 原始形式或对偶形式来训练一个模型？
假设你用 RBF 核来训练一个 SVM 分类器，如果对训练集欠拟合：你应该增大或者减小γ吗？调整参数C呢？

1、支持向量机背后的基本目标是在训练实例中分隔两个类的决策边界之间具有最大可能的间隔。当执行软间隔分类时，SVM在完全分离两个类和具有尽可能宽的街道之间搜索折衷。另一个关键思想是在训练非线性数据集时使用核技巧。

2、在训练SVM之后，支持向量是位于“街道”上的任何实例，包括其边界。决策边界完全由支持向量决定。任何不是支持向量的实例（即街道外）都没有任何影响; 你可以删除它们，添加更多实例或移动它们，只要它们离开街道它们就不会影响决策边界。计算预测仅涉及支持向量，而不是整个训练集。

3、SVM尝试适应类之间最大可能的“街道”，因此如果训练集未缩放，SVM将倾向于忽略数值小的特征。

4、SVM分类器可以输出测试实例与决策边界之间的距离，您可以将其用作置信度分数。但是，这个分数不能直接转换为类概率的估计。如果在Scikit-Learn中创建SVM时设置probability = True，则在训练之后，它将使用SVM分数的Logistic回归校准概率（通过对训练数据进行额外的5折交叉验证进行训练）。这会将predict_proba（）和predict_log_proba（）方法添加到SVM。

5、此问题仅适用于线性SVM，因为进行核化kernelized需要用到对偶形式(对偶问题就是为了可以进行核化)。 SVM问题的原始形式的计算复杂度与训练实例m的数量成比例，而对偶形式的计算复杂度与m2和m3之间的数量成比例。所以如果有数百万个实例，你肯定应该使用原始形式，因为对偶形式会太慢。

6、如果使用RBF内核训练的SVM分类器不适合训练集，则可能存在正则化强度过高。要减少它，您需要增加gamma或C（或两者）。