机器学习实战：朴素贝叶斯和Logistic回归

￡神魔★判官ぃ 2023-06-17 11:55 47阅读 0赞

# 机器学习实战 #

### 文章目录 ###

*  机器学习实战
 *   *  一、基于概率论的分类方法：朴素贝叶斯
     *   *  1、朴素贝叶斯 概述
         *  2、贝叶斯决策理论 & 条件概率
         *   *  贝叶斯决策理论
             *  条件概率
             *  使用条件概率来分类
         *  3、朴素贝叶斯 原理
         *   *  工作原理
             *  开发流程
             *  算法特点
         *  4、朴素贝叶斯 项目案例
         *   *  项目案例1: 屏蔽社区留言板的侮辱性言论
             *  项目案例2: 使用朴素贝叶斯过滤垃圾邮件
             *  项目案例3: 使用朴素贝叶斯分类器从个人广告中获取区域倾向
         *  5、朴素贝叶斯 小结
     *  二、Logistic回归
     *   *  1、Logistic 回归 概述
         *   *  须知概念
         *  2、Logistic 回归 原理
         *   *  工作原理
             *  开发流程
             *  算法特点
         *  3、Logistic 回归 项目案例
         *   *  项目案例1: 使用 Logistic 回归在简单数据集上的分类
             *  项目案例2: 从疝气病症预测病马的死亡率
         *  4、Logistic回归 小结
 *  资料来源

## 一、基于概率论的分类方法：朴素贝叶斯 ##

### 1、朴素贝叶斯 概述 ###

**贝叶斯分类**是一类分类算法的总称，这类算法均以贝叶斯决策理论为基础，故统称为贝叶斯分类。  
我们通过实例来讨论贝叶斯分类的中最简单的一种: 朴素贝叶斯分类。我们称之为“朴素”，是因为整个形式化过程只做原始、 简单的假设。

### 2、贝叶斯决策理论 & 条件概率 ###

#### 贝叶斯决策理论 ####

我们现在有一个数据集，它由两类数据组成，数据分布如下图所示：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70]  
我们现在用 p1(x,y) 表示数据点 (x,y) 属于类别1（图中用圆点表 示的类别）的概率，用 p2(x,y) 表示数据点 (x,y) 属于类别2（图中用三角形表示的类别）的概率， 那么对于一个新数据点 (x,y)，可以用下面的规则来判断它的类别：

*  如果 p1(x,y) > p2(x,y) ，那么类别为1
 *  如果 p2(x,y) > p1(x,y) ，那么类别为2

也就是说，我们会选择高概率对应的类别。这就是贝叶斯决策理论的核心思想，即选择具有最高概率的决策。

#### 条件概率 ####

如果你对 p(x,y|c1) 符号很熟悉，那么可以跳过本小节。

有一个装了 7 块石头的罐子，其中 3 块是白色的，4 块是黑色的。如果从罐子中随机取出一块石头，那么是白色石头的可能性是多少？由于取石头有 7 种可能，其中 3 种为白色，所以取出白色石头的概率为 3/7 。那么取到黑色石头的概率又是多少呢？很显然，是 4/7 。我们使用 P(white) 来表示取到白色石头的概率，其概率值可以通过白色石头数目除以总的石头数目来得到。  
![在这里插入图片描述][20191128212448540.png]  
如果这 7 块石头如下图所示，放在两个桶中，那么上述概率应该如何计算？  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 1]  
计算 P(white) 或者 P(black) ，如果事先我们知道石头所在桶的信息是会改变结果的。这就是所谓的条件概率（conditional probablity）。假定计算的是从 B 桶取到白色石头的概率，这个概率可以记作 P(white|bucketB) ，我们称之为“在已知石头出自 B 桶的条件下，取出白色石头的概率”。很容易得到，P(white|bucketA) 值为 2/4 ，P(white|bucketB) 的值为 1/3 。

条件概率的计算公式如下：

P(white|bucketB) = P(white and bucketB) / P(bucketB)

首先，我们用 B 桶中白色石头的个数除以两个桶中总的石头数，得到 P(white and bucketB) = 1/7 .其次，由于 B 桶中有 3 块石头，而总石头数为 7 ，于是 P(bucketB) 就等于 3/7 。于是又 P(white|bucketB) = P(white and bucketB) / P(bucketB) = (1/7) / (3/7) = 1/3 。

另外一种有效计算条件概率的方法称为贝叶斯准则。贝叶斯准则告诉我们如何交换条件概率中的条件与结果，即如果已知 P(x|c)，要求 P(c|x)，那么可以使用下面的计算方法：  
![在这里插入图片描述][20191128212615487.png]

#### 使用条件概率来分类 ####

上面我们提到贝叶斯决策理论要求计算两个概率 p1(x, y) 和 p2(x, y):

*  如果 p1(x, y) > p2(x, y), 那么属于类别 1;
 *  如果 p2(x, y) > p1(X, y), 那么属于类别 2.

这并不是贝叶斯决策理论的所有内容。使用 p1() 和 p2() 只是为了尽可能简化描述，而真正需要计算和比较的是 p(c1|x, y) 和 p(c2|x, y) .这些符号所代表的具体意义是: 给定某个由 x、y 表示的数据点，那么该数据点来自类别 c1 的概率是多少？数据点来自类别 c2 的概率又是多少？注意这些概率与概率 p(x, y|c1) 并不一样，不过可以使用贝叶斯准则来交换概率中条件与结果。具体地，应用贝叶斯准则得到:  
![在这里插入图片描述][20191128212720251.png]  
使用上面这些定义，可以定义贝叶斯分类准则为:

*  如果 P(c1|x, y) > P(c2|x, y), 那么属于类别 c1;
 *  如果 P(c2|x, y) > P(c1|x, y), 那么属于类别 c2.

在文档分类中，整个文档（如一封电子邮件）是实例，而电子邮件中的某些元素则构成特征。我们可以观察文档中出现的词，并把每个词作为一个特征，而每个词的出现或者不出现作为该特征的值，这样得到的特征数目就会跟词汇表中的词的数目一样多。

我们假设特征之间 **相互独立** 。所谓 **独立(independence)** 指的是统计意义上的独立，即一个特征或者单词出现的可能性与它和其他单词相邻没有关系，比如说，“我们”中的“我”和“们”出现的概率与这两个字相邻没有任何关系，虽然我们知道这个假设并不正确。这个假设正是朴素贝叶斯分类器中 朴素(naive) 一词的含义。朴素贝叶斯分类器中的另一个假设是，**每个特征同等重要**。其实这个假设也有问题。 如果要判断留言板的留言是否得当，那么可能不需要看 完所有的1000个单词，而只需要看10～20个特征就足以做出判断了。尽管上述假设存在一些小的 瑕疵，但朴素贝叶斯的实际效果却很好。

**Note**: 朴素贝叶斯分类器通常有两种实现方式: 一种基于伯努利模型实现，一种基于多项式模型实现。这里采用前一种实现方式。该实现方式中并不考虑词在文档中出现的次数，只考虑出不出现，因此在这个意义上相当于假设词是等权重的。

### 3、朴素贝叶斯 原理 ###

#### 工作原理 ####

提取所有文档中的词条并进行去重
    获取文档的所有类别
    计算每个类别中的文档数目
    对每篇训练文档: 
        对每个类别: 
            如果词条出现在文档中-->增加该词条的计数值（for循环或者矩阵相加）
            增加所有词条的计数值（此类别下词条总数）
    对每个类别: 
        对每个词条: 
            将该词条的数目除以总词条数目得到的条件概率（P(词条|类别)）
    返回该文档属于每个类别的条件概率（P(类别|文档的所有词条)）

#### 开发流程 ####

> 1.  收集数据: 可以使用任何方法。
> 2.  准备数据: 需要数值型或者布尔型数据。
> 3.  分析数据: 有大量特征时，绘制特征作用不大，此时使用直方图效果更好。
> 4.  训练算法: 计算不同的独立特征的条件概率。
> 5.  测试算法: 计算错误率。
> 6.  使用算法: 一个常见的朴素贝叶斯应用是文档分类。可以在任意的分类场景中使用朴素贝叶斯分类器，不一定非要是文本。

#### 算法特点 ####

>  *  优点：在数据较少的情况下仍然有效，可以处理多类别问题。
>  *  缺点：对于输入数据的准备方式较为敏感。
>  *  适用数据类型：标称型数据。

### 4、朴素贝叶斯 项目案例 ###

完整代码地址：[https://blog.csdn.net/qq\_45556599/article/details/103315806][https_blog.csdn.net_qq_45556599_article_details_103315806]

#### 项目案例1: 屏蔽社区留言板的侮辱性言论 ####

**项目概述**

构建一个快速过滤器来屏蔽在线社区留言板上的侮辱性言论。如果某条留言使用了负面或者侮辱性的语言，那么就将该留言标识为内容不当。对此问题建立两个类别: 侮辱类和非侮辱类，使用 1 和 0 分别表示。

**开发流程**

> 收集数据: 可以使用任何方法  
> 准备数据: 从文本中构建词向量  
> 分析数据: 检查词条确保解析的正确性  
> 训练算法: 从词向量计算概率  
> 测试算法: 根据现实情况修改分类器  
> 使用算法: 对社区留言板言论进行分类

*收集数据: 可以使用任何方法*

我们将把文本看成单词向量或者词条向量，也就是说将句子转换为向量。考虑出现在所有文 档中的所有单词，再决定将哪些词纳入词汇表或者说所要的词汇集合，然后必须要将每一篇文档 转换为词汇表上的向量。

本例是我们自己构造的词表:

def loadDataSet():
        """
        创建数据集
        :return: 单词列表postingList, 所属类别classVec
        """
        postingList = [['my', 'dog', 'has', 'flea', 'problems', 'help', 'please'], #[0,0,1,1,1......]
                       ['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

*准备数据: 从文本中构建词向量*

def createVocabList(dataSet):
        """
        获取所有单词的集合
        :param dataSet:数据集
        :return:所有单词的集合(即不含重复元素的单词列表)
        """
        vocabSet = set([])  # 创建空集合
        for document in dataSet:
            # 操作符 | 用于求两个集合的并集
            vocabSet = vocabSet | set(document)
        return list(vocabSet)
    
    
    def setOfWords2Vec(vocabList, inputSet):
        """
        遍历查看该单词是否出现，出现该单词则将该单词置1
        :param vocabList: 所有单词集合列表
        :param inputSet: 输入数据集
        :return:匹配列表[0,1,0,1...]，其中 1与0 表示词汇表中的单词是否出现在输入的数据集中
        """
        # 创建一个和词汇表等长的向量，并将其元素都设置为0
        returnVec = [0] * len(vocabList)
        # 遍历文档中的所有单词，如果出现了词汇表中的单词，则将输出的文档向量中的对应值设为1
        for word in inputSet:
            if word in vocabList:
                returnVec[vocabList.index(word)] = 1
            else:
                print("the word: %s is not in my Vocabulary!" % word)
        return returnVec

*分析数据: 检查词条确保解析的正确性*

检查函数执行情况，检查词表，不出现重复单词，需要的话，可以对其进行排序。  
![在这里插入图片描述][20191128221233424.png]  
检查函数有效性。例如：myVocabList 中索引为 2 的元素是什么单词？应该是是 help 。该单词在第一篇文档中出现了，现在检查一下看看它是否出现在第四篇文档中。  
![在这里插入图片描述][20191128221454740.png]  
*训练算法: 从词向量计算概率*

现在已经知道了一个词是否出现在一篇文档中，也知道该文档所属的类别。接下来我们重写贝叶斯准则，将之前的 x, y 替换为 w. 粗体的 w 表示这是一个向量，即它由多个值组成。在这个例子中，数值个数与词汇表中的词个数相同。  
![重写贝叶斯准则][20191128221609785.png]  
我们使用上述公式，对每个类计算该值，然后比较这两个概率值的大小。

问: 上述代码实现中，为什么没有计算P(w)？

答：根据上述公式可知，我们右边的式子等同于左边的式子，由于对于每个ci，P(w)是固定的。并且我们只需要比较左边式子值的大小来决策分类，那么我们就可以简化为通过比较右边分子值得大小来做决策分类。

首先可以通过类别 i (侮辱性留言或者非侮辱性留言)中的文档数除以总的文档数来计算概率 p(ci) 。接下来计算 p(w | ci) ，这里就要用到朴素贝叶斯假设。如果将 w 展开为一个个独立特征，那么就可以将上述概率写作 p(w0, w1, w2…wn | ci) 。这里假设所有词都互相独立，该假设也称作条件独立性假设（例如 A 和 B 两个人抛骰子，概率是互不影响的，也就是相互独立的，A 抛 2点的同时 B 抛 3 点的概率就是 1/6 \* 1/6），它意味着可以使用 p(w0 | ci)p(w1 | ci)p(w2 | ci)…p(wn | ci) 来计算上述概率，这样就极大地简化了计算的过程。

朴素贝叶斯分类器训练函数

def _trainNB0(trainMatrix, trainCategory):
        """
        训练数据原版
        :param trainMatrix: 文件单词矩阵 [[1,0,1,1,1....],[],[]...]
        :param trainCategory: 文件对应的类别[0,1,1,0....]，列表长度等于单词矩阵数，其中的1代表对应的文件是侮辱性文件，0代表不是侮辱性矩阵
        :return:
        """
        # 文件数
        numTrainDocs = len(trainMatrix)
        # 单词数
        numWords = len(trainMatrix[0])
        # 侮辱性文件的出现概率，即trainCategory中所有的1的个数，
        # 代表的就是多少个侮辱性文件，与文件的总数相除就得到了侮辱性文件的出现概率
        pAbusive = sum(trainCategory) / float(numTrainDocs)
        # 构造单词出现次数列表
        p0Num = zeros(numWords)  # [0,0,0,.....]
        p1Num = zeros(numWords)  # [0,0,0,.....]
    
        # 整个数据集单词出现总数
        p0Denom = 0.0
        p1Denom = 0.0
        for i in range(numTrainDocs):
            # 遍历所有的文件，如果是侮辱性文件，就计算此侮辱性文件中出现的侮辱性单词的个数
            if trainCategory[i] == 1:
                p1Num += trainMatrix[i]  # [0,1,1,....]->[0,1,1,...]
                p1Denom += sum(trainMatrix[i])
            else:
                # 如果不是侮辱性文件，则计算非侮辱性文件中出现的侮辱性单词的个数
                p0Num += trainMatrix[i]
                p0Denom += sum(trainMatrix[i])
        # 类别1，即侮辱性文档的[P(F1|C1),P(F2|C1),P(F3|C1),P(F4|C1),P(F5|C1)....]列表
        # 即 在1类别下，每个单词出现次数的占比
        p1Vect = p1Num / p1Denom  # [1,2,3,5]/90->[1/90,...]
        # 类别0，即正常文档的[P(F1|C0),P(F2|C0),P(F3|C0),P(F4|C0),P(F5|C0)....]列表
        # 即 在0类别下，每个单词出现次数的占比
        p0Vect = p0Num / p0Denom
        return p0Vect, p1Vect, pAbusive

*测试算法: 根据现实情况修改分类器*

在利用贝叶斯分类器对文档进行分类时，要计算多个概率的乘积以获得文档属于某个类别的概率，即计算 p(w0|1) \* p(w1|1) \* p(w2|1)。如果其中一个概率值为 0，那么最后的乘积也为 0。为降低这种影响，可以将所有词的出现数初始化为 1，并将分母初始化为 2 （取1 或 2 的目的主要是为了保证分子和分母不为0，大家可以根据业务需求进行更改）。

另一个遇到的问题是下溢出，这是由于太多很小的数相乘造成的。当计算乘积 p(w0|ci) \* p(w1|ci) \* p(w2|ci)… p(wn|ci) 时，由于大部分因子都非常小，所以程序会下溢出或者得到不正确的答案。（用 Python 尝试相乘许多很小的数，最后四舍五入后会得到 0）。一种解决办法是对乘积取自然对数。在代数中有 ln(a \* b) = ln(a) + ln(b), 于是通过求对数可以避免下溢出或者浮点数舍入导致的错误。同时，采用自然对数进行处理不会有任何损失。

下图给出了函数 f(x) 与 ln(f(x)) 的曲线。可以看出，它们在相同区域内同时增加或者减少，并且在相同点上取到极值。它们的取值虽然不同，但不影响最终结果。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 2]

def trainNB0(trainMatrix, trainCategory):
        """
        训练数据优化版本
        :param trainMatrix: 文件单词矩阵
        :param trainCategory: 文件对应的类别
        :return:
        """
        # 总文件数
        numTrainDocs = len(trainMatrix)
        # 总单词数
        numWords = len(trainMatrix[0])
        # 侮辱性文件的出现概率
        pAbusive = sum(trainCategory) / float(numTrainDocs)
        # 构造单词出现次数列表
        # p0Num 正常的统计
        # p1Num 侮辱的统计
        # 避免单词列表中的任何一个单词为0，而导致最后的乘积为0，所以将每个单词的出现次数初始化为 1
        p0Num = ones(numWords)  # [0,0......]->[1,1,1,1,1.....]
        p1Num = ones(numWords)
    
        # 整个数据集单词出现总数，2.0根据样本/实际调查结果调整分母的值（2主要是避免分母为0，当然值可以调整）
        # p0Denom 正常的统计
        # p1Denom 侮辱的统计
        p0Denom = 2.0
        p1Denom = 2.0
        for i in range(numTrainDocs):
            if trainCategory[i] == 1:
                # 累加辱骂词的频次
                p1Num += trainMatrix[i]
                # 对每篇文章的辱骂的频次 进行统计汇总
                p1Denom += sum(trainMatrix[i])
            else:
                p0Num += trainMatrix[i]
                p0Denom += sum(trainMatrix[i])
        # 类别1，即侮辱性文档的[log(P(F1|C1)),log(P(F2|C1)),log(P(F3|C1)),log(P(F4|C1)),log(P(F5|C1))....]列表
        p1Vect = log(p1Num / p1Denom)
        # 类别0，即正常文档的[log(P(F1|C0)),log(P(F2|C0)),log(P(F3|C0)),log(P(F4|C0)),log(P(F5|C0))....]列表
        p0Vect = log(p0Num / p0Denom)
        return p0Vect, p1Vect, pAbusive

*使用算法: 对社区留言板言论进行分类*

朴素贝叶斯分类函数

def classifyNB(vec2Classify, p0Vec, p1Vec, pClass1):
        """
        使用算法：
            # 将乘法转换为加法
            乘法：P(C|F1F2...Fn) = P(F1F2...Fn|C)P(C)/P(F1F2...Fn)
            加法：P(F1|C)*P(F2|C)....P(Fn|C)P(C) -> log(P(F1|C))+log(P(F2|C))+....+log(P(Fn|C))+log(P(C))
        :param vec2Classify: 待测数据[0,1,1,1,1...]，即要分类的向量
        :param p0Vec: 类别0，即正常文档的[log(P(F1|C0)),log(P(F2|C0)),log(P(F3|C0)),log(P(F4|C0)),log(P(F5|C0))....]列表
        :param p1Vec: 类别1，即侮辱性文档的[log(P(F1|C1)),log(P(F2|C1)),log(P(F3|C1)),log(P(F4|C1)),log(P(F5|C1))....]列表
        :param pClass1: 类别1，侮辱性文件的出现概率
        :return: 类别1 or 0
        """
        # 计算公式  log(P(F1|C))+log(P(F2|C))+....+log(P(Fn|C))+log(P(C))
        # 使用 NumPy 数组来计算两个向量相乘的结果，这里的相乘是指对应元素相乘，即先将两个向量中的第一个元素相乘，然后将第2个元素相乘，以此类推。
        # 我的理解是：这里的 vec2Classify * p1Vec 的意思就是将每个词与其对应的概率相关联起来
        # 可以理解为 1.单词在词汇表中的条件下，文件是good 类别的概率 也可以理解为 2.在整个空间下，文件既在词汇表中又是good类别的概率
        p1 = sum(vec2Classify * p1Vec) + log(pClass1)
        p0 = sum(vec2Classify * p0Vec) + log(1.0 - pClass1)
        if p1 > p0:
            return 1
        else:
            return 0
    
    
    def testingNB():
        """
        测试朴素贝叶斯算法
        """
        # 1. 加载数据集
        listOPosts, listClasses = loadDataSet()
        # 2. 创建单词集合
        myVocabList = createVocabList(listOPosts)
        # 3. 计算单词是否出现并创建数据矩阵
        trainMat = []
        for postinDoc in listOPosts:
            # 返回m*len(myVocabList)的矩阵， 记录的都是0，1信息
            trainMat.append(setOfWords2Vec(myVocabList, postinDoc))
        # 4. 训练数据
        p0V, p1V, pAb = trainNB0(array(trainMat), array(listClasses))
        # 5. 测试数据
        testEntry = ['love', 'my', 'dalmation']
        thisDoc = array(setOfWords2Vec(myVocabList, testEntry))
        print(testEntry, 'classified as: ', classifyNB(thisDoc, p0V, p1V, pAb))
        testEntry = ['stupid', 'garbage']
        thisDoc = array(setOfWords2Vec(myVocabList, testEntry))
        print(testEntry, 'classified as: ', classifyNB(thisDoc, p0V, p1V, pAb))

#### 项目案例2: 使用朴素贝叶斯过滤垃圾邮件 ####

**项目概述**  
完成朴素贝叶斯的一个最著名的应用: 电子邮件垃圾过滤。

**开发流程**  
使用朴素贝叶斯对电子邮件进行分类

> 收集数据: 提供文本文件  
> 准备数据: 将文本文件解析成词条向量  
> 分析数据: 检查词条确保解析的正确性  
> 训练算法: 使用我们之前建立的 trainNB() 函数  
> 测试算法: 使用朴素贝叶斯进行交叉验证  
> 使用算法: 构建一个完整的程序对一组文档进行分类，将错分的文档输出到屏幕上

*收集数据: 提供文本文件*

文本文件内容如下:

> Hi Peter,
> 
> With Jose out of town, do you want to meet once in a while to keep  
> things going and do some interesting stuff?
> 
> Let me know Eugene

*准备数据: 将文本文件解析成词条向量*

def textParse(bigString):
        """
        Desc:
            接收一个大字符串并将其解析为字符串列表
        Args:
            bigString -- 大字符串
        Returns:
            去掉少于 2 个字符的字符串，并将所有字符串转换为小写，返回字符串列表
        """
        import re
        # 使用正则表达式来切分句子，其中分隔符是除单词、数字外的任意字符串
        listOfTokens = re.split(r'\W+', bigString)
        return [tok.lower() for tok in listOfTokens if len(tok) > 2]

*分析数据: 检查词条确保解析的正确性*

*训练算法: 使用我们之前建立的 trainNB0() 函数*

*测试算法: 使用朴素贝叶斯进行交叉验证*

完整的垃圾邮件测试函数

def spamTest():
        """
        Desc:
            对贝叶斯垃圾邮件分类器进行自动化处理。
        Returns:
            对测试集中的每封邮件进行分类，若邮件分类错误，则错误数加 1，最后返回总的错误百分比。
        """
        docList = []
        classList = []
        fullText = []
        for i in range(1, 26):
            # 切分，解析数据，并归类为 1 类别
            wordList = textParse(open('data/email/spam/%d.txt' % i).read())
            docList.append(wordList)
            classList.append(1)
            # 切分，解析数据，并归类为 0 类别
            wordList = textParse(open('data/email/ham/%d.txt' % i).read())
            docList.append(wordList)
            fullText.extend(wordList)
            classList.append(0)
        # 创建词汇表
        vocabList = createVocabList(docList)
        trainingSet = list(range(50))
        testSet = []
        # 随机取 10 个邮件用来测试
        for i in range(10):
            # random.uniform(x, y) 随机生成一个范围为 x - y 的实数
            randIndex = int(random.uniform(0, len(trainingSet)))
            testSet.append(trainingSet[randIndex])
            del (trainingSet[randIndex])
        trainMat = []
        trainClasses = []
        for docIndex in trainingSet:
            trainMat.append(setOfWords2Vec(vocabList, docList[docIndex]))
            trainClasses.append(classList[docIndex])
        p0V, p1V, pSpam = trainNB0(array(trainMat), array(trainClasses))
        errorCount = 0
        for docIndex in testSet:
            wordVector = setOfWords2Vec(vocabList, docList[docIndex])
            if classifyNB(array(wordVector), p0V, p1V, pSpam) != classList[docIndex]:
                errorCount += 1
        print('the errorCount is: ', errorCount)
        print('the testSet length is :', len(testSet))
        print('the error rate is :', float(errorCount) / len(testSet))

*使用算法: 构建一个完整的程序对一组文档进行分类*

![在这里插入图片描述][20191129164311554.png]

#### 项目案例3: 使用朴素贝叶斯分类器从个人广告中获取区域倾向 ####

**项目概述**

广告商往往想知道关于一个人的一些特定人口统计信息，以便能更好地定向推销广告。

我们将分别从美国的两个城市中选取一些人，通过分析这些人发布的信息，来比较这两个城市的人们在广告用词上是否不同。如果结论确实不同，那么他们各自常用的词是哪些，从人们的用词当中，我们能否对不同城市的人所关心的内容有所了解。

**开发流程**

> 收集数据: 从 RSS 源收集内容，这里需要对 RSS 源构建一个接口  
> 准备数据: 将文本文件解析成词条向量  
> 分析数据: 检查词条确保解析的正确性  
> 训练算法: 使用我们之前建立的 trainNB0() 函数  
> 测试算法: 观察错误率，确保分类器可用。可以修改切分程序，以降低错误率，提高分类结果  
> 使用算法: 构建一个完整的程序，封装所有内容。给定两个 RSS 源，改程序会显示最常用的公共词

*收集数据: 从 RSS 源收集内容，这里需要对 RSS 源构建一个接口*

也就是导入 RSS 源，使用 pip 下载RSS阅读器，即 feedparser。  
同时，由于文中提到的 RSS源 不可用，所以替换了差不多的两个`https://newyork.craigslist.org/search/res?format=rss`  
`https://sfbay.craigslist.org/search/apa?format=rss`

*准备数据: 将文本文件解析成词条向量*

**文档词袋模型**  
我们将每个词的出现与否作为一个特征，这可以被描述为 **词集模型(set-of-words model)**。如果一个词在文档中出现不止一次，这可能意味着包含该词是否出现在文档中所不能表达的某种信息，这种方法被称为词袋模型(bag-of-words model)。在词袋中，每个单词可以出现多次，而在词集中，每个词只能出现一次。为适应词袋模型，需要对函数 setOfWords2Vec() 稍加修改，修改后的函数为 bagOfWords2Vec() 。

如下给出了基于词袋模型的朴素贝叶斯代码。它与函数 setOfWords2Vec() 几乎完全相同，唯一不同的是每当遇到一个单词时，它会增加词向量中的对应值，而不只是将对应的数值设为 1 。

def setOfWords2VecMN(vocabList, inputSet):
        returnVec = [0] * len(vocabList)  # 创建一个其中所含元素都为0的向量
        for word in inputSet:
            if word in vocabList:
                returnVec[vocabList.index(word)] += 1
        return returnVec
    
    
    # 文件解析
    def textParse(bigString):
        import re
        listOfTokens = re.split(r'\W+', bigString)  # 注意这里*改为+号，否则会将单词分成一个一个的字母，详尽可查看正则表达式
        return [tok.lower() for tok in listOfTokens if len(tok) > 2]

*分析数据: 检查词条确保解析的正确性*

*训练算法: 使用我们之前建立的 trainNB0() 函数*

*测试算法: 观察错误率，确保分类器可用。可以修改切分程序，以降低错误率，提高分类结果*

# RSS源分类器及高频词去除函数
    def calcMostFreq(vocabList, fullText):
        import operator
        freqDict = {
        }
        for token in vocabList:  # 遍历词汇表中的每个词
            freqDict[token] = fullText.count(token)  # 统计每个词在文本中出现的次数
        sortedFreq = sorted(freqDict.items(), key=operator.itemgetter(1), reverse=True)  # 根据每个词出现的次数从高到底对字典进行排序
        return sortedFreq[:30]  # 返回出现次数最高的30个单词
    
    
    def localWords(feed1, feed0):
        import random
        import feedparser
        docList = []
        classList = []
        fullText = []
        minLen = min(len(feed1['entries']), len(feed0['entries']))
        for i in range(minLen):
            wordList = textParse(feed1['entries'][i]['summary'])  # 每次访问一条RSS源
            docList.append(wordList)
            fullText.extend(wordList)
            classList.append(1)
            wordList = textParse(feed0['entries'][i]['summary'])
            docList.append(wordList)
            fullText.extend(wordList)
            classList.append(0)
        vocabList = createVocabList(docList)
        top30Words = calcMostFreq(vocabList, fullText)
        for pairW in top30Words:
            if pairW[0] in vocabList:
                vocabList.remove(pairW[0])  # 去掉出现次数最高的那些词
        trainingSet = list(range(2 * minLen))
        testSet = []
        for i in range(5):
            randIndex = int(random.uniform(0, len(trainingSet)))
            testSet.append(trainingSet[randIndex])
            del (trainingSet[randIndex])
        trainMat = []
        trainClasses = []
        for docIndex in trainingSet:
            trainMat.append(bagOfWords2VecMN(vocabList, docList[docIndex]))
            trainClasses.append(classList[docIndex])
        p0V, p1V, pSpam = trainNB0(array(trainMat), array(trainClasses))
        errorCount = 0
        for docIndex in testSet:
            wordVector = bagOfWords2VecMN(vocabList, docList[docIndex])
            if classifyNB(array(wordVector), p0V, p1V, pSpam) != classList[docIndex]:
                errorCount += 1
        print('the error rate is:', float(errorCount) / len(testSet))
        return vocabList, p0V, p1V

*使用算法: 构建一个完整的程序，封装所有内容。给定两个 RSS 源，改程序会显示最常用的公共词*

函数 localWords() 使用了两个 RSS 源作为参数，RSS 源要在函数外导入，这样做的原因是 RSS 源会随时间而改变，重新加载 RSS 源就会得到新的数据  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 3]  
为了得到错误率的精确估计，应该多次进行上述实验，然后取平均值

接下来，我们要分析一下数据，显示地域相关的用词

可以先对向量pSF与pNY进行排序，然后按照顺序打印出来，将下面的代码添加到文件中：

# 最具表征性的词汇显示函数
    def getTopWords(ny, sf):
        import operator
        vocabList, p0V, p1V = localWords(ny, sf)
        topNY = []
        topSF = []
        for i in range(len(p0V)):
            if p0V[i] > -6.0:
                topSF.append((vocabList[i], p0V[i]))
            if p1V[i] > -6.0:
                topNY.append((vocabList[i], p1V[i]))
        sortedSF = sorted(topSF, key=lambda pair: pair[1], reverse=True)
        print("SF**SF**SF**SF**SF**SF**SF**SF**SF**SF**SF**SF**SF**SF**")
        for item in sortedSF:
            print(item[0])
        sortedNY = sorted(topNY, key=lambda pair: pair[1], reverse=True)
        print("NY**NY**NY**NY**NY**NY**NY**NY**NY**NY**NY**NY**NY**NY**")
        for item in sortedNY:
            print(item[0])

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 4]  
从最后输出的单词，可以看出程序输出了大量的停用词，可以移除固定的停用词看看结果如何，这样做的话，分类错误率也会降低。

### 5、朴素贝叶斯 小结 ###

> 对于分类而言，使用概率有时要比使用硬规则更为有效。贝叶斯概率及贝叶斯准则提供了一 种利用已知值来估计未知概率的有效方法。  
> 可以通过特征之间的条件独立性假设，降低对数据量的需求。独立性假设是指一个词的出现 概率并不依赖于文档中的其他词。当然我们也知道这个假设过于简单。这就是之所以称为朴素贝 叶斯的原因。尽管条件独立性假设并不正确，但是朴素贝叶斯仍然是一种有效的分类器。  
> 利用现代编程语言来实现朴素贝叶斯时需要考虑很多实际因素。下溢出就是其中一个问题， 它可以通过对概率取对数来解决。词袋模型在解决文档分类问题上比词集模型有所提高。还有其 他一些方面的改进，比如说移除停用词，当然也可以花大量时间对切分器进行优化。

## 二、Logistic回归 ##

### 1、Logistic 回归 概述 ###

> Logistic 回归 或者叫逻辑回归 虽然名字有回归，但是它是用来做分类的。其主要思想是: 根据现有数据对分类边界线(Decision Boundary)建立回归公式，以此进行分类。

#### 须知概念 ####

**Sigmoid 函数**

回归 概念

假设现在有一些数据点，我们用一条直线对这些点进行拟合（这条直线称为最佳拟合直线），这个拟合的过程就叫做回归。进而可以得到对这些点的拟合直线方程，那么我们根据这个回归方程，怎么进行分类呢？请看下面。

二值型输出分类函数

我们想要的函数应该是: 能接受所有的输入然后预测出类别。例如，在两个类的情况下，上述函数输出 0 或 1.或许你之前接触过具有这种性质的函数，该函数称为 海维塞得阶跃函数(Heaviside step function)，或者直接称为 单位阶跃函数。然而，海维塞得阶跃函数的问题在于: 该函数在跳跃点上从 0 瞬间跳跃到 1，这个瞬间跳跃过程有时很难处理。幸好，另一个函数也有类似的性质（可以输出 0 或者 1 的性质），且数学上更易处理，这就是 Sigmoid 函数。 Sigmoid 函数具体的计算公式如下:  
![Sigmoid 函数计算公式][Sigmoid]  
下图给出了 Sigmoid 函数在不同坐标尺度下的两条曲线图。当 x 为 0 时，Sigmoid 函数值为 0.5 。随着 x 的增大，对应的 Sigmoid 值将逼近于 1 ; 而随着 x 的减小， Sigmoid 值将逼近于 0 。如果横坐标刻度足够大， Sigmoid 函数看起来很像一个阶跃函数。  
![Sigmoid 函数在不同坐标下的图片][Sigmoid 1]  
因此，为了实现 Logistic 回归分类器，我们可以在每个特征上都乘以一个回归系数（如下公式所示），然后把所有结果值相加，将这个总和代入 Sigmoid 函数中，进而得到一个范围在 0~1 之间的数值。任何大于 0.5 的数据被分入 1 类，小于 0.5 即被归入 0 类。所以，Logistic 回归也是一种概率估计，比如这里Sigmoid 函数得出的值为0.5，可以理解为给定数据和参数，数据被分入 1 类的概率为0.5。

**基于最优化方法的回归系数确定**  
Sigmoid 函数的输入记为 z ，由下面公式得到:  
![Sigmoid 函数计算公式][Sigmoid 2]  
如果采用向量的写法，上述公式可以写成 Sigmoid 函数计算公式向量形式 ，它表示将这两个数值向量对应元素相乘然后全部加起来即得到 z 值。其中的向量 x 是分类器的输入数据，向量 w 也就是我们要找到的最佳参数（系数），从而使得分类器尽可能地精确。为了寻找该最佳参数，需要用到最优化理论的一些知识。我们这里使用的是——梯度上升法（Gradient Ascent）。

**梯度上升法**  
梯度的介绍

> 向量 = 值 + 方向  
> 梯度 = 向量  
> 梯度 = 梯度值 + 梯度方向

梯度上升法的思想

要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。如果梯度记为 ▽ ，则函数 f(x, y) 的梯度由下式表示:  
![梯度上升计算公式][20191129173513504.png]  
这个梯度意味着要沿 x 的方向移动 ![f(x, y)对x求偏导][f_x_ y_x] ，沿 y 的方向移动 ![f(x, y)对y求偏导][f_x_ y_y] 。其中，函数f(x, y) 必须要在待计算的点上有定义并且可微。下图是一个具体的例子。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 5]  
上图展示的，梯度上升算法到达每个点后都会重新估计移动的方向。从 P0 开始，计算完该点的梯度，函数就根据梯度移动到下一点 P1。在 P1 点，梯度再次被重新计算，并沿着新的梯度方向移动到 P2 。如此循环迭代，直到满足停止条件。迭代过程中，梯度算子总是保证我们能选取到最佳的移动方向。

上图中的梯度上升算法沿梯度方向移动了一步。可以看到，梯度算子总是指向函数值增长最快的方向。这里所说的是移动方向，而未提到移动量的大小。该量值称为步长，记作 α 。用向量来表示的话，梯度上升算法的迭代公式如下:  
![在这里插入图片描述][20191129173715281.png]  
该公式将一直被迭代执行，直至达到某个停止条件为止，比如迭代次数达到某个指定值或者算法达到某个可以允许的误差范围。

介绍一下几个相关的概念：

> 例如：y = w0 + w1x1 + w2x2 + … + wnxn 梯度：参考上图的例子，二维图像，x方向代表第一个系数，也就是  
> w1，y方向代表第二个系数也就是 w2，这样的向量就是梯度。 α：上面的梯度算法的迭代公式中的阿尔法，这个代表的是移动步长（step  
> length）。移动步长会影响最终结果的拟合程度，最好的方法就是随着迭代次数更改移动步长。  
> 步长通俗的理解，100米，如果我一步走10米，我需要走10步；如果一步走20米，我只需要走5步。这里的一步走多少米就是步长的意思。  
> ▽f(w)：代表沿着梯度变化的方向。

问：有人会好奇为什么有些书籍上说的是梯度下降法（Gradient Decent）?

答： 其实这个两个方法在此情况下本质上是相同的。关键在于代价函数（cost function）或者叫目标函数（objective function）。如果目标函数是损失函数，那就是最小化损失函数来求函数的最小值，就用梯度下降。 如果目标函数是似然函数（Likelihood function），就是要最大化似然函数来求函数的最大值，那就用梯度上升。在逻辑回归中， 损失函数和似然函数无非就是互为正负关系。

只需要在迭代公式中的加法变成减法。因此，对应的公式可以写成  
![在这里插入图片描述][20191129173827306.png]  
**局部最优现象 （Local Optima）**  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 6]  
上图表示参数 θ 与误差函数 J(θ) 的关系图 (这里的误差函数是损失函数，所以我们要最小化损失函数)，红色的部分是表示 J(θ) 有着比较高的取值，我们需要的是，能够让 J(θ) 的值尽量的低。也就是深蓝色的部分。θ0，θ1 表示 θ 向量的两个维度（此处的θ0，θ1是x0和x1的系数，也对应的是上文w0和w1）。

可能梯度下降的最终点并非是全局最小点，可能是一个局部最小点，如我们上图中的右边的梯度下降曲线，描述的是最终到达一个局部最小点，这是我们重新选择了一个初始点得到的。

看来我们这个算法将会在很大的程度上被初始点的选择影响而陷入局部最小点。

### 2、Logistic 回归 原理 ###

#### 工作原理 ####

> 每个回归系数初始化为 1 重复 R 次:  
> 计算整个数据集的梯度  
> 使用 步长 x 梯度 更新回归系数的向量 返回回归系数

#### 开发流程 ####

> 收集数据: 采用任意方法收集数据  
> 准备数据: 由于需要进行距离计算，因此要求数据类型为数值型。另外，结构化数据格式则最佳。  
> 分析数据: 采用任意方法对数据进行分析。  
> 训练算法: 大部分时间将用于训练，训练的目的是为了找到最佳的分类回归系数。  
> 测试算法: 一旦训练步骤完成，分类将会很快。  
> 使用算法: 首先，我们需要输入一些数据，并将其转换成对应的结构化数值；接着，基于训练好的回归系数就可以对这些数值进行简单的回归计算，判定它们属于哪个类别；在这之后，我们就可以在输出的类别上做一些其他分析工作。

#### 算法特点 ####

> 优点: 计算代价不高，易于理解和实现。  
> 缺点: 容易欠拟合，分类精度可能不高。  
> 适用数据类型: 数值型和标称型数据。

### 3、Logistic 回归 项目案例 ###

完整代码地址：[https://blog.csdn.net/qq\_45556599/article/details/103315793][https_blog.csdn.net_qq_45556599_article_details_103315793]

#### 项目案例1: 使用 Logistic 回归在简单数据集上的分类 ####

**项目概述**  
在一个简单的数据集上，采用梯度上升法找到 Logistic 回归分类器在此数据集上的最佳回归系数

**开发流程**

> 收集数据: 可以使用任何方法  
> 准备数据: 由于需要进行距离计算，因此要求数据类型为数值型。另外，结构化数据格式则最佳  
> 分析数据: 画出决策边界  
> 训练算法: 使用梯度上升找到最佳参数  
> 测试算法: 使用 Logistic 回归进行分类  
> 使用算法: 对简单数据集中数据进行分类

*收集数据: 可以使用任何方法*

我们采用存储在 TestSet.txt 文本文件中的数据，存储格式如下:

-0.017612	14.053064	0
    -1.395634	4.662541	1
    -0.752157	6.538620	0
    -1.322371	7.152853	0
    0.423363	11.054677	0

绘制在图中，如下图所示:  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 7]  
*准备数据: 由于需要进行距离计算，因此要求数据类型为数值型。另外，结构化数据格式则最佳*

def loadDataSet(file_name):
        """
        Desc:
            加载并解析数据
        Args:
            file_name -- 文件名称，要解析的文件所在磁盘位置
        Returns:
            dataMat -- 原始数据的特征
            labelMat -- 原始数据的标签，也就是每条样本对应的类别
        """
        # dataMat为原始数据， labelMat为原始数据的标签
        dataMat = []
        labelMat = []
        fr = open(file_name)
        for line in fr.readlines():
            lineArr = line.strip().split()
            if len(lineArr) == 1:
                continue    # 这里如果就一个空的元素，则跳过本次循环
            # 为了方便计算，我们将 X0 的值设为 1.0 ，也就是在每一行的开头添加一个 1.0 作为 X0
            dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
            labelMat.append(int(lineArr[2]))
        return dataMat, labelMat

*训练算法: 使用梯度上升找到最佳参数*

# sigmoid跳跃函数
    def sigmoid(inX):
        # return 1.0 / (1 + exp(-inX))
    
        # Tanh是Sigmoid的变形，与 sigmoid 不同的是，tanh 是0均值的。因此，实际应用中，tanh 会比 sigmoid 更好。
        return 2 * 1.0/(1+exp(-2*inX)) - 1
    
    
    # 正常的处理方案
    # 两个参数：第一个参数==> dataMatIn 是一个2维NumPy数组，每列分别代表每个不同的特征，每行则代表每个训练样本。
    # 第二个参数==> classLabels 是类别标签，它是一个 1*100 的行向量。为了便于矩阵计算，需要将该行向量转换为列向量，做法是将原向量转置，再将它赋值给labelMat。
    def gradAscent(dataMatIn, classLabels):
        """
        Desc:
            正常的梯度上升法
        Args:
            dataMatIn -- 输入的 数据的特征 List
            classLabels -- 输入的数据的类别标签
        Returns:
            array(weights) -- 得到的最佳回归系数
        """
    
        # 转化为矩阵[[1,1,2],[1,1,2]....]
        dataMatrix = mat(dataMatIn)  # 转换为 NumPy 矩阵
        # 转化为矩阵[[0,1,0,1,0,1.....]]，并转制[[0],[1],[0].....]
        # transpose() 行列转置函数
        # 将行向量转化为列向量   =>  矩阵的转置
        labelMat = mat(classLabels).transpose()  # 首先将数组转换为 NumPy 矩阵，然后再将行向量转置为列向量
        # m->数据量，样本数 n->特征数
        m, n = shape(dataMatrix)
        # print m, n, '__'*10, shape(dataMatrix.transpose()), '__'*100
        # alpha代表向目标移动的步长
        alpha = 0.001
        # 迭代次数
        maxCycles = 500
        # 生成一个长度和特征数相同的矩阵，此处n为3 -> [[1],[1],[1]]
        # weights 代表回归系数， 此处的 ones((n,1)) 创建一个长度和特征数相同的矩阵，其中的数全部都是 1
        weights = ones((n, 1))
        for k in range(maxCycles):  # heavy on matrix operations
            # m*3 的矩阵 * 3*1 的单位矩阵 ＝ m*1的矩阵
            # 那么乘上单位矩阵的意义，就代表：通过公式得到的理论值
            # 参考地址： 矩阵乘法的本质是什么？ https://www.zhihu.com/question/21351965/answer/31050145
            # print 'dataMatrix====', dataMatrix
            # print 'weights====', weights
            # n*3   *  3*1  = n*1
            h = sigmoid(dataMatrix * weights)  # 矩阵乘法
            # print 'hhhhhhh====', h
            # labelMat是实际值
            error = (labelMat - h)  # 向量相减
            # 0.001* (3*m)*(m*1) 表示在每一个列上的一个误差情况，最后得出 x1,x2,xn的系数的偏移量
            weights = weights + alpha * dataMatrix.transpose() * error  # 矩阵乘法，最后得到回归系数
        return array(weights)

*分析数据: 采用任意方法对数据进行分析*

画出数据集和 Logistic 回归最佳拟合直线的函数

# 可视化展示
    def plotBestFit(dataArr, labelMat, weights):
        """
            Desc:
                将我们得到的数据可视化展示出来
            Args:
                dataArr:样本数据的特征
                labelMat:样本数据的类别标签，即目标变量
                weights:回归系数
            Returns:
                None
        """
        n = shape(dataArr)[0]
        xcord1 = []
        ycord1 = []
        xcord2 = []
        ycord2 = []
        for i in range(n):
            if int(labelMat[i]) == 1:
                xcord1.append(dataArr[i, 1])
                ycord1.append(dataArr[i, 2])
            else:
                xcord2.append(dataArr[i, 1])
                ycord2.append(dataArr[i, 2])
        fig = plt.figure()
        ax = fig.add_subplot(111)
        ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
        ax.scatter(xcord2, ycord2, s=30, c='green')
        x = arange(-3.0, 3.0, 0.1)
        """
        y的由来，卧槽，是不是没看懂？
        首先理论上是这个样子的。
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        w0*x0+w1*x1+w2*x2=f(x)
        x0最开始就设置为1叻， x2就是我们画图的y值，而f(x)被我们磨合误差给算到w0,w1,w2身上去了
        所以： w0+w1*x+w2*y=0 => y = (-w0-w1*x)/w2   
        """
        y = (-weights[0] - weights[1] * x) / weights[2]
        ax.plot(x, y)
        plt.xlabel('X')
        plt.ylabel('Y')
        plt.show()

*测试算法: 使用 Logistic 回归进行分类*

def simpleTest():
        # 1.收集并准备数据
        dataMat, labelMat = loadDataSet("data/TestSet.txt")
    
        # print dataMat, '---\n', labelMat
        # 2.训练模型，  f(x)=a1*x1+b2*x2+..+nn*xn中 (a1,b2, .., nn).T的矩阵值
        # 因为数组没有是复制n份， array的乘法就是乘法
        dataArr = array(dataMat)
        # print dataArr
        weights = gradAscent(dataArr, labelMat)
        # weights = stocGradAscent0(dataArr,
        # labelMat)
        # weights = stocGradAscent1(dataArr, labelMat)
        # print '*'*30, weights
    
        # 数据可视化
        plotBestFit(dataArr, labelMat, weights)

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 8]  
*训练算法：随机梯度上升*  
**注意**  
梯度上升算法在每次更新回归系数时都需要遍历整个数据集，该方法在处理 100 个左右的数据集时尚可，但如果有数十亿样本和成千上万的特征，那么该方法的计算复杂度就太高了。一种改进方法是一次仅用一个样本点来更新回归系数，该方法称为 随机梯度上升算法。由于可以在新样本到来时对分类器进行增量式更新，因而随机梯度上升算法是一个在线学习(online learning)算法。与 “在线学习” 相对应，一次处理所有数据被称作是 “批处理” （batch） 。

随机梯度上升算法可以写成如下的伪代码:

所有回归系数初始化为 1
    对数据集中每个样本
        计算该样本的梯度
        使用 alpha x gradient 更新回归系数值
    返回回归系数值

以下是随机梯度上升算法的实现代码:

# 随机梯度下降
    # 梯度下降优化算法在每次更新数据集时都需要遍历整个数据集，计算复杂都较高
    # 随机梯度下降一次只用一个样本点来更新回归系数
    def stocGradAscent0(dataMatrix, classLabels):
        """
        Desc:
            随机梯度下降，只使用一个样本点来更新回归系数
        Args:
            dataMatrix -- 输入数据的数据特征（除去最后一列）
            classLabels -- 输入数据的类别标签（最后一列数据）
        Returns:
            weights -- 得到的最佳回归系数
        """
        m, n = shape(dataMatrix)
        alpha = 0.01
        # n*1的矩阵
        # 函数ones创建一个全1的数组
        weights = ones(n)  # 初始化长度为n的数组，元素全部为 1
        for i in range(m):
            # sum(dataMatrix[i]*weights)为了求 f(x)的值， f(x)=a1*x1+b2*x2+..+nn*xn,此处求出的 h 是一个具体的数值，而不是一个矩阵
            h = sigmoid(sum(dataMatrix[i] * weights))
            # print 'dataMatrix[i]===', dataMatrix[i]
            # 计算真实类别与预测类别之间的差值，然后按照该差值调整回归系数
            error = classLabels[i] - h
            # 0.01*(1*1)*(1*n)
            # print weights, "*" * 10, dataMatrix[i], "*" * 10, error
            weights = weights + alpha * error * dataMatrix[i]
        return weights

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 9]  
可以看到，随机梯度上升算法与梯度上升算法在代码上很相似，但也有一些区别: 第一，后者的变量 h 和误差 error 都是向量，而前者则全是数值；第二，前者没有矩阵的转换过程，所有变量的数据类型都是 NumPy 数组。

判断优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到了稳定值，是否还会不断地变化？下图展示了随机梯度上升算法在 200 次迭代过程中回归系数的变化情况。其中的系数2，也就是 X2 只经过了 50 次迭代就达到了稳定值，但系数 1 和 0 则需要更多次的迭代。如下图所示:  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 10]

针对这个问题，我们改进了之前的随机梯度上升算法，如下:

# 随机梯度下降算法（随机化）
    def stocGradAscent1(dataMatrix, classLabels, numIter=150):
        """
        Desc:
            改进版的随机梯度下降，使用随机的一个样本来更新回归系数
        Args:
            dataMatrix -- 输入数据的数据特征（除去最后一列数据）
            classLabels -- 输入数据的类别标签（最后一列数据）
            numIter=150 --  迭代次数
        Returns:
            weights -- 得到的最佳回归系数
        """
        m, n = shape(dataMatrix)
        weights = ones(n)  # 创建与列数相同的矩阵的系数矩阵，所有的元素都是1
        # 随机梯度, 循环150,观察是否收敛
        for j in range(numIter):
            # [0, 1, 2 .. m-1]
            dataIndex = list(range(m))
            for i in range(m):
                # i和j的不断增大，导致alpha的值不断减少，但是不为0
                alpha = 4 / (
                    1.0 + j + i
                ) + 0.0001  # alpha 会随着迭代不断减小，但永远不会减小到0，因为后边还有一个常数项0.0001
                # 随机产生一个 0～len()之间的一个值
                # random.uniform(x, y) 方法将随机生成下一个实数，它在[x,y]范围内,x是这个范围内的最小值，y是这个范围内的最大值。
                randIndex = int(random.uniform(0, len(dataIndex)))
                # sum(dataMatrix[i]*weights)为了求 f(x)的值， f(x)=a1*x1+b2*x2+..+nn*xn
                h = sigmoid(sum(dataMatrix[dataIndex[randIndex]] * weights))
                error = classLabels[dataIndex[randIndex]] - h
                # print weights, '__h=%s' % h, '__'*20, alpha, '__'*20, error, '__'*20, dataMatrix[randIndex]
                weights = weights + alpha * error * dataMatrix[dataIndex[randIndex]]
                del (dataIndex[randIndex])
        return weights

上面的改进版随机梯度上升算法，我们修改了两处代码。

第一处改进为 alpha 的值。alpha 在每次迭代的时候都会调整，这回缓解上面波动图的数据波动或者高频波动。另外，虽然 alpha 会随着迭代次数不断减少，但永远不会减小到 0，因为我们在计算公式中添加了一个常数项。

第二处修改为 randIndex 更新，这里通过随机选取样本拉来更新回归系数。这种方法将减少周期性的波动。这种方法每次随机从列表中选出一个值，然后从列表中删掉该值（再进行下一次迭代）。

程序运行之后能看到类似于下图的结果图。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 11]

#### 项目案例2: 从疝气病症预测病马的死亡率 ####

**项目概述**  
使用 Logistic 回归来预测患有疝病的马的存活问题。疝病是描述马胃肠痛的术语。然而，这种病不一定源自马的胃肠问题，其他问题也可能引发马疝病。这个数据集中包含了医院检测马疝病的一些指标，有的指标比较主观，有的指标难以测量，例如马的疼痛级别。

**开发流程**

> 收集数据: 给定数据文件  
> 准备数据: 用 Python 解析文本文件并填充缺失值  
> 分析数据: 可视化并观察数据 训练算法:使用优化算法，找到最佳的系数  
> 测试算法: 为了量化回归的效果，需要观察错误率。根据错误率决定是否回退到训练阶段，通过改变迭代的次数和步长的参数来得到更好的回归系数  
> 使用算法: 实现一个简单的命令行程序来收集马的症状并输出预测结果并非难事，这可以作为留给大家的一道习题

*收集数据: 给定数据文件*  
病马的训练数据已经给出来了，如下形式存储在文本文件中:

1.000000	1.000000	39.200000	88.000000	20.000000	0.000000	0.000000	4.000000	1.000000	3.000000	4.000000	2.000000	0.000000	0.000000	0.000000	4.000000	2.000000	50.000000	85.000000	2.000000	2.000000	0.000000
    2.000000	1.000000	38.300000	40.000000	24.000000	1.000000	1.000000	3.000000	1.000000	3.000000	3.000000	1.000000	0.000000	0.000000	0.000000	1.000000	1.000000	33.000000	6.700000	0.000000	0.000000	1.000000

*准备数据: 用 Python 解析文本文件并填充缺失值*

处理数据中的缺失值

假设有100个样本和20个特征，这些数据都是机器收集回来的。若机器上的某个传感器损坏导致一个特征无效时该怎么办？此时是否要扔掉整个数据？这种情况下，另外19个特征怎么办？ 它们是否还可以用？答案是肯定的。因为有时候数据相当昂贵，扔掉和重新获取都是不可取的，所以必须采用一些方法来解决这个问题。

下面给出了一些可选的做法：

*  使用可用特征的均值来填补缺失值；
 *  使用特殊值来填补缺失值，如 -1；
 *  忽略有缺失值的样本；
 *  使用有相似样本的均值添补缺失值；
 *  使用另外的机器学习算法预测缺失值。

现在，我们对下一节要用的数据集进行预处理，使其可以顺利地使用分类算法。在预处理需要做两件事:

*  所有的缺失值必须用一个实数值来替换，因为我们使用的 NumPy 数据类型不允许包含缺失值。我们这里选择实数 0 来替换所有缺失值，恰好能适用于 Logistic 回归。这样做的直觉在于，我们需要的是一个在更新时不会影响系数的值。回归系数的更新公式如下:  
    weights = weights + alpha \* error \* dataMatrix\[dataIndex\[randIndex\]\]  
    如果 dataMatrix 的某个特征对应值为 0，那么该特征的系数将不做更新，即:  
    weights = weights  
    另外，由于 Sigmoid(0) = 0.5 ，即它对结果的预测不具有任何倾向性，因此我们上述做法也不会对误差造成任何影响。基于上述原因，将缺失值用 0 代替既可以保留现有数据，也不需要对优化算法进行修改。此外，该数据集中的特征取值一般不为 0，因此在某种意义上说它也满足 “特殊值” 这个要求。
 *  如果在测试数据集中发现了一条数据的类别标签已经缺失，那么我们的简单做法是将该条数据丢弃。这是因为类别标签与特征不同，很难确定采用某个合适的值来替换。采用 Logistic 回归进行分类时这种做法是合理的，而如果采用类似 kNN 的方法，则保留该条数据显得更加合理。

原始的数据集经过预处理后，保存成两个文件: horseColicTest.txt 和 horseColicTraining.txt 。

*分析数据: 可视化并观察数据*

将数据使用 MatPlotlib 打印出来，观察数据是否是我们想要的格式

*训练算法: 使用优化算法，找到最佳的系数*

*测试算法：用 Logistic回归进行分类*

# 分类函数，根据回归系数和特征向量来计算 Sigmoid的值
    def classifyVector(inX, weights):
        """
        Desc:
            最终的分类函数，根据回归系数和特征向量来计算 Sigmoid 的值，大于0.5函数返回1，否则返回0
        Args:
            inX -- 特征向量，features
            weights -- 根据梯度下降/随机梯度下降 计算得到的回归系数
        Returns:
            如果 prob 计算大于 0.5 函数返回 1
            否则返回 0
        """
        prob = sigmoid(sum(inX * weights))
        if prob > 0.5:
            return 1.0
        else:
            return 0.0
    
    
    # 打开测试集和训练集,并对数据进行格式化处理
    def colicTest():
        """
        Desc:
            打开测试集和训练集，并对数据进行格式化处理
        Args:
            None
        Returns:
            errorRate -- 分类错误率
        """
        frTrain = open('data/horseColicTraining.txt')
        frTest = open('data/horseColicTest.txt')
        trainingSet = []
        trainingLabels = []
        # 解析训练数据集中的数据特征和Labels
        # trainingSet 中存储训练数据集的特征，trainingLabels 存储训练数据集的样本对应的分类标签
        for line in frTrain.readlines():
            currLine = line.strip().split('\t')
            lineArr = []
            for i in range(21):
                lineArr.append(float(currLine[i]))
            trainingSet.append(lineArr)
            trainingLabels.append(float(currLine[21]))
        # 使用 改进后的 随机梯度下降算法 求得在此数据集上的最佳回归系数 trainWeights
        trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 500)
        # trainWeights = stocGradAscent0(array(trainingSet), trainingLabels)
        errorCount = 0
        numTestVec = 0.0
        # 读取 测试数据集 进行测试，计算分类错误的样本条数和最终的错误率
        for line in frTest.readlines():
            numTestVec += 1.0
            currLine = line.strip().split('\t')
            lineArr = []
            for i in range(21):
                lineArr.append(float(currLine[i]))
            if int(classifyVector(array(lineArr), trainWeights)) != int(
                    currLine[21]):
                errorCount += 1
        errorRate = (float(errorCount) / numTestVec)
        print("the error rate of this test is: %f" % errorRate)
        return errorRate
        
    
    # 调用 colicTest() 10次并求结果的平均值
    def multiTest():
        numTests = 10
        errorSum = 0.0
        for k in range(numTests):
            errorSum += colicTest()
        print("after %d iterations the average error rate is: %f" % (numTests, errorSum / float(numTests)))

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 12]  
*使用算法: 实现一个简单的命令行程序来收集马的症状并输出预测结果并非难事，这可以作为留给大家的一道习题*

### 4、Logistic回归 小结 ###

> Logistic回归的目的是寻找一个非线性函数Sigmoid的佳拟合参数，求解过程可以由优化 算法来完成。在优化算法中，常用的就是梯度上升算法，而梯度上升算法又可以简化为随机 梯度上升算法。  
> 随机梯度上升算法与梯度上升算法的效果相当，但占用更少的计算资源。此外，随机梯度上 升是一个在线算法，它可以在新数据到来时就完成参数更新，而不需要重新读取整个数据集来进 行批处理运算。  
> 机器学习的一个重要问题就是如何处理缺失数据。这个问题没有标准答案，取决于实际应用 中的需求。现有一些解决方案，每种方案都各有优缺点。

# 资料来源 #

https://github.com/apachecn/AiLearning  
https://github.com/apachecn/AiLearning/tree/master/docs/ml  
机器学习实战（作者: Peter Harrington 出版社: 人民邮电出版社原作名: Machine Learning in Action译者: 李锐 / 李鹏 / 曲亚东 / 王斌 ）  
矫正错误提供的帮助：https://blog.csdn.net/weixin\_43744799/article/details/86087542

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70]: https://img-blog.csdnimg.cn/20191128212137909.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5,size_16,color_FFFFFF,t_70
[20191128212448540.png]: https://img-blog.csdnimg.cn/20191128212448540.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 1]: https://img-blog.csdnimg.cn/20191128212502582.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5,size_16,color_FFFFFF,t_70
[20191128212615487.png]: https://img-blog.csdnimg.cn/20191128212615487.png
[20191128212720251.png]: https://img-blog.csdnimg.cn/20191128212720251.png
[https_blog.csdn.net_qq_45556599_article_details_103315806]: https://blog.csdn.net/qq_45556599/article/details/103315806
[20191128221233424.png]: https://img-blog.csdnimg.cn/20191128221233424.png
[20191128221454740.png]: https://img-blog.csdnimg.cn/20191128221454740.png
[20191128221609785.png]: https://img-blog.csdnimg.cn/20191128221609785.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 2]: https://img-blog.csdnimg.cn/20191128222147554.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5,size_16,color_FFFFFF,t_70
[20191129164311554.png]: https://img-blog.csdnimg.cn/20191129164311554.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 3]: https://img-blog.csdnimg.cn/20191129164702390.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 4]: https://img-blog.csdnimg.cn/20191129164510185.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5,size_16,color_FFFFFF,t_70
[Sigmoid]: https://img-blog.csdnimg.cn/2019112917310339.png
[Sigmoid 1]: https://img-blog.csdnimg.cn/20191129173139823.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5,size_16,color_FFFFFF,t_70
[Sigmoid 2]: https://img-blog.csdnimg.cn/20191129173339169.png
[20191129173513504.png]: https://img-blog.csdnimg.cn/20191129173513504.png
[f_x_ y_x]: https://img-blog.csdnimg.cn/20191129173551549.png
[f_x_ y_y]: https://img-blog.csdnimg.cn/20191129173610556.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 5]: https://img-blog.csdnimg.cn/20191129173636622.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5,size_16,color_FFFFFF,t_70
[20191129173715281.png]: https://img-blog.csdnimg.cn/20191129173715281.png
[20191129173827306.png]: https://img-blog.csdnimg.cn/20191129173827306.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 6]: https://img-blog.csdnimg.cn/20191129173856500.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5,size_16,color_FFFFFF,t_70
[https_blog.csdn.net_qq_45556599_article_details_103315793]: https://blog.csdn.net/qq_45556599/article/details/103315793
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 7]: https://img-blog.csdnimg.cn/20191129174455365.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 8]: https://img-blog.csdnimg.cn/20191129180415323.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 9]: https://img-blog.csdnimg.cn/20191129180534981.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 10]: https://img-blog.csdnimg.cn/20191129180117204.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 11]: https://img-blog.csdnimg.cn/20191129180617945.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5_size_16_color_FFFFFF_t_70 12]: https://img-blog.csdnimg.cn/20191129181336644.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQ1NTU2NTk5,size_16,color_FFFFFF,t_70