机器学习十大经典算法——环境搭建与算例分析

机器学习十大经典算法——环境搭建与算例分析

根据对机器学习相关算法资料的收集情况，我们总结出了十四个经典常用的机器学习算法，分别在Windows环境下和Linux环境下搭建了开发环境，对每一个算法进行了原理分析与代码实现。同时，对于每一个算法我们也给出了相应的算例分析，希望达到在学习的同时也达到进一步运用的目的。十个算法如下：

• 1、K近邻算法
• 2、决策树
• 3、随机森林
• 4、逻辑回归
• 5、支持向量机（SVM）
• 6、朴素贝叶斯
• 7、K均值算法（聚类算法）
• 8、集成算法(Adaboost)
• 9、线性回归（最小二乘法）
• 10、神经网络

环境搭建

考虑到python在科学计算与数据分析方面的方便与快捷，我们首选了python语言作为本文所有算法实现的第一语言，那么我们必须搭建起一个合适的开发环境。
众所周知，Anaconda集成了最新版本的python环境，同时也包含了conda、python等180多个科学包及其依赖项，因此它几乎是每一个用python做数据分析的人的首选开发环境。毫无疑问，我们也选了Anaconda做为我们的开发环境。

Anaconda的安装

Anaconda是开源的软件，它的安装也极其简单，参考文档,下面给出了安装的步骤：

下载

Windows系统下：

官网下载

选择Windows系统，推荐下载集成了python3版本Anaconda, 毕竟未来python2是要停止维护的

Linux系统下：

官网下载

选择Linux系统，推荐下载集成了python3版本Anaconda, 毕竟未来python2是要停止维护的

安装

按照安装程序提示一步步安装就好了, 安装完成之后会多几个应用：

• Anaconda Navigtor ：用于管理工具包和环境的图形用户界面，后续涉及的众多管理命令也可以在 Navigator 中手工实现。
• Jupyter notebook ：基于web的交互式计算环境，可以编辑易于人们阅读的文档，用于展示数据分析的过程。
• qtconsole ：一个可执行 IPython 的仿终端图形界面程序，相比 Python Shell 界面，qtconsole 可以直接显示代码生成的图形，实现多行代码输入执行，以及内置许多有用的功能和函数。
• spyder ：一个使用Python语言、跨平台的、科学运算集成开发环境。

配置环境变量

Windows系统

• windows的话需要去： 控制面板\系统和安全\系统\高级系统设置\环境变量\用户变量\PATH 中添加 anaconda的安装目录的Scripts文件夹, 比如我的路径是C:\Anaconda\Scripts, 看个人安装路径不同需要自己调整.

• 之后就可以打开命令行(最好用管理员模式打开) 输入 ‘conda --version’ ，如果输出conda 4.5.6之类的就说明环境变量设置成功了.

• 为了避免可能发生的错误, 我们在命令行输入’conda upgrade --all’ , 先把所有工具包进行升级.

• 环境搭建完成

Linux系统
参考文档

• 1.下载installer

  推荐清华镜像下载地址

  我下载的版本是Anaconda3-5.3.0-Linux-x86_64.sh

• 2.进入下载好文件的文件夹，在命令行中输入：

  bash Anaconda3-5.3.0-Linux-x86_64.sh
  

• 3.接下来根据提示进行安装，一般默认为default，之后还可以修改。

  在安装过程中，如果之前已经安装了pyhton，numpy等依赖包，它会给你一个提示，选择将Anaconda路径添加上就行，不需要将以前的删掉。

• 4.重新打开命令行，并进行测试

  conda list
  

  如果显示很多python包，说明安装成功了，接下来就可以正常使用了

一、K-近邻法(k-nearest neighbor, k-NN)

1.1 算法简述

k近邻法(k-nearest neighbor, k-NN)是1967年由Cover T和Hart P提出的一种基本分类与回归方法。它的工作原理是：存在一个样本数据集合，也称作为训练样本集，并且样本集中每个数据都存在标签，即我们知道样本集中每一个数据与所属分类的对应关系。输入没有标签的新数据后，将新的数据的每个特征与样本集中数据对应的特征进行比较，然后算法提取样本最相似数据(最近邻)的分类标签。一般来说，我们只选择样本数据集中前k个最相似的数据，这就是k-近邻算法中k的出处，通常k是不大于20的整数。最后，选择k个最相似数据中出现次数最多的分类，作为新数据的分类。

1.2 算法的数学模型

针对于不同维度的距离求解，我们可以用到欧氏距离公式，公式如下：

1.3 算法的实现步骤

• 1、计算已知类别数据集中的点与当前点之间的距离；
• 2、按照距离递增次序排序；
• 3、选取与当前点距离最小的k个点；
• 4、确定前k个点所在类别的出现频率；
• 5、返回前k个点所出现频率最高的类别作为当前点的预测分类。

1.4 算法的运用实例

1.4.1 算法实现的一般流程

• 收集数据：可以使用爬虫进行数据的收集，也可以使用第三方提供的免费或收费的数据。一般来讲，数据放在txt文本文件中，按照一定的格式进行存储，便于解析及处理。
• 准备数据：使用Python解析、预处理数据。
• 分析数据：可以使用很多方法对数据进行分析，例如使用Matplotlib将数据可视化。
• 测试算法：计算错误率。
• 使用算法：错误率在可接受范围内，就可以运行k-近邻算法进行分类。

1.4.2 问题描述

海伦女士一直使用在线约会网站寻找适合自己的约会对象。尽管约会网站会推荐不同的任选，但她并不是喜欢每一个人。经过一番总结，她发现自己交往过的人可以进行如下分类：

• 不喜欢的人
• 魅力一般的人
• 极具魅力的人

海伦收集约会数据已经有了一段时间，她把这些数据存放在文本文件data/datingTestSet.txt中，每个样本数据占据一行，总共有1000行。

海伦收集的样本数据主要包含以下3种特征：

• 每年获得的飞行常客里程数
• 玩视频游戏所消耗时间百分比
• 每周消费的冰淇淋公升数

请利用数据训练模型，用户给一组相关数据，能够较为准确的预测该数据属于哪一类？

1.4.3 数据准备：数据解析

在将上述特征数据输入到分类器前，必须将待处理的数据的格式改变为分类器可以接收的格式，即要将数据分类两部分：特征矩阵和对应的分类标签向量。简而言之，这一步主要是将目标数据加载到内存中，已形成我们想要的数据结构，以便于进一步分析。

在kNN.py文件中创建名为file2matrix的函数，以此来处理输入格式问题。 将datingTestSet.txt放到与kNN.py相同目录下的data文件夹中，编写代码如下：

import numpy as np
"""
函数说明:打开并解析文件，对数据进行分类：1代表不喜欢,2代表魅力一般,3代表极具魅力
Parameters:
    filename - 文件名
Returns:
    returnMat - 特征矩阵
    classLabelVector - 分类Label向量
"""
def file2matrix(filename):
    #打开文件
    fr = open(filename)
    #读取文件所有内容
    arrayOLines = fr.readlines()
    #得到文件行数
    numberOfLines = len(arrayOLines)
    #返回的NumPy矩阵,解析完成的数据:numberOfLines行,3列
    returnMat = np.zeros((numberOfLines,3))
    #返回的分类标签向量
    classLabelVector = []
    #行的索引值
    index = 0
    for line in arrayOLines:
        #s.strip(rm)，当rm空时,默认删除空白符(包括'\n','\r','\t',' ')
        line = line.strip()
        #使用s.split(str="",num=string,cout(str))将字符串根据'\t'分隔符进行切片。
        listFromLine = line.split('\t')
        #将数据前三列提取出来,存放到returnMat的NumPy矩阵中,也就是特征矩阵
        returnMat[index,:] = listFromLine[0:3]
        #根据文本中标记的喜欢的程度进行分类,1代表不喜欢,2代表魅力一般,3代表极具魅力
        if listFromLine[-1] == 'didntLike':
            classLabelVector.append(1)
        elif listFromLine[-1] == 'smallDoses':
            classLabelVector.append(2)
        elif listFromLine[-1] == 'largeDoses':
            classLabelVector.append(3)
        index += 1
    return returnMat, classLabelVector

我们已经顺利导入数据，并对数据进行解析，格式化为分类器需要的数据格式。接着我们需要了解数据的真正含义。可以通过友好、直观的图形化的方式观察数据。

1.4.4 分析数据：数据可视化

在kNN.py文件中编写名为showdatas的函数，用来将数据可视化。编写代码如下：

from matplotlib.font_manager import FontProperties
import matplotlib.lines as mlines
import matplotlib.pyplot as plt
import numpy as np

"""
函数说明:可视化数据

Parameters:
    datingDataMat - 特征矩阵
    datingLabels - 分类Label
Returns:
    无
"""
def showdatas(datingDataMat, datingLabels):
    #设置汉字格式
    font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
    #将fig画布分隔成1行1列,不共享x轴和y轴,fig画布的大小为(13,8)
    #当nrow=2,nclos=2时,代表fig画布被分为四个区域,axs[0][0]表示第一行第一个区域
    fig, axs = plt.subplots(nrows=2, ncols=2,sharex=False, sharey=False, figsize=(13,8))

    numberOfLabels = len(datingLabels)
    LabelsColors = []
    for i in datingLabels:
        if i == 1:
            LabelsColors.append('black')
        if i == 2:
            LabelsColors.append('orange')
        if i == 3:
            LabelsColors.append('red')
    #画出散点图,以datingDataMat矩阵的第一(飞行常客例程)、第二列(玩游戏)数据画散点数据,散点大小为15,透明度为0.5
    axs[0][0].scatter(x=datingDataMat[:,0], y=datingDataMat[:,1], color=LabelsColors,s=15, alpha=.5)
    #设置标题,x轴label,y轴label
    axs0_title_text = axs[0][0].set_title(u'每年获得的飞行常客里程数与玩视频游戏所消耗时间占比',FontProperties=font)
    axs0_xlabel_text = axs[0][0].set_xlabel(u'每年获得的飞行常客里程数',FontProperties=font)
    axs0_ylabel_text = axs[0][0].set_ylabel(u'玩视频游戏所消耗时间占',FontProperties=font)
    plt.setp(axs0_title_text, size=9, weight='bold', color='red') 
    plt.setp(axs0_xlabel_text, size=7, weight='bold', color='black') 
    plt.setp(axs0_ylabel_text, size=7, weight='bold', color='black')

    #画出散点图,以datingDataMat矩阵的第一(飞行常客例程)、第三列(冰激凌)数据画散点数据,散点大小为15,透明度为0.5
    axs[0][1].scatter(x=datingDataMat[:,0], y=datingDataMat[:,2], color=LabelsColors,s=15, alpha=.5)
    #设置标题,x轴label,y轴label
    axs1_title_text = axs[0][1].set_title(u'每年获得的飞行常客里程数与每周消费的冰激淋公升数',FontProperties=font)
    axs1_xlabel_text = axs[0][1].set_xlabel(u'每年获得的飞行常客里程数',FontProperties=font)
    axs1_ylabel_text = axs[0][1].set_ylabel(u'每周消费的冰激淋公升数',FontProperties=font)
    plt.setp(axs1_title_text, size=9, weight='bold', color='red') 
    plt.setp(axs1_xlabel_text, size=7, weight='bold', color='black') 
    plt.setp(axs1_ylabel_text, size=7, weight='bold', color='black')

    #画出散点图,以datingDataMat矩阵的第二(玩游戏)、第三列(冰激凌)数据画散点数据,散点大小为15,透明度为0.5
    axs[1][0].scatter(x=datingDataMat[:,1], y=datingDataMat[:,2], color=LabelsColors,s=15, alpha=.5)
    #设置标题,x轴label,y轴label
    axs2_title_text = axs[1][0].set_title(u'玩视频游戏所消耗时间占比与每周消费的冰激淋公升数',FontProperties=font)
    axs2_xlabel_text = axs[1][0].set_xlabel(u'玩视频游戏所消耗时间占比',FontProperties=font)
    axs2_ylabel_text = axs[1][0].set_ylabel(u'每周消费的冰激淋公升数',FontProperties=font)
    plt.setp(axs2_title_text, size=9, weight='bold', color='red') 
    plt.setp(axs2_xlabel_text, size=7, weight='bold', color='black') 
    plt.setp(axs2_ylabel_text, size=7, weight='bold', color='black')
    #设置图例
    didntLike = mlines.Line2D([], [], color='black', marker='.',
                      markersize=6, label='didntLike')
    smallDoses = mlines.Line2D([], [], color='orange', marker='.',
                      markersize=6, label='smallDoses')
    largeDoses = mlines.Line2D([], [], color='red', marker='.',
                      markersize=6, label='largeDoses')
    #添加图例
    axs[0][0].legend(handles=[didntLike,smallDoses,largeDoses])
    axs[0][1].legend(handles=[didntLike,smallDoses,largeDoses])
    axs[1][0].legend(handles=[didntLike,smallDoses,largeDoses])
    #显示图片
    plt.show()

可视化结果如下图所示：

通过可视化数据可很直观的表达出数据的规律，观察图中的三幅图可知，海伦喜欢娱乐和工作能够平衡的男士。为什么这么说呢？每年获得的飞行常客里程数表明，海伦喜欢工作能力强的男人，但是不能经常坐飞机，疲于奔波，满世界飞。同时，这个男人也要玩视频游戏，并且占一定时间比例。

1.4.5 准备数据：数据归一化

如果直接用距离计算公式，我们会发祥每年获取的飞行常客里程数对于计算结果的影响将远远大于其他两个特征-玩视频游戏所耗时间占比和每周消费冰淇淋公斤数的影响。而产生这种现象的唯一原因，仅仅是因为飞行常客里程数远大于其他特征值。但海伦认为这三种特征是同等重要的，因此作为三个等权重的特征之一，飞行常客里程数并不应该如此严重地影响到计算结果。
所以，我们通常采用的方法是将数值归一化，如将取值范围处理为０到１或者-１到１之间。下面的公式可以将任意取值范围的特征值转化为０到１区间内的值：

newValue = (oldValue - min) / (max - min)

其中min和max分别是数据集中的最小特征值和最大特征值.在kNN.py文件中编写名为autoNorm的函数，用该函数自动将数据归一化。代码如下：

"""
函数说明:对数据进行归一化

Parameters:
    dataSet - 特征矩阵
Returns:
    normDataSet - 归一化后的特征矩阵
    ranges - 数据范围
    minVals - 数据最小值

"""
def autoNorm(dataSet):
    #获得数据的最小值
    minVals = dataSet.min(0)
    maxVals = dataSet.max(0)
    #最大值和最小值的范围
    ranges = maxVals - minVals
    #shape(dataSet)返回dataSet的矩阵行列数
    normDataSet = np.zeros(np.shape(dataSet))
    #返回dataSet的行数
    m = dataSet.shape[0]
    #原始值减去最小值
    normDataSet = dataSet - np.tile(minVals, (m, 1))
    #除以最大和最小值的差,得到归一化数据
    normDataSet = normDataSet / np.tile(ranges, (m, 1))
    #返回归一化数据结果,数据范围,最小值
    return normDataSet, ranges, minVals

通过这一步我们已经顺利将数据归一化了，并且求出了数据的取值范围和数据的最小值，这两个值是在分类的时候需要用到的，直接先求解出来，也算是对数据预处理了。

1.4.6 测试算法：验证分类器

机器学习算法一个很重要的工作就是评估算法的正确率，通常我们只提供已有数据的90%作为训练样本来训练分类器，而使用其余的10%数据去测试分类器，检测分类器的正确率。需要注意的是，10%的测试数据应该是随机选择的，由于海伦提供的数据并没有按照特定目的来排序，所以我们可以随意选择10%数据而不影响其随机性。
为了测试分类器效果，在kNN.py文件中创建分类器函数classify0函数与验证函数datingClassTest，编写代码如下：

"""
函数说明:kNN算法,分类器

Parameters:
    inX - 用于分类的数据(测试集)
    dataSet - 用于训练的数据(训练集)
    labes - 分类标签
    k - kNN算法参数,选择距离最小的k个点
Returns:
    sortedClassCount[0][0] - 分类结果


"""
def classify0(inX, dataSet, labels, k):
    #numpy函数shape[0]返回dataSet的行数
    dataSetSize = dataSet.shape[0]
    #在列向量方向上重复inX共1次(横向),行向量方向上重复inX共dataSetSize次(纵向)
    diffMat = np.tile(inX, (dataSetSize, 1)) - dataSet
    #二维特征相减后平方
    sqDiffMat = diffMat**2
    #sum()所有元素相加,sum(0)列相加,sum(1)行相加
    sqDistances = sqDiffMat.sum(axis=1)
    #开方,计算出距离
    distances = sqDistances**0.5
    #返回distances中元素从小到大排序后的索引值
    sortedDistIndices = distances.argsort()
    #定一个记录类别次数的字典
    classCount = {}
    for i in range(k):
        #取出前k个元素的类别
        voteIlabel = labels[sortedDistIndices[i]]
        #dict.get(key,default=None),字典的get()方法,返回指定键的值,如果值不在字典中返回默认值。
        #计算类别次数
        classCount[voteIlabel] = classCount.get(voteIlabel,0) + 1
    #python3中用items()替换python2中的iteritems()
    #key=operator.itemgetter(1)根据字典的值进行排序
    #key=operator.itemgetter(0)根据字典的键进行排序
    #reverse降序排序字典
    sortedClassCount = sorted(classCount.items(),key=operator.itemgetter(1),reverse=True)
    #返回次数最多的类别,即所要分类的类别
    return sortedClassCount[0][0]


"""
函数说明:分类器测试函数

Parameters:
    无
Returns:
    normDataSet - 归一化后的特征矩阵
    ranges - 数据范围
    minVals - 数据最小值

"""
def datingClassTest():
    #打开的文件名
    filename = "data/datingTestSet.txt"
    #将返回的特征矩阵和分类向量分别存储到datingDataMat和datingLabels中
    datingDataMat, datingLabels = file2matrix(filename)
    #取所有数据的百分之十
    hoRatio = 0.10
    #数据归一化,返回归一化后的矩阵,数据范围,数据最小值
    normMat, ranges, minVals = autoNorm(datingDataMat)
    #获得normMat的行数
    m = normMat.shape[0]
    #百分之十的测试数据的个数
    numTestVecs = int(m * hoRatio)
    #分类错误计数
    errorCount = 0.0

    for i in range(numTestVecs):
        #前numTestVecs个数据作为测试集,后m-numTestVecs个数据作为训练集
        classifierResult = classify0(normMat[i,:], normMat[numTestVecs:m,:],
            datingLabels[numTestVecs:m], 4)
        print("分类结果:%d\t真实类别:%d" % (classifierResult, datingLabels[i]))
        if classifierResult != datingLabels[i]:
            errorCount += 1.0
    print("错误率:%f%%" %(errorCount/float(numTestVecs)*100))

运行结果为：

从测试结果可以看出，错误率是4%，这是一个相当不错的结果，我们可以改变函数datingClassTest内变量hoRatio和分类器k的值，检测错误率是否随着变量值的变化而增加。依赖于分类算法、数据集和程序设置，分类器的输出结果可能有很大的不同。

1.4.7 使用算法：构建完整可用系统

我们可以给海伦一个小段程序，通过该程序海伦会在约会网站上找到某个人并输入他的信息。程序会给出她对男方喜欢程度的预测值。
在kNN.py文件中创建函数classifyPerson，整个程序预测的代码如下：

# -*- coding: UTF-8 -*-

import numpy as np
import operator

"""
函数说明:打开并解析文件，对数据进行分类：1代表不喜欢,2代表魅力一般,3代表极具魅力

Parameters:
    filename - 文件名
Returns:
    returnMat - 特征矩阵
    classLabelVector - 分类Label向量

"""
def file2matrix(filename):
    #打开文件
    fr = open(filename)
    #读取文件所有内容
    arrayOLines = fr.readlines()
    #得到文件行数
    numberOfLines = len(arrayOLines)
    #返回的NumPy矩阵,解析完成的数据:numberOfLines行,3列
    returnMat = np.zeros((numberOfLines,3))
    #返回的分类标签向量
    classLabelVector = []
    #行的索引值
    index = 0
    for line in arrayOLines:
        #s.strip(rm)，当rm空时,默认删除空白符(包括'\n','\r','\t',' ')
        line = line.strip()
        #使用s.split(str="",num=string,cout(str))将字符串根据'\t'分隔符进行切片。
        listFromLine = line.split('\t')
        #将数据前三列提取出来,存放到returnMat的NumPy矩阵中,也就是特征矩阵
        returnMat[index,:] = listFromLine[0:3]
        #根据文本中标记的喜欢的程度进行分类,1代表不喜欢,2代表魅力一般,3代表极具魅力
        if listFromLine[-1] == 'didntLike':
            classLabelVector.append(1)
        elif listFromLine[-1] == 'smallDoses':
            classLabelVector.append(2)
        elif listFromLine[-1] == 'largeDoses':
            classLabelVector.append(3)
        index += 1
    return returnMat, classLabelVector
"""
函数说明:通过输入一个人的三维特征,进行分类输出

Parameters:
    无
Returns:
    无
"""
def classifyPerson():
    #输出结果
    resultList = ['讨厌','有些喜欢','非常喜欢']
    #三维特征用户输入
    ffMiles = float(input("每年获得的飞行常客里程数:"))
    precentTats = float(input("玩视频游戏所耗时间百分比:"))
    iceCream = float(input("每周消费的冰激淋公升数:"))
    #打开的文件名
    filename = "data/datingTestSet.txt"
    #打开并处理数据
    datingDataMat, datingLabels = file2matrix(filename)
    #训练集归一化
    normMat, ranges, minVals = autoNorm(datingDataMat)
    #生成NumPy数组,测试集
    inArr = np.array([ ffMiles,precentTats,iceCream])
    #测试集归一化
    norminArr = (inArr - minVals) / ranges
    #返回分类结果
    classifierResult = classify0(norminArr, normMat, datingLabels, 3)
    #打印结果
    print("你可能%s这个人" % (resultList[classifierResult-1]))

"""
函数说明:main函数

Parameters:
    无
Returns:
    无

"""
if __name__ == '__main__':
    classifyPerson()

运行程序，并输入数据(12,44000,0.5)，预测结果是"你可能有些喜欢这个人"，也就是这个人魅力一般。一共有三个档次：讨厌、有些喜欢、非常喜欢，对应着不喜欢的人、魅力一般的人、极具魅力的人。结果如图所示：

二、K均值算法

2.1 算法简述

K-means算法是很典型的基于距离的聚类算法，采用距离作为相似性的评价指标，即认为两个对象的距离越近，其相似度就越大。该算法认为簇是由距离靠近的对象组成的，因此把得到紧凑且独立的簇作为最终目标。

选择K个点作为初始质心
repeat
	将每个点指派到最近的质心，形成K个簇
	重新计算每个簇的质心
until 簇不发生变化或达到最大迭代次数

2.2 算法的数学模型

计算相似度的公式（这里采用欧氏距离）：

Kmeans聚类算法的优化目标函数：

μi是第k个类的重心位置。成本函数是各个类畸变程度(distortions)之和。

2.3 算法的流程

输入：聚类个数k，数据集X。

• 首先从n个数据选择k个初始中心点，例如c[0]=X[0],…,c[k-1]=X[k-1]作为初始聚类中心；

• 对于所剩下其它对象，则根据它们与这些聚类中心的相似度（距离），分别将它们分配给与其最相似的（聚类中心所代表的）聚类；分别与c[0]…c[k-1]比较，假定与c[i]差值最少，就标记i；

• 然后再计算每个所获新聚类的聚类中心（该聚类中所有对象的均值）；

• 不断重复这一过程直到直到所有c[i]值的变化小于给定阈值或者达到最大迭代次数。

输出：满足方差最小标准的k个聚类。

注意：一般都采用均方差作为标准测度函数。k个聚类具有以下特点：各聚类本身尽可能的紧凑，而各聚类之间尽可能的分开。

2.4 算法的运用实例

2.4.1 算法实现的重点：

K值的确定
初始质心的选取
距离的度量
聚类效果评估

2.4.2 问题描述

对于给定的数据集，要求把其分为4类；

数据的读入

def DataInput(file_path):
    dataSet = []
    with open(file_path) as fileIn:
        for line in fileIn.readlines():
            lineArr = line.strip().split('\t')
            dataSet.append([float(lineArr[0]), float(lineArr[1])])
    return dataSet

K值的确定

这里只讨论K值已知的情况，根据目标分类的个数可知:

k = 4

初始质心的选取

这里选用随机选取初始质心的方法

# 随机选取初始质心
def initCentroids(dataSet, k):
	numSamples, dim = dataSet.shape
	centroids = zeros((k, dim))
	for i in range(k):
	    # 从一个均匀分布[low,high)中随机采样,选取初始质心
		index = int(random.uniform(0, numSamples))
		centroids[i, :] = dataSet[index, :]
	return centroids

距离的度量

采用欧氏距离

# 计算欧氏距离
def euclDistance(vector1, vector2):
	return sqrt(sum(power(vector2 - vector1, 2)))

聚类过程

# k-means cluster
def kmeans(dataSet, k):
	numSamples = dataSet.shape[0]
	# first column stores which cluster this sample belongs to,
	# second column stores the error between this sample and its centroid
	clusterAssment = mat(zeros((numSamples, 2)))
	clusterChanged = True
 
	## step 1: init centroids
	centroids = initCentroids(dataSet, k)
 
	while clusterChanged:
		clusterChanged = False
		## for each sample
		for i in range(numSamples):
			minDist  = 100000.0
			minIndex = 0
			## for each centroid
			## step 2: find the centroid who is closest
			for j in range(k):
				distance = euclDistance(centroids[j, :], dataSet[i, :])
				if distance < minDist:
					minDist  = distance
					minIndex = j
			
			## step 3: update its cluster
			if clusterAssment[i, 0] != minIndex:
				clusterChanged = True
				clusterAssment[i, :] = minIndex, minDist**2
 
		## step 4: update centroids
		for j in range(k):
			pointsInCluster = dataSet[nonzero(clusterAssment[:, 0].A == j)[0]]
			centroids[j, :] = mean(pointsInCluster, axis = 0)
 
	print('Congratulations, cluster complete!')
	return centroids, clusterAssment

聚类结果的可视化

import matplotlib.pyplot as plt
def showCluster(dataSet, k, centroids, clusterAssment):
	numSamples, dim = dataSet.shape
	if dim != 2:
		print("Sorry! I can not draw because the dimension of your data is not 2!")
		return 1
 
	mark = ['or', 'ob', 'og', 'ok', '^r', '+r', 'sr', 'dr', '<r', 'pr']
	if k > len(mark):
		print("Sorry! Your k is too large! please contact Zouxy")
		return 1
 
	# draw all samples
	for i in range(numSamples):
		markIndex = int(clusterAssment[i, 0])
		plt.plot(dataSet[i, 0], dataSet[i, 1], mark[markIndex])
 
	mark = ['Dr', 'Db', 'Dg', 'Dk', '^b', '+b', 'sb', 'db', '<b', 'pb']
	# draw the centroids
	for i in range(k):
		plt.plot(centroids[i, 0], centroids[i, 1], mark[i], markersize = 12)
	plt.show()
showCluster(dataSet, k, centroids, clusterAssment)

运行测试

from numpy import *
import time
import matplotlib.pyplot as plt

···

if __name__ == "__main__":
    # 数据读入
    file_path = r'data/k-means.txt'
    dataSet = DataInput(file_path)
    
    # 设置k值
    k = 4
    dataSet = mat(dataSet)
    
    # 执行聚类
    centroids, clusterAssment = kmeans(dataSet, k)
    
    # 聚类结果可视化
    showCluster(dataSet, k, centroids, clusterAssment)

运行结果为：

K-Means优化算法

为了克服K-Means算法收敛于局部最小值的问题，提出了一种二分K-均值(bisecting K-means)

• K-Means优化算法的伪代码如下：

  将所有的点看成是一个簇
  当簇小于数目k时
  对于每一个簇
      计算总误差
      在给定的簇上进行K-均值聚类,k值为2
      计算将该簇划分成两个簇后总误差
  选择是的误差最小的那个簇进行划分
  
  

三、决策树算法

3.1 算法简述

决策树(decision tree)是一种基本的分类与回归方法。

3.2 算法的数学模型

香农熵

熵定义为信息的期望值。在信息论与概率统计中，熵是表示随机变量不确定性的度量。如果待分类的事物可能划分在多个分类之中，则符号xi的信息定义为 ：

其中p(xi)是选择该分类的概率。上述式中的对数以2为底，也可以e为底(自然对数)。

通过上式，我们可以得到所有类别的信息。为了计算熵，我们需要计算所有类别所有可能值包含的信息期望值(数学期望)，通过下面的公式得到：

期中n是分类的数目。熵越大，随机变量的不确定性就越大。

经验熵

当熵中的概率由数据估计(特别是最大似然估计)得到时，所对应的熵称为经验熵(empirical entropy)。
我们定义样本数据表中的数据为训练数据集D，则训练数据集D的经验熵为H(D)，|D|表示其样本容量，及样本个数。设有K个类Ck, = 1,2,3,…,K,|Ck|为属于类Ck的样本个数，因此经验熵公式就可以写为：

条件熵

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

这里，

同理，当条件熵中的概率由数据估计(特别是极大似然估计)得到时，所对应的条件熵成为条件经验熵(empirical conditional entropy)。

信息增益

明确了条件熵和经验条件熵的概念。接下来，让我们说说信息增益。前面也提到了，信息增益是相对于特征而言的。所以，特征A对训练数据集D的信息增益g(D,A)，定义为集合D的经验熵H(D)与特征A给定条件下D的经验条件熵H(D|A)之差，即：

一般地，熵H(D)与条件熵H(D|A)之差成为互信息(mutual information)。决策树学习中的信息增益等价于训练数据集中类与特征的互信息。
设特征A有n个不同的取值{a1,a2,···,an}，根据特征A的取值将D划分为n个子集{D1,D2，···,Dn}，|Di|为Di的样本个数。记子集Di中属于Ck的样本的集合为Dik，即Dik = Di ∩ Ck，|Dik|为Dik的样本个数。于是经验条件熵的公式可以些为：

3.3 算法的实现步骤

• 特征选择
• 决策树的生成
• 决策树的修剪

3.4 算法的运用实例

3.4.2 问题描述

训练data/DecisionTree.txt所给的数据构造决策树，实现由用户提供相关数据，可得出可信的预测。

数据准备

"""
函数说明:读入数据，序列化
 
Parameters:
    
Returns:
    lenses_pd - 序列化后的数据集
"""
import pandas as pd
from sklearn.preprocessing import LabelEncoder

import pydotplus
from sklearn.externals.six import StringIO

def creatDataSet():
    with open('data/DecisionTree.txt', 'r') as fr:                                        #加载文件
        lenses = [inst.strip().split('\t') for inst in fr.readlines()]        #处理文件
    lenses_target = []                                                        #提取每组数据的类别，保存在列表里
    for each in lenses:
        lenses_target.append(each[-1])
 
    lensesLabels = ['age', 'prescript', 'astigmatic', 'tearRate']            #特征标签       
    lenses_list = []                                                        #保存lenses数据的临时列表
    lenses_dict = {}                                                        #保存lenses数据的字典，用于生成pandas
    for each_label in lensesLabels:                                            #提取信息，生成字典
        for each in lenses:
            lenses_list.append(each[lensesLabels.index(each_label)])
        lenses_dict[each_label] = lenses_list
        lenses_list = []
    lenses_pd = pd.DataFrame(lenses_dict)                                    #生成pandas.DataFrame
    print(lenses_pd)
    le = LabelEncoder()                                                        #创建LabelEncoder()对象，用于序列化            
    for col in lenses_pd.columns:                                            #为每一列序列化
        lenses_pd[col] = le.fit_transform(lenses_pd[col])
    print(lenses_pd)

特征选择

"""
函数说明:计算给定数据集的经验熵(香农熵)
 
Parameters:
    dataSet - 数据集
Returns:
    shannonEnt - 经验熵(香农熵)
"""
def calcShannonEnt(dataSet):
    numEntires = len(dataSet)                        #返回数据集的行数
    labelCounts = {}                                #保存每个标签(Label)出现次数的字典
    for featVec in dataSet:                            #对每组特征向量进行统计
        currentLabel = featVec[-1]                    #提取标签(Label)信息
        if currentLabel not in labelCounts.keys():    #如果标签(Label)没有放入统计次数的字典,添加进去
            labelCounts[currentLabel] = 0
        labelCounts[currentLabel] += 1                #Label计数
    shannonEnt = 0.0                                #经验熵(香农熵)
    for key in labelCounts:                            #计算香农熵
        prob = float(labelCounts[key]) / numEntires    #选择该标签(Label)的概率
        shannonEnt -= prob * log(prob, 2)            #利用公式计算
    return shannonEnt                                #返回经验熵(香农熵)

构建决策树

def createTree(dataSet, labels, featLabels):
    classList = [example[-1] for example in dataSet]            #取分类标签(是否放贷:yes or no)
    if classList.count(classList[0]) == len(classList):            #如果类别完全相同则停止继续划分
        return classList[0]
    if len(dataSet[0]) == 1:                                    #遍历完所有特征时返回出现次数最多的类标签
        return majorityCnt(classList)
    bestFeat = chooseBestFeatureToSplit(dataSet)                #选择最优特征
    bestFeatLabel = labels[bestFeat]                            #最优特征的标签
    featLabels.append(bestFeatLabel)
    myTree = {bestFeatLabel:{}}                                    #根据最优特征的标签生成树
    del(labels[bestFeat])                                        #删除已经使用特征标签
    featValues = [example[bestFeat] for example in dataSet]        #得到训练集中所有最优特征的属性值
    uniqueVals = set(featValues)                                #去掉重复的属性值
    for value in uniqueVals:                                    #遍历特征，创建决策树。                       
        myTree[bestFeatLabel][value] = createTree(splitDataSet(dataSet, bestFeat, value), labels, featLabels)
    return myTree

学习结果以及决策树可视化

# -*- coding: UTF-8 -*-
from math import log
import operator
import pandas as pd
from sklearn.preprocessing import LabelEncoder
from matplotlib.font_manager import FontProperties
import matplotlib.pyplot as plt
 
"""
函数说明:计算给定数据集的经验熵(香农熵)
 
Parameters:
    dataSet - 数据集
Returns:
    shannonEnt - 经验熵(香农熵)
"""
def calcShannonEnt(dataSet):
    numEntires = len(dataSet)                        #返回数据集的行数
    labelCounts = {}                                #保存每个标签(Label)出现次数的字典
    for featVec in dataSet:                            #对每组特征向量进行统计
        currentLabel = featVec[-1]                    #提取标签(Label)信息
        if currentLabel not in labelCounts.keys():    #如果标签(Label)没有放入统计次数的字典,添加进去
            labelCounts[currentLabel] = 0
        labelCounts[currentLabel] += 1                #Label计数
    shannonEnt = 0.0                                #经验熵(香农熵)
    for key in labelCounts:                            #计算香农熵
        prob = float(labelCounts[key]) / numEntires    #选择该标签(Label)的概率
        shannonEnt -= prob * log(prob, 2)            #利用公式计算
    return shannonEnt                                #返回经验熵(香农熵)
 
"""
函数说明:创建测试数据集
 
Parameters:
    无
Returns:
    dataSet - 数据集
    labels - 特征标签
"""

def createDataSet():
    with open('data/DecisionTree.txt', 'r') as fr:                                        #加载文件
        lenses = [inst.strip().split('\t') for inst in fr.readlines()]        #处理文件
    lenses_target = []                                                        #提取每组数据的类别，保存在列表里
    for each in lenses:
        lenses_target.append(each[-1])
 
    lensesLabels = ['age', 'prescript', 'astigmatic', 'tearRate','class']            #特征标签       
    lenses_list = []                                                        #保存lenses数据的临时列表
    lenses_dict = {}                                                        #保存lenses数据的字典，用于生成pandas
    for each_label in lensesLabels:                                            #提取信息，生成字典
        for each in lenses:
            lenses_list.append(each[lensesLabels.index(each_label)])
        lenses_dict[each_label] = lenses_list
        lenses_list = []
    dataSet = pd.DataFrame(lenses_dict)                                    #生成pandas.DataFrame
    # print(dataSet)
    le = LabelEncoder()                                                        #创建LabelEncoder()对象，用于序列化            
    for col in dataSet.columns:                                            #为每一列序列化
        dataSet[col] = le.fit_transform(dataSet[col])
    labels = ['年龄', '症状', '是否散光', '眼泪数量']
    dataSet = dataSet.values.tolist()
    return dataSet, labels                            #返回数据集和分类属性
    

"""
函数说明:按照给定特征划分数据集
 
Parameters:
    dataSet - 待划分的数据集
    axis - 划分数据集的特征
    value - 需要返回的特征的值
Returns:
    无
"""
def splitDataSet(dataSet, axis, value):       
    retDataSet = []                                        #创建返回的数据集列表
    for featVec in dataSet:                             #遍历数据集
        if featVec[axis] == value:
            reducedFeatVec = featVec[:axis]                #去掉axis特征
            reducedFeatVec.extend(featVec[axis+1:])     #将符合条件的添加到返回的数据集
            retDataSet.append(reducedFeatVec)
    return retDataSet                                      #返回划分后的数据集
 
"""
函数说明:选择最优特征
 
Parameters:
    dataSet - 数据集
Returns:
    bestFeature - 信息增益最大的(最优)特征的索引值
"""
def chooseBestFeatureToSplit(dataSet):
    numFeatures = len(dataSet[0]) - 1                    #特征数量
    baseEntropy = calcShannonEnt(dataSet)                 #计算数据集的香农熵
    bestInfoGain = 0.0                                  #信息增益
    bestFeature = -1                                    #最优特征的索引值
    for i in range(numFeatures):                         #遍历所有特征
        #获取dataSet的第i个所有特征
        featList = [example[i] for example in dataSet]
        uniqueVals = set(featList)                         #创建set集合{},元素不可重复
        newEntropy = 0.0                                  #经验条件熵
        for value in uniqueVals:                         #计算信息增益
            subDataSet = splitDataSet(dataSet, i, value)         #subDataSet划分后的子集
            prob = len(subDataSet) / float(len(dataSet))           #计算子集的概率
            newEntropy += prob * calcShannonEnt(subDataSet)     #根据公式计算经验条件熵
        infoGain = baseEntropy - newEntropy                     #信息增益
        # print("第%d个特征的增益为%.3f" % (i, infoGain))            #打印每个特征的信息增益
        if (infoGain > bestInfoGain):                             #计算信息增益
            bestInfoGain = infoGain                             #更新信息增益，找到最大的信息增益
            bestFeature = i                                     #记录信息增益最大的特征的索引值
    return bestFeature                                             #返回信息增益最大的特征的索引值
 
 
"""
函数说明:统计classList中出现此处最多的元素(类标签)
 
Parameters:
    classList - 类标签列表
Returns:
    sortedClassCount[0][0] - 出现此处最多的元素(类标签)
"""
def majorityCnt(classList):
    classCount = {}
    for vote in classList:                                        #统计classList中每个元素出现的次数
        if vote not in classCount.keys():classCount[vote] = 0   
        classCount[vote] += 1
    sortedClassCount = sorted(classCount.items(), key = operator.itemgetter(1), reverse = True)        #根据字典的值降序排序
    return sortedClassCount[0][0]                                #返回classList中出现次数最多的元素
 
"""
函数说明:创建决策树
 
Parameters:
    dataSet - 训练数据集
    labels - 分类属性标签
    featLabels - 存储选择的最优特征标签
Returns:
    myTree - 决策树
"""
def createTree(dataSet, labels, featLabels):
    print("dataset:",dataSet)
    print('******')
    classList = [example[-1] for example in dataSet]            #取分类标签(是否放贷:yes or no)
    print('classList',classList)
    if classList.count(classList[0]) == len(classList):            #如果类别完全相同则停止继续划分
        return classList[0]
    if len(dataSet[0]) == 1 or len(labels) == 0:                                    #遍历完所有特征时返回出现次数最多的类标签
        return majorityCnt(classList)
    bestFeat = chooseBestFeatureToSplit(dataSet)                #选择最优特征
    print(labels)
    print(bestFeat)
    bestFeatLabel = labels[bestFeat]                            #最优特征的标签
    featLabels.append(bestFeatLabel)
    myTree = {bestFeatLabel:{}}                                    #根据最优特征的标签生成树
    print('myTree:',myTree)
    del(labels[bestFeat])                                        #删除已经使用特征标签
    featValues = [example[bestFeat] for example in dataSet]        #得到训练集中所有最优特征的属性值
    uniqueVals = set(featValues)                                #去掉重复的属性值
    for value in uniqueVals:                                    #遍历特征，创建决策树。                       
        myTree[bestFeatLabel][value] = createTree(splitDataSet(dataSet, bestFeat, value), labels, featLabels)
    return myTree

"""
函数说明:获取决策树叶子结点的数目
 
Parameters:
    myTree - 决策树
Returns:
    numLeafs - 决策树的叶子结点的数目
"""
def getNumLeafs(myTree):
    numLeafs = 0                                                #初始化叶子
    firstStr = next(iter(myTree))                                #python3中myTree.keys()返回的是dict_keys,不在是list,所以不能使用myTree.keys()[0]的方法获取结点属性，可以使用list(myTree.keys())[0]
    secondDict = myTree[firstStr]                                #获取下一组字典
    for key in secondDict.keys():
        if type(secondDict[key]).__name__=='dict':                #测试该结点是否为字典，如果不是字典，代表此结点为叶子结点
            numLeafs += getNumLeafs(secondDict[key])
        else:   numLeafs +=1
    return numLeafs
 
"""
函数说明:获取决策树的层数
 
Parameters:
    myTree - 决策树
Returns:
    maxDepth - 决策树的层数
"""
def getTreeDepth(myTree):
    maxDepth = 0                                                #初始化决策树深度
    firstStr = next(iter(myTree))                                #python3中myTree.keys()返回的是dict_keys,不在是list,所以不能使用myTree.keys()[0]的方法获取结点属性，可以使用list(myTree.keys())[0]
    secondDict = myTree[firstStr]                                #获取下一个字典
    for key in secondDict.keys():
        if type(secondDict[key]).__name__=='dict':                #测试该结点是否为字典，如果不是字典，代表此结点为叶子结点
            thisDepth = 1 + getTreeDepth(secondDict[key])
        else:   thisDepth = 1
        if thisDepth > maxDepth: maxDepth = thisDepth            #更新层数
    return maxDepth
 
"""
函数说明:绘制结点
 
Parameters:
    nodeTxt - 结点名
    centerPt - 文本位置
    parentPt - 标注的箭头位置
    nodeType - 结点格式
Returns:
    无
"""
def plotNode(nodeTxt, centerPt, parentPt, nodeType):
    arrow_args = dict(arrowstyle="<-")                                            #定义箭头格式
    font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)        #设置中文字体
    createPlot.ax1.annotate(nodeTxt, xy=parentPt,  xycoords='axes fraction',    #绘制结点
        xytext=centerPt, textcoords='axes fraction',
        va="center", ha="center", bbox=nodeType, arrowprops=arrow_args, FontProperties=font)
 
"""
函数说明:标注有向边属性值
 
Parameters:
    cntrPt、parentPt - 用于计算标注位置
    txtString - 标注的内容
Returns:
    无
"""
def plotMidText(cntrPt, parentPt, txtString):
    xMid = (parentPt[0]-cntrPt[0])/2.0 + cntrPt[0]                                            #计算标注位置                   
    yMid = (parentPt[1]-cntrPt[1])/2.0 + cntrPt[1]
    createPlot.ax1.text(xMid, yMid, txtString, va="center", ha="center", rotation=30)
 
"""
函数说明:绘制决策树
 
Parameters:
    myTree - 决策树(字典)
    parentPt - 标注的内容
    nodeTxt - 结点名
Returns:
    无
"""
def plotTree(myTree, parentPt, nodeTxt):
    decisionNode = dict(boxstyle="sawtooth", fc="0.8")                                        #设置结点格式
    leafNode = dict(boxstyle="round4", fc="0.8")                                            #设置叶结点格式
    numLeafs = getNumLeafs(myTree)                                                          #获取决策树叶结点数目，决定了树的宽度
    depth = getTreeDepth(myTree)                                                            #获取决策树层数
    firstStr = next(iter(myTree))                                                            #下个字典                                                 
    cntrPt = (plotTree.xOff + (1.0 + float(numLeafs))/2.0/plotTree.totalW, plotTree.yOff)    #中心位置
    plotMidText(cntrPt, parentPt, nodeTxt)                                                    #标注有向边属性值
    plotNode(firstStr, cntrPt, parentPt, decisionNode)                                        #绘制结点
    secondDict = myTree[firstStr]                                                            #下一个字典，也就是继续绘制子结点
    plotTree.yOff = plotTree.yOff - 1.0/plotTree.totalD                                        #y偏移
    for key in secondDict.keys():                               
        if type(secondDict[key]).__name__=='dict':                                            #测试该结点是否为字典，如果不是字典，代表此结点为叶子结点
            plotTree(secondDict[key],cntrPt,str(key))                                        #不是叶结点，递归调用继续绘制
        else:                                                                                #如果是叶结点，绘制叶结点，并标注有向边属性值                                             
            plotTree.xOff = plotTree.xOff + 1.0/plotTree.totalW
            plotNode(secondDict[key], (plotTree.xOff, plotTree.yOff), cntrPt, leafNode)
            plotMidText((plotTree.xOff, plotTree.yOff), cntrPt, str(key))
    plotTree.yOff = plotTree.yOff + 1.0/plotTree.totalD

"""
函数说明:创建绘制面板
 
Parameters:
    inTree - 决策树(字典)
Returns:
    无
"""
def createPlot(inTree):
    fig = plt.figure(1, facecolor='white')                                                    #创建fig
    fig.clf()                                                                                #清空fig
    axprops = dict(xticks=[], yticks=[])
    createPlot.ax1 = plt.subplot(111, frameon=False, **axprops)                                #去掉x、y轴
    plotTree.totalW = float(getNumLeafs(inTree))                                            #获取决策树叶结点数目
    plotTree.totalD = float(getTreeDepth(inTree))                                            #获取决策树层数
    plotTree.xOff = -0.5/plotTree.totalW; plotTree.yOff = 1.0;                                #x偏移
    plotTree(inTree, (0.5,1.0), '')                                                            #绘制决策树
    plt.show()                                                                                 #显示绘制结果 
 
if __name__ == '__main__':
    dataSet, labels = createDataSet()
    #print(dataSet)
    featLabels = []
    myTree = createTree(dataSet, labels, featLabels)
    print(featLabels)
    print(myTree)
    createPlot(myTree)

可视化结果：

使用决策树执行分类

"""
函数说明:使用决策树分类
 
Parameters:
    inputTree - 已经生成的决策树
    featLabels - 存储选择的最优特征标签
    testVec - 测试数据列表，顺序对应最优特征标签
Returns:
    classLabel - 分类结果
"""
def classify(inputTree, featLabels, testVec):
    firstStr = next(iter(inputTree))                                                        #获取决策树结点
    secondDict = inputTree[firstStr]                                                        #下一个字典
    featIndex = featLabels.index(firstStr)                                               
    for key in secondDict.keys():
        if testVec[featIndex] == key:
            if type(secondDict[key]).__name__ == 'dict':
                classLabel = classify(secondDict[key], featLabels, testVec)
            else: classLabel = secondDict[key]
    return classLabel
if __name__ == '__main__':
    dataSet, labels = createDataSet()
    #print(dataSet)
    featLabels = []
    myTree = createTree(dataSet, labels, featLabels)
    print(featLabels)
    print(myTree)
    createPlot(myTree)
    testVec = [0,0,1,1]                                        #测试数据
    result = classify(myTree, featLabels, testVec)
    if result == 0:
        print('hard')
    if result == 1:
        print('no lenese')
    if result == 2:
        print('soft')
        
# 输出“no lenese”

四、朴素贝叶斯算法

4.1 算法简述

朴素贝叶斯算法是有监督的学习算法，解决的是分类问题，如客户是否流失、是否值得投资、信用等级评定等多分类问题。该算法的优点在于简单易懂、学习效率高、在某些领域的分类问题中能够与决策树、神经网络相媲美。但由于该算法以自变量之间的独立（条件特征独立）性和连续变量的正态性假设为前提，就会导致算法精度在某种程度上受影响。

4.2 算法的数学模型

条件概率公式

全概率公式

贝叶斯推断

我们把P(A)称为"先验概率"（Prior probability），即在B事件发生之前，我们对A事件概率的一个判断。

P(A|B)称为"后验概率"（Posterior probability），即在B事件发生之后，我们对A事件概率的重新评估。

P(B|A)/P(B)称为"可能性函数"（Likelyhood），这是一个调整因子，使得预估概率更接近真实概率。

朴素贝叶斯推断

贝叶斯和朴素贝叶斯的概念是不同的，区别就在于“朴素”二字，朴素贝叶斯对条件个概率分布做了条件独立性的假设。 比如下面的公式，假设有n个特征：

由于每个特征都是独立的，我们可以进一步拆分公式 ：

4.4 算法的运用实例

使用朴素贝叶斯解决一些现实生活中的问题时，需要先从文本内容得到字符串列表，然后生成词向量。朴素贝叶斯的一个最著名的应用：电子邮件垃圾过滤。

4.4.1 算法实现的一般流程

• 收集数据：提供文本文件。
• 准备数据：将文本文件解析成词条向量。
• 分析数据：检查词条确保解析的正确性。
• 训练算法：使用我们之前建立的trainNB0()函数。
• 测试算法：使用classifyNB()，并构建一个新的测试函数来计算文档集的错误率。
• 使用算法：构建一个完整的程序对一组文档进行分类，将错分的文档输出到屏幕上。

算法的实现

• 收集数据

  有两个文件夹ham和spam，spam文件下的txt文件为垃圾邮件。

• 数据基本处理

# -*- coding: UTF-8 -*-
import numpy as np
import random
import re

"""
函数说明:将切分的实验样本词条整理成不重复的词条列表，也就是词汇表

Parameters:
    dataSet - 整理的样本数据集
Returns:
    vocabSet - 返回不重复的词条列表，也就是词汇表
"""
def createVocabList(dataSet):
    vocabSet = set([])                      #创建一个空的不重复列表
    for document in dataSet:               
        vocabSet = vocabSet | set(document) #取并集
    return list(vocabSet)

"""
函数说明:根据vocabList词汇表，将inputSet向量化，向量的每个元素为1或0

Parameters:
    vocabList - createVocabList返回的列表
    inputSet - 切分的词条列表
Returns:
    returnVec - 文档向量,词集模型
"""
def setOfWords2Vec(vocabList, inputSet):
    returnVec = [0] * len(vocabList)                                    #创建一个其中所含元素都为0的向量
    for word in inputSet:                                                #遍历每个词条
        if word in vocabList:                                            #如果词条存在于词汇表中，则置1
            returnVec[vocabList.index(word)] = 1
        else: print("the word: %s is not in my Vocabulary!" % word)
    return returnVec                                                    #返回文档向量

• 构建词袋模型

"""
函数说明:根据vocabList词汇表，构建词袋模型

Parameters:
    vocabList - createVocabList返回的列表
    inputSet - 切分的词条列表
Returns:
    returnVec - 文档向量,词袋模型
"""
def bagOfWords2VecMN(vocabList, inputSet):
    returnVec = [0]*len(vocabList)                                        #创建一个其中所含元素都为0的向量
    for word in inputSet:                                                #遍历每个词条
        if word in vocabList:                                            #如果词条存在于词汇表中，则计数加一
            returnVec[vocabList.index(word)] += 1
    return returnVec                                                    #返回词袋模型

• 训练朴素贝叶斯分类器

"""
函数说明:朴素贝叶斯分类器训练函数

Parameters:
    trainMatrix - 训练文档矩阵，即setOfWords2Vec返回的returnVec构成的矩阵
    trainCategory - 训练类别标签向量，即loadDataSet返回的classVec
Returns:
    p0Vect - 非侮辱类的条件概率数组
    p1Vect - 侮辱类的条件概率数组
    pAbusive - 文档属于侮辱类的概率
"""
def trainNB0(trainMatrix,trainCategory):
    numTrainDocs = len(trainMatrix)                            #计算训练的文档数目
    numWords = len(trainMatrix[0])                            #计算每篇文档的词条数
    pAbusive = sum(trainCategory)/float(numTrainDocs)        #文档属于侮辱类的概率
    p0Num = np.ones(numWords); p1Num = np.ones(numWords)    #创建numpy.ones数组,词条出现数初始化为1，拉普拉斯平滑
    p0Denom = 2.0; p1Denom = 2.0                            #分母初始化为2,拉普拉斯平滑
    for i in range(numTrainDocs):
        if trainCategory[i] == 1:                            #统计属于侮辱类的条件概率所需的数据，即P(w0|1),P(w1|1),P(w2|1)···
            p1Num += trainMatrix[i]
            p1Denom += sum(trainMatrix[i])
        else:                                                #统计属于非侮辱类的条件概率所需的数据，即P(w0|0),P(w1|0),P(w2|0)···
            p0Num += trainMatrix[i]
            p0Denom += sum(trainMatrix[i])
    p1Vect = np.log(p1Num/p1Denom)                            #取对数，防止下溢出         
    p0Vect = np.log(p0Num/p0Denom)         
    return p0Vect,p1Vect,pAbusive                            #返回属于侮辱类的条件概率数组，属于非侮辱类的条件概率数组，文档属于侮辱类的概率

• 执行分类

"""
函数说明:朴素贝叶斯分类器分类函数

Parameters:
    vec2Classify - 待分类的词条数组
    p0Vec - 非侮辱类的条件概率数组
    p1Vec -侮辱类的条件概率数组
    pClass1 - 文档属于侮辱类的概率
Returns:
    0 - 属于非侮辱类
    1 - 属于侮辱类
"""
def classifyNB(vec2Classify, p0Vec, p1Vec, pClass1):
    p1 = sum(vec2Classify * p1Vec) + np.log(pClass1)        #对应元素相乘。logA * B = logA + logB，所以这里加上log(pClass1)
    p0 = sum(vec2Classify * p0Vec) + np.log(1.0 - pClass1)
    if p1 > p0:
        return 1
    else:
        return 0

"""
函数说明:接收一个大字符串并将其解析为字符串列表

Parameters:
    无
Returns:
    无
"""
def textParse(bigString):                                                   #将字符串转换为字符列表
    listOfTokens = re.split(r'\W*', bigString)                              #将特殊符号作为切分标志进行字符串切分，即非字母、非数字
    return [tok.lower() for tok in listOfTokens if len(tok) > 2]            #除了单个字母，例如大写的I，其它单词变成小写

• 测试分类情况

"""
函数说明:测试朴素贝叶斯分类器

Parameters:
    无
Returns:
    无
"""
def spamTest():
    docList = []; classList = []; fullText = []
    for i in range(1, 26):                                                  #遍历25个txt文件
        wordList = textParse(open('data/Naive_Bayes/spam/%d.txt' % i, 'r').read())     #读取每个垃圾邮件，并字符串转换成字符串列表
        docList.append(wordList)
        fullText.append(wordList)
        classList.append(1)                                                 #标记垃圾邮件，1表示垃圾文件
        wordList = textParse(open('data/Naive_Bayes/ham/%d.txt' % i, 'r').read())      #读取每个非垃圾邮件，并字符串转换成字符串列表
        docList.append(wordList)
        fullText.append(wordList)
        classList.append(0)                                                 #标记非垃圾邮件，1表示垃圾文件   
    vocabList = createVocabList(docList)                                    #创建词汇表，不重复
    trainingSet = list(range(50)); testSet = []                             #创建存储训练集的索引值的列表和测试集的索引值的列表                       
    for i in range(10):                                                     #从50个邮件中，随机挑选出40个作为训练集,10个做测试集
        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(np.array(trainMat), np.array(trainClasses))  #训练朴素贝叶斯模型
    errorCount = 0                                                          #错误分类计数
    for docIndex in testSet:                                                #遍历测试集
        wordVector = setOfWords2Vec(vocabList, docList[docIndex])           #测试集的词集模型
        if classifyNB(np.array(wordVector), p0V, p1V, pSpam) != classList[docIndex]:    #如果分类错误
            errorCount += 1                                                 #错误计数加1
            print("分类错误的测试集：",docList[docIndex])
    print('错误率：%.2f%%' % (float(errorCount) / len(testSet) * 100))


if __name__ == '__main__':
    spamTest()

运行结果：

五、Logistic回归

5.1 算法简述

Logistic回归是分类方法，它利用的是Sigmoid函数阈值在[0,1]这个特性。Logistic回归进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，以此进行分类。其实，Logistic本质上是一个基于条件概率的判别模型(Discriminative Model)。

就比如下面这个散点图，我们要将其分为两类，很显然我们可以通过画一条直线来实现。

所以我们的目的是去计算出这条直线的方程，下面我们从数学的角度来分析这个问题。

5.2 算法的数学模型

要求拟合直线，假设是二维的情况下，我们可以假设直线的方程：W0+W1·X1+W2·X2=0，这里的x2相当于y值。

逻辑回归法便是要用已知的x值取训练处这里的待定系数w,那么问题来了，我们怎么训练呢？

于是我们引入了Sigmoid函数，我们也可以称它为Logistic函数。它的公式如下：

整合公式得出了如下公式：

下面这张图，为我们展示了Sigmoid函数的样子：

z是一个矩阵，θ是参数列向量(要求解的)，x是样本列向量(给定的数据集)。θ^T表示θ的转置。g(z)函数实现了任意实数到[0,1]的映射，这样我们的数据集([x0,x1,…,xn])，不管是大于1或者小于0，都可以映射到[0,1]区间进行分类。

如果我们有合适的参数列向量θ([θ0,θ1,…θn]^T)，以及样本列向量x([x0,x1,…,xn])，那么我们对样本x分类就可以通过上述公式计算出一个概率，如果这个概率大于0.5，我们就可以说样本是正样本，否则样本是负样本。

那么问题来了！如何得到合适的参数向量θ?

根据sigmoid函数的特性，我们可以做出如下的假设：

我们可以把上述两个概率公式合二为一：

合并出来的Cost，我们称之为代价函数(CostFunction)。当y等于1时，(1-y)项(第二项)为0；当y等于0时，y项(第一项)为0。为了简化问题，我们对整个表达式求对数，(将指数问题对数化是处理数学问题常见的方法)：

这个代价函数，是对于一个样本而言的。给定一个样本，我们就可以通过这个代价函数求出，样本所属类别的概率，而这个概率越大越好，所以也就是求解这个代价函数的最大值。既然概率出来了，那么最大似然估计也该出场了。假定样本与样本之间相互独立，那么整个样本集生成的概率即为所有样本生成概率的乘积，再将公式对数化，便可得到如下公式：

其中，m为样本的总数，y(i)表示第i个样本的类别，x(i)表示第i个样本，需要注意的是θ是多维向量，x(i)也是多维向量。

综上所述，满足J(θ)的最大的θ值即是我们需要求解的模型。

怎么求解使J(θ)最大的θ值呢？因为是求最大值，所以我们需要使用梯度上升算法。如果面对的问题是求解使J(θ)最小的θ值，那么我们就需要使用梯度下降算法。

梯度上升算法

就像爬坡一样，一点一点逼近极值。这种寻找最佳拟合参数的方法，就是最优化算法。爬坡这个动作用数学公式表达即为：

其中，α为步长，也就是学习速率，控制更新的幅度。

根据推到，得出该函数的梯度上升迭代公式为：

将上述公式矢量化：

算法的实现

实现sigmoid函数

"""
函数说明:sigmoid函数
 
Parameters:
    inX - 数据
Returns:
    sigmoid函数
"""
def sigmoid(inX):
    return 1.0 / (1 + np.exp(-inX))

实现梯度上升法

"""
函数说明:梯度上升算法
 
Parameters:
    dataMatIn - 数据集
    classLabels - 数据标签
Returns:
    weights.getA() - 求得的权重数组(最优参数)
"""
def gradAscent(dataMatIn, classLabels):
    dataMatrix = np.mat(dataMatIn)                                        #转换成numpy的mat
    labelMat = np.mat(classLabels).transpose()                            #转换成numpy的mat,并进行转置
    m, n = np.shape(dataMatrix)                                            #返回dataMatrix的大小。m为行数,n为列数。
    alpha = 0.001                                                        #移动步长,也就是学习速率,控制更新的幅度。
    maxCycles = 500                                                        #最大迭代次数
    weights = np.ones((n,1))
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights)                                #梯度上升矢量化公式
        error = labelMat - h
        weights = weights + alpha * dataMatrix.transpose() * error
    return weights.getA()                                                #将矩阵转换为数组，返回权重数组

绘制决策边界

我们已经解出了一组回归系数，它确定了不同类别数据之间的分隔线。现在开始绘制这个分隔线，编写代码如下：

"""
函数说明:绘制数据集
 
Parameters:
    weights - 权重参数数组
Returns:
    无
"""
def plotBestFit(weights):
    dataMat, labelMat = loadDataSet()                                    #加载数据集
    dataArr = np.array(dataMat)                                            #转换成numpy的array数组
    n = np.shape(dataMat)[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])    #1为正样本
        else:
            xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2])    #0为负样本
    fig = plt.figure()
    ax = fig.add_subplot(111)                                            #添加subplot
    ax.scatter(xcord1, ycord1, s = 20, c = 'red', marker = 's',alpha=.5)#绘制正样本
    ax.scatter(xcord2, ycord2, s = 20, c = 'green',alpha=.5)            #绘制负样本
    x = np.arange(-3.0, 3.0, 0.1)
    y = (-weights[0] - weights[1] * x) / weights[2]
    ax.plot(x, y)
    plt.title('BestFit')                                                #绘制title
    plt.xlabel('X1'); plt.ylabel('X2')                                    #绘制label
    plt.show()       

算法的改进

该方法处理100个左右的数据集时尚可，但如果有数十亿样本和成千上万的特征，那么该方法的计算复杂度就太高了。因此，需要对算法进行改进，我们每次更新回归系数(最优参数)的时候，能不能不用所有样本呢？一次只用一个样本点去更新回归系数(最优参数)？这样就可以有效减少计算量了，这种方法就叫做随机梯度上升算法。

随机梯度上升算法

"""
函数说明:改进的随机梯度上升算法
 
Parameters:
    dataMatrix - 数据数组
    classLabels - 数据标签
    numIter - 迭代次数
Returns:
    weights - 求得的回归系数数组(最优参数)

"""
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    m,n = np.shape(dataMatrix)                                                #返回dataMatrix的大小。m为行数,n为列数。
    weights = np.ones(n)                                                       #参数初始化
    for j in range(numIter):                                           
        dataIndex = list(range(m))
        for i in range(m):           
            alpha = 4/(1.0+j+i)+0.01                                            #降低alpha的大小，每次减小1/(j+i)。
            randIndex = int(random.uniform(0,len(dataIndex)))                #随机选取样本
            h = sigmoid(sum(dataMatrix[randIndex]*weights))                    #选择随机选取的一个样本，计算h
            error = classLabels[randIndex] - h                                 #计算误差
            weights = weights + alpha * error * dataMatrix[randIndex]       #更新回归系数
            del(dataIndex[randIndex])                                         #删除已经使用的样本
    return weights                                                      #返回

该算法第一个改进之处在于，alpha在每次迭代的时候都会调整，并且，虽然alpha会随着迭代次数不断减小，但永远不会减小到0，因为这里还存在一个常数项。必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。如果需要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获得更大的回归系数。另一点值得注意的是，在降低alpha的函数中，alpha每次减少1/(j+i)，其中j是迭代次数，i是样本点的下标。第二个改进的地方在于更新回归系数(最优参数)时，只使用一个样本点，并且选择的样本点是随机的，每次迭代不使用已经用过的样本点。这样的方法，就有效地减少了计算量，并保证了回归效果。

六、支持向量机算法

6.1 算法简述

支持向量机（support vector machines）是一种二分类模型，它的目的是寻找一个超平面来对样本进行分割，分割的原则是间隔最大化，最终转化为一个凸二次规划问题来求解。由简至繁的模型包括：

当训练样本线性可分时，通过硬间隔最大化，学习一个线性可分支持向量机；

当训练样本近似线性可分时，通过软间隔最大化，学习一个线性支持向量机；

当训练样本线性不可分时，通过核技巧和软间隔最大化，学习一个非线性支持向量机；

6.2 算法的数学模型

不管是训练样本是线性可分还是线性不可分的，SVM算法的目的就是寻找一个超平面来将样本分开，所以求解这个“决策面”的过程就是寻找最优化解的过程，在数学上我们求解最优解的思想是：首先，找出目标函数；其次，明确优化对象，也就是你希望通过哪些因数来使你的目标函数达到最优。在线性SVM中，目标函数显然是“分类间隔”，而优化对象则是“决策面”。所以，我们先对上述的两个要点进行数学建模。

按照我们的常规思维，我们先在二维空间上建模，然后再推广到多维。这里需要注意的是我们的思维不要总想着在二维坐标上横轴是X，纵轴是y,其实它们只是代号而已，我们正确的理解是每个变量它只是代表一个维度的变化。而由这些维度的变化我们希望得到的结果就可以用y表示。

详细推导过程,参照文档

线性SVM

决策面的方程

二维情况下的决策面方程

将二维情况推广到多维后方程不变，只是系数w和x变了：

我们已经顺利推导出了"决策面"方程，它就是我们的超平面方程，之后，我们统称其为超平面方程。

分类间隔方程

我们依然对于一个二维平面的简单例子进行推导。

支持向量到决策面的距离：

其中：

于是我们可以将其推广到多维情况。

我们目的是为了找出一个分类效果好的超平面作为分类器。分类器的好坏的评定依据是分类间隔W=2d的大小，即分类间隔w越大，我们认为这个超平面的分类效果越好。此时，求解超平面的问题就变成了求解分类间隔W最大化的为题。W的最大化也就是d最大化的。

约束条件

我们已经顺利获得了目标函数的数学形式。但是为了求解w的最大值。我们不得不面对如下问题：

我们如何判断超平面是否将样本点正确分类？

我们知道要求距离d的最大值，我们首先需要找到支持向量上的点，怎么在众多的点中选出支持向量上的点呢？

用数学公式描述：

对每个样本点xi加上一个类别标签yi：

假设决策面正好处于间隔区域的中轴线上，并且相应的支持向量对应的样本点到决策面的距离为d，那么公式进一步写成：

对于所有分类标签为1和-1样本点，它们到直线的距离都大于等于d(支持向量上的样本点到超平面的距离)。公式两边都除以d.
两个公式，都是描述一条直线，数学模型代表的意义是一样的。现在，让我们对wd和γd重新起个名字，就叫它们w和γ。因此，我们就可以说：“对于存在分类间隔的两类样本点，我们一定可以找到一些超平面，使其对于所有的样本点均满足下面的条件：”

进一步优化得出目标函数和约束条件

现在我们就可以将我们的目标函数进一步化简：

因为，我们只关心支持向量上的点。随后我们求解d的最大化问题变成了||w||的最小化问题。进而||w||的最小化问题等效于

为什么要做这样的等效呢？这是为了在进行最优化的过程中对目标函数求导时比较方便，但这绝对不影响最优化问题最后的求解。我们将最终的目标函数和约束条件放在一起进行描述：

这里n是样本点的总个数，缩写s.t.表示"Subject to"，是"服从某某条件"的意思。上述公式描述的是一个典型的不等式约束条件下的二次型函数优化问题，同时也是支持向量机的基本数学模型。

求解目标函数

要想求解上面所得的数学模型我们需要引入拉格朗日函数，为什么要这样做呢？

我们知道我们要求解的是最小化问题，所以一个直观的想法是如果我能够构造一个函数，使得该函数在可行解区域内与原目标函数完全一致，而在可行解区域外的数值非常大，甚至是无穷大，那么这个没有约束条件的新目标函数的优化问题就与原来有约束条件的原始目标函数的优化问题是等价的问题。这就是使用拉格朗日方程的目的，它将约束条件放到目标函数中，从而将有约束优化问题转换为无约束优化问题。

• 将有约束的原始目标函数转换为无约束的新构造的拉格朗日目标函数

其中αi是拉格朗日乘子，αi大于等于0，于是我们构造新目标函数时引入的系数变量(我们自己设置)。现在我们令：

此时，显然θ(w)为原目标函数本身。我们将上述两种情况结合一下，就得到了新的目标函数：

此时，再看我们的初衷，就是为了建立一个在可行解区域内与原目标函数相同，在可行解区域外函数值趋近于无穷大的新函数，现在我们做到了。

现在，我们的问题变成了求新目标函数的最小值，即：

这里用p*表示这个问题的最优值，且和最初的问题是等价的。

• 将不易求解的优化问题转化为易求解的优化

我们看一下我们的新目标函数，先求最大值，再求最小值。这样的话，我们首先就要面对带有需要求解的参数w和b的方程，而αi又是不等式约束，这个求解过程不好做。所以，我们需要使用拉格朗日函数对偶性，将最小和最大的位置交换一下，这样就变成了：

交换以后的新问题是原始问题的对偶问题，这个新问题的最优值用d来表示。而且d<=p*。我们关心的是d=p的时候，这才是我们要的解。需要什么条件才能让d=p呢？

首先必须满足这个优化问题是凸优化问题。
其次，需要满足KKT条件。
凸优化问题的定义是：求取最小值的目标函数为凸函数的一类优化问题。目标函数是凸函数我们已经知道，这个优化问题又是求最小值。所以我们的最优化问题就是凸优化问题。

接下里，就是探讨是否满足KKT条件了。

• KKT条件

我们已经使用拉格朗日函数对我们的目标函数进行了处理，生成了一个新的目标函数。通过一些条件，可以求出最优值的必要条件，这个条件就是接下来要说的KKT条件。一个最优化模型能够表示成下列标准形式：

根据推导目标函数满足KKT条件，所以问题转换成了对偶问题。而求解这个对偶学习问题，可以分为三个步骤：首先要让L(w,b,α)关于w和b最小化，然后求对α的极大，最后利用SMO算法求解对偶问题中的拉格朗日乘子。现在，我们继续推导。

• 对偶问题求解
  根据上述推导得出：

所以我们要先求出内侧的min,根据拉格朗日函数求解最优解的方法：

首先固定α，要让L(w,b,α)关于w和b最小化，我们分别对w和b偏导数，令其等于0，即：

将上述结果带回L(w,b,α)得到：

从上面的最后一个式子，我们可以看出，此时的L(w,b,α)函数只含有一个变量，即αi。

然后求解外侧的max,将上面的min带入得:

现在我们的优化问题变成了如上的形式。对于这个问题，我们有更高效的优化算法，即序列最小优化（SMO）算法。我们通过这个优化算法能得到α，再根据α，我们就可以求解出w和b，进而求得我们最初的目的：找到超平面，即"决策平面"。我们上述所有推导的目的就是将目标函数转化为只含有一个未知数的目标函数，因为这样我们才可以用SMO算法求解。

• SMO算法

SMO算法的目标是求出一系列alpha和b，一旦求出了这些alpha，就很容易计算出权重向量w并得到分隔超平面。

先来定义特征到结果的输出函数为：

将上面推出的w带入输出函数：

将目标函数变形，在前面增加一个符号，将最大值问题转换成最小值问题：

实际上，对于上述目标函数，是存在一个假设的，即数据100%线性可分。但是，目前为止，我们知道几乎所有数据都不那么"干净"。这时我们就可以通过引入所谓的松弛变量(slack variable)，来允许有些数据点可以处于超平面的错误的一侧。这样我们的优化目标就能保持仍然不变，但是此时我们的约束条件有所改变：

• 经过推导我们可以归纳出SMO算法的步骤：

  • 步骤1：计算误差：

    

  • 步骤2：计算上下界L和H：

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Wzc5SLS5-1577762322622)(http://www.cuijiahua.com/wp-content/uploads/2017/11/ml_8_107.png)]

  • 步骤3：计算η：

    

  • 步骤4：更新αj：

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-mjgiM54i-1577762322625)(http://www.cuijiahua.com/wp-content/uploads/2017/11/ml_8_109.png)]

  • 步骤5：根据取值范围修剪αj：

    

  • 步骤6：更新αi：

    

  • 步骤7：更新b1和b2：

    

  • 步骤8：根据b1和b2更新b：

    

SMO算法的优化

在几百个点组成的小规模数据集上，简化版SMO算法的运行是没有什么问题的，但是在更大的数据集上的运行速度就会变慢。简化版SMO算法的第二个α的选择是随机的，针对这一问题，我们可以使用启发式选择第二个α值，来达到优化效果。

下面这两个公式想必已经不再陌生：

在实现SMO算法的时候，先计算η，再更新αj。为了加快第二个αj乘子的迭代速度，需要让直线的斜率增大，对于αj的更新公式，其中η值没有什么文章可做，于是只能令:

因此，我们可以明确自己的优化方法了：

最外层循环，首先在样本中选择违反KKT条件的一个乘子作为最外层循环，然后用"启发式选择"选择另外一个乘子并进行这两个乘子的优化
在非边界乘子中寻找使得 |Ei - Ej| 最大的样本
如果没有找到，则从整个样本中随机选择一个样本
接下来，让我们看看完整版SMO算法如何实现。

算法的实现

需要引用的包

import matplotlib.pyplot as plt
import numpy as np
import random

建立存储数据的的数据结构

class optStruct:
    """
    数据结构，维护所有需要操作的值
    Parameters：
        dataMatIn - 数据矩阵
        classLabels - 数据标签
        C - 松弛变量
        toler - 容错率
    """
    def __init__(self, dataMatIn, classLabels, C, toler):
        self.X = dataMatIn                                #数据矩阵
        self.labelMat = classLabels                        #数据标签
        self.C = C                                         #松弛变量
        self.tol = toler                                 #容错率
        self.m = np.shape(dataMatIn)[0]                 #数据矩阵行数
        self.alphas = np.mat(np.zeros((self.m,1)))         #根据矩阵行数初始化alpha参数为0   
        self.b = 0                                         #初始化b参数为0
        self.eCache = np.mat(np.zeros((self.m,2)))         #根据矩阵行数初始化虎误差缓存，第一列为是否有效的标志位，第二列为实际的误差E的值。

加载数据集到内存

def loadDataSet(fileName):
    """
    读取数据
    Parameters:
        fileName - 文件名
    Returns:
        dataMat - 数据矩阵
        labelMat - 数据标签
    """
    dataMat = []; labelMat = []
    fr = open(fileName)
    for line in fr.readlines():                                     #逐行读取，滤除空格等
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])      #添加数据
        labelMat.append(float(lineArr[2]))                          #添加标签
    return dataMat,labelMat

计算误差

def calcEk(oS, k):
    """
    计算误差
    Parameters：
        oS - 数据结构
        k - 标号为k的数据
    Returns:
        Ek - 标号为k的数据误差
    """
    fXk = float(np.multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T) + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek

随机选择alpha_j的索引值

def selectJrand(i, m):
    """
    函数说明:随机选择alpha_j的索引值
 
    Parameters:
        i - alpha_i的索引值
        m - alpha参数个数
    Returns:
        j - alpha_j的索引值
    """
    j = i                                 #选择一个不等于i的j
    while (j == i):
        j = int(random.uniform(0, m))
    return j

以内循环启发方式选择alpha_j

def selectJ(i, oS, Ei):
    """
    内循环启发方式2
    Parameters：
        i - 标号为i的数据的索引值
        oS - 数据结构
        Ei - 标号为i的数据误差
    Returns:
        j, maxK - 标号为j或maxK的数据的索引值
        Ej - 标号为j的数据误差
    """
    maxK = -1; maxDeltaE = 0; Ej = 0                         #初始化
    oS.eCache[i] = [1,Ei]                                      #根据Ei更新误差缓存
    validEcacheList = np.nonzero(oS.eCache[:,0].A)[0]        #返回误差不为0的数据的索引值
    if (len(validEcacheList)) > 1:                            #有不为0的误差
        for k in validEcacheList:                           #遍历,找到最大的Ek
            if k == i: continue                             #不计算i,浪费时间
            Ek = calcEk(oS, k)                                #计算Ek
            deltaE = abs(Ei - Ek)                            #计算|Ei-Ek|
            if (deltaE > maxDeltaE):                        #找到maxDeltaE
                maxK = k; maxDeltaE = deltaE; Ej = Ek
        return maxK, Ej                                        #返回maxK,Ej
    else:                                                   #没有不为0的误差
        j = selectJrand(i, oS.m)                            #随机选择alpha_j的索引值
        Ej = calcEk(oS, j)                                    #计算Ej
    return j, Ej                                             #j,Ej

计算误差Ek，并更新误差缓存

def updateEk(oS, k):
    """
    计算Ek,并更新误差缓存
    Parameters：
        oS - 数据结构
        k - 标号为k的数据的索引值
    Returns:
        无
    """
    Ek = calcEk(oS, k)                                        #计算Ek
    oS.eCache[k] = [1,Ek]                                    #更新误差缓存

修剪alpha_j

ef clipAlpha(aj,H,L):
    """
    修剪alpha_j
    Parameters:
        aj - alpha_j的值
        H - alpha上限
        L - alpha下限
    Returns:
        aj - 修剪后的alpah_j的值
    """
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

优化的SMO算法

def innerL(i, oS):
    """
    优化的SMO算法
    Parameters：
        i - 标号为i的数据的索引值
        oS - 数据结构
    Returns:
        1 - 有任意一对alpha值发生变化
        0 - 没有任意一对alpha值发生变化或变化太小
    """
    #步骤1：计算误差Ei
    Ei = calcEk(oS, i)
    #优化alpha,设定一定的容错率。
    if ((oS.labelMat[i] * Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i] * Ei > oS.tol) and (oS.alphas[i] > 0)):
        #使用内循环启发方式2选择alpha_j,并计算Ej
        j,Ej = selectJ(i, oS, Ei)
        #保存更新前的aplpha值，使用深拷贝
        alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy();
        #步骤2：计算上下界L和H
        if (oS.labelMat[i] != oS.labelMat[j]):
            L = max(0, oS.alphas[j] - oS.alphas[i])
            H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
        else:
            L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
            H = min(oS.C, oS.alphas[j] + oS.alphas[i])
        if L == H:
            print("L==H")
            return 0
        #步骤3：计算eta
        eta = 2.0 * oS.X[i,:] * oS.X[j,:].T - oS.X[i,:] * oS.X[i,:].T - oS.X[j,:] * oS.X[j,:].T
        if eta >= 0:
            print("eta>=0")
            return 0
        #步骤4：更新alpha_j
        oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej)/eta
        #步骤5：修剪alpha_j
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        #更新Ej至误差缓存
        updateEk(oS, j)
        if (abs(oS.alphas[j] - alphaJold) < 0.00001):
            print("alpha_j变化太小")
            return 0
        #步骤6：更新alpha_i
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])
        #更新Ei至误差缓存
        updateEk(oS, i)
        #步骤7：更新b_1和b_2
        b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[i,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[i,:]*oS.X[j,:].T
        b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[j,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[j,:]*oS.X[j,:].T
        #步骤8：根据b_1和b_2更新b
        if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1
        elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2
        else: oS.b = (b1 + b2)/2.0
        return 1
    else:
        return 0

完整的SMO算法

def smoP(dataMatIn, classLabels, C, toler, maxIter):
    """
    完整的线性SMO算法
    Parameters：
        dataMatIn - 数据矩阵
        classLabels - 数据标签
        C - 松弛变量
        toler - 容错率
        maxIter - 最大迭代次数
    Returns:
        oS.b - SMO算法计算的b
        oS.alphas - SMO算法计算的alphas
    """
    oS = optStruct(np.mat(dataMatIn), np.mat(classLabels).transpose(), C, toler)                    #初始化数据结构
    iter = 0                                                                                         #初始化当前迭代次数
    entireSet = True; alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):                            #遍历整个数据集都alpha也没有更新或者超过最大迭代次数,则退出循环
        alphaPairsChanged = 0
        if entireSet:                                                                                #遍历整个数据集                           
            for i in range(oS.m):       
                alphaPairsChanged += innerL(i,oS)                                                    #使用优化的SMO算法
                print("全样本遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter,i,alphaPairsChanged))
            iter += 1
        else:                                                                                         #遍历非边界值
            nonBoundIs = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]                        #遍历不在边界0和C的alpha
            for i in nonBoundIs:
                alphaPairsChanged += innerL(i,oS)
                print("非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter,i,alphaPairsChanged))
            iter += 1
        if entireSet:                                                                                #遍历一次后改为非边界遍历
            entireSet = False
        elif (alphaPairsChanged == 0):                                                                #如果alpha没有更新,计算全样本遍历
            entireSet = True 
        print("迭代次数: %d" % iter)
    return oS.b,oS.alphas                                                                             #返回SMO算法计算的b和alphas

分类结果可视化

def showClassifer(dataMat, classLabels, w, b):
    """
    分类结果可视化
    Parameters:
        dataMat - 数据矩阵
        w - 直线法向量
        b - 直线解决
    Returns:
        无
    """
    #绘制样本点
    data_plus = []                                  #正样本
    data_minus = []                                 #负样本
    for i in range(len(dataMat)):
        if classLabels[i] > 0:
            data_plus.append(dataMat[i])
        else:
            data_minus.append(dataMat[i])
    data_plus_np = np.array(data_plus)              #转换为numpy矩阵
    data_minus_np = np.array(data_minus)            #转换为numpy矩阵
    plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1], s=30, alpha=0.7)   #正样本散点图
    plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1], s=30, alpha=0.7) #负样本散点图
    #绘制直线
    x1 = max(dataMat)[0]
    x2 = min(dataMat)[0]
    a1, a2 = w
    b = float(b)
    a1 = float(a1[0])
    a2 = float(a2[0])
    y1, y2 = (-b- a1*x1)/a2, (-b - a1*x2)/a2
    plt.plot([x1, x2], [y1, y2])
    #找出支持向量点
    for i, alpha in enumerate(alphas):
        if abs(alpha) > 0:
            x, y = dataMat[i]
            plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
    plt.show()

计算W值

def calcWs(alphas,dataArr,classLabels):
    """
    计算w
    Parameters:
        dataArr - 数据矩阵
        classLabels - 数据标签
        alphas - alphas值
    Returns:
        w - 计算得到的w
    """
    X = np.mat(dataArr); labelMat = np.mat(classLabels).transpose()
    m,n = np.shape(X)
    w = np.zeros((n,1))
    for i in range(m):
        w += np.multiply(alphas[i]*labelMat[i],X[i,:].T)
    return w

main函数

if __name__ == '__main__':
    dataArr, classLabels = loadDataSet('testSet.txt')
    b, alphas = smoP(dataArr, classLabels, 0.6, 0.001, 40)
    w = calcWs(alphas,dataArr, classLabels)
    showClassifer(dataArr, classLabels, w, b)

运行结果

非线性SVM

我们已经了解到，SVM如何处理线性可分的情况，而对于非线性的情况，SVM的处理方式就是选择一个核函数。简而言之：在线性不可分的情况下，SVM通过某种事先选择的非线性映射（核函数）将输入变量映到一个高维特征空间，将其变成在高维空间线性可分，在这个高维空间中构造最优分类超平面。

数学模型

核技巧

根据线性SVM的推导，我们可知线性可分的情况下，可知最终的超平面方程为：

为了表述方便我们将上述公式用内积来表示：

对于线性不可分，我们使用一个非线性映射，将数据映射到特征空间，在特征空间中使用线性学习器，分类函数变形如下：

其中ϕ从输入空间(X)到某个特征空间(F)的映射，这意味着建立非线性学习器分为两步：

首先使用一个非线性映射将数据变换到一个特征空间F；
然后在特征空间使用线性学习器分类。
如果有一种方法可以在特征空间中直接计算内积 <ϕ(xi),ϕ(x)>，就像在原始输入点的函数中一样，就有可能将两个步骤融合到一起建立一个分线性的学习器，这样直接计算的方法称为核函数方法。
这种将内积替换成核函数的方式被称为核技巧(kernel trick)。

核技巧的实现

通过核技巧的转变，我们的分类函数变为：

所以，我们的对偶问题变成了：

这样，我们就避开了高纬度空间中的计算。当然，我们刚刚的例子是非常简单的，我们可以手动构造出来对应映射的核函数出来，如果对于任意一个映射，要构造出对应的核函数就很困难了。因此，通常，人们会从一些常用的核函数中进行选择，根据问题和数据的不同，选择不同的参数，得到不同的核函数。接下来，要介绍的就是一个非常流行的核函数，那就是径向基核函数。

径向基核函数是SVM中常用的一个核函数。径向基核函数采用向量作为自变量的函数，能够基于向量举例运算输出一个标量。径向基核函数的高斯版本的公式如下：

其中，σ是用户自定义的用于确定到达率(reach)或者说函数值跌落到0的速度参数。上述高斯核函数将数据从原始空间映射到无穷维空间。关于无穷维空间，我们不必太担心。高斯核函数只是一个常用的核函数，使用者并不需要确切地理解数据到底是如何表现的，而且使用高斯核函数还会得到一个理想的结果。如果σ选得很大的话，高次特征上的权重实际上衰减得非常快，所以实际上（数值上近似一下）相当于一个低维的子空间；反过来，如果σ选得很小，则可以将任意的数据映射为线性可分——当然，这并不一定是好事，因为随之而来的可能是非常严重的过拟合问题。不过，总的来说，通过调控参数σ，高斯核实际上具有相当高的灵活性，也是使用最广泛的核函数之一。

算法实现

相比于线性SVM只有小部分的改变

建立数据结构，维护所需要操作的值

class optStruct:
    """
    数据结构，维护所有需要操作的值
    Parameters：
        dataMatIn - 数据矩阵
        classLabels - 数据标签
        C - 松弛变量
        toler - 容错率
        kTup - 包含核函数信息的元组,第一个参数存放核函数类别，第二个参数存放必要的核函数需要用到的参数
    """
    def __init__(self, dataMatIn, classLabels, C, toler, kTup):
        self.X = dataMatIn                                #数据矩阵
        self.labelMat = classLabels                        #数据标签
        self.C = C                                         #松弛变量
        self.tol = toler                                 #容错率
        self.m = np.shape(dataMatIn)[0]                 #数据矩阵行数
        self.alphas = np.mat(np.zeros((self.m,1)))         #根据矩阵行数初始化alpha参数为0   
        self.b = 0                                         #初始化b参数为0
        self.eCache = np.mat(np.zeros((self.m,2)))         #根据矩阵行数初始化虎误差缓存，第一列为是否有效的标志位，第二列为实际的误差E的值。
        self.K = np.mat(np.zeros((self.m,self.m)))        #初始化核K
        for i in range(self.m):                            #计算所有数据的核K
            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)

核函数的实现

def kernelTrans(X, A, kTup):
    """
    通过核函数将数据转换更高维的空间
    Parameters：
        X - 数据矩阵
        A - 单个数据的向量
        kTup - 包含核函数信息的元组
    Returns:
        K - 计算的核K
    """
    m,n = np.shape(X)
    K = np.mat(np.zeros((m,1)))
    if kTup[0] == 'lin': K = X * A.T                       #线性核函数,只进行内积。
    elif kTup[0] == 'rbf':                                 #高斯核函数,根据高斯核函数公式进行计算
        for j in range(m):
            deltaRow = X[j,:] - A
            K[j] = deltaRow*deltaRow.T
        K = np.exp(K/(-1*kTup[1]**2))                     #计算高斯核K
    else: raise NameError('核函数无法识别')
    return K                                             #返回计算的核K

非线性SVM的测试函数

def testRbf(k1 = 1.2):
    """
    测试函数
    Parameters:
        k1 - 使用高斯核函数的时候表示到达率
    Returns:
        无
    """
    dataArr,labelArr = loadDataSet('testSetRBF.txt')                        #加载训练集
    b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 100, ('rbf', k1))        #根据训练集计算b和alphas
    datMat = np.mat(dataArr); labelMat = np.mat(labelArr).transpose()
    svInd = np.nonzero(alphas.A > 0)[0]                                        #获得支持向量
    sVs = datMat[svInd]                                                     
    labelSV = labelMat[svInd];
    print("支持向量个数:%d" % np.shape(sVs)[0])
    m,n = np.shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))                #计算各个点的核
        predict = kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b     #根据支持向量的点，计算超平面，返回预测结果
        if np.sign(predict) != np.sign(labelArr[i]): errorCount += 1        #返回数组中各元素的正负符号，用1和-1表示，并统计错误个数
    print("训练集错误率: %.2f%%" % ((float(errorCount)/m)*100))   #打印错误率
    dataArr,labelArr = loadDataSet('testSetRBF2.txt')            #加载测试集
    errorCount = 0
    datMat = np.mat(dataArr); labelMat = np.mat(labelArr).transpose()         
    m,n = np.shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))   #计算各个点的核           
        predict=kernelEval.T * np.multiply(labelSV,alphas[svInd]) + b         #根据支持向量的点，计算超平面，返回预测结果
        if np.sign(predict) != np.sign(labelArr[i]): errorCount += 1        #返回数组中各元素的正负符号，用1和-1表示，并统计错误个数
    print("测试集错误率: %.2f%%" % ((float(errorCount)/m)*100))             #打印错误率

测试结果

可以看到，训练集错误率为1%，测试集错误率都是8%。可以尝试更换不同的K1参数以观察测试错误率、训练错误率、支持向量个数随k1的变化情况。你会发现K1过大，会出现过拟合的情况，即训练集错误率低，但是测试集错误率高。

七、Adaboost算法

7.1 算法简述

AdaBoost算法是基于Boosting思想的机器学习算法，AdaBoost是adaptive boosting（自适应boosting）的缩写

Boosting的思路是采用重赋权（re-weighting）法迭代地训练基分类器，主要思想：

每一轮的训练数据样本赋予一个权重，并且每一轮样本的权值分布依赖上一轮的分类结果。

基分类器之间采用序列式的线性加权方式进行组合。

7.2 算法的数学模型

AdaBoost可拆分为如下几个过程：

• 计算样本权值

  对已知的训练样本赋予权重，即样本权重，用向量D表示，这些权重都初始化成相等值。假设有n个样本的训练集：

  

  开始迭代计算时，设定每个样本的权重都是相等的，即1/n。

• 计算错误率

  利用第一个弱学习算法h1对其进行学习，学习完成后进行错误率ε的统计：

  

• 计算弱学习算法权重

  每次迭代后会有一个错误率（最低错误率），我们利用这个错误率计算出下一次迭代所需的a

  

• 更新样本权重

  在每一次学习完成后，需要重新调整样本的权重，以使得在前一次分类中被错分的样本的权重，在接下来的学习中可以重点对其进行学习：

  

  其中，ht(xi) = yi表示对第i个样本训练正确，不等于则表示分类错误。Zt是一个归一化因子：

  

  这个公式我们可以继续化简，将两个公式进行合并，化简如下：

  

7.4 算法的运用实例

算法实现的一般流程

• 数据加载

  加载后的数据用于训练分类器，得出单层决策树；

• 寻找数据集上的最佳单层决策树（弱分类器）

  因为我们不知道数据集中的哪一维（特征）的数据对分类有显著的效果，也不知道在该特征下哪个值（阈值）可以将数据分为两类。所以我们对每一维度的特征数据使用小步长迭代的方法，去找到一个最佳的分类阈值。那么怎样才算找到呢？我们有一个评判标准，即当分类误差（每一组特征向量赋予权值，在每个阈值的改变后得到的分类结果，找到其中分类错误的数据集，计算出其权值的和记为误差）最小时所记录的阈值、分类维度、分类方式（比阈值大为正或者比阈值小为正）即为单层决策树的构造结果。

• Adaboost算法使用迭代寻找不同权值alpha下的单层决策树（弱分类器）

  每次迭代后更新样本权重，该权重的计算上面已经给出类数学模型，累加每次由上一次单层决策树的构造结果的最小误差计算出的权重alpha与分类结果的乘积aggClassEst；再根据AdaBoost误差计算，当误差值为0时退出循环，此时的aggClassEst序列值的符号即为分类结果；

• 对输入的数据进行预测

  利用训练集的多次迭代得到的多个弱分类器的结果对输入数据执行分类，累加每次的权重alpha与弱分类器的结果的积得到aggClassEst预测结果，取预测结果的符号作为最终的预测情况。

• 调度函数

  起调度的作用；

算法的实现

• 用于训练弱分类器的数据的构造

def loadSimpData():
    """
    创建单层决策树的数据集
    Parameters:
        无
    Returns:
        dataMat - 数据矩阵
        classLabels - 数据标签
    """
    datMat = np.matrix([[ 1. ,  2.1],
        [ 1.5,  1.6],
        [ 1.3,  1. ],
        [ 1. ,  1. ],
        [ 2. ,  1. ]])
    classLabels = [1.0, 1.0, -1.0, -1.0, 1.0]
    return datMat,classLabels

• 单层决策树的构造

def stumpClassify(dataMatrix,dimen,threshVal,threshIneq):
    """
    单层决策树分类函数
    Parameters:
        dataMatrix - 数据矩阵
        dimen - 第dimen列，也就是第几个特征
        threshVal - 阈值
        threshIneq - 标志
    Returns:
        retArray - 分类结果
    """
    retArray = np.ones((np.shape(dataMatrix)[0],1))                #初始化retArray为1
    if threshIneq == 'lt':
        retArray[dataMatrix[:,dimen] <= threshVal] = -1.0   #如果小于阈值,则赋值为-1
    else:
        retArray[dataMatrix[:,dimen] > threshVal] = -1.0    #如果大于阈值,则赋值为-1
    return retArray
    
def buildStump(dataArr,classLabels,D):
    """
    找到数据集上最佳的单层决策树
    Parameters:
        dataArr - 数据矩阵
        classLabels - 数据标签
        D - 样本权重
    Returns:
        bestStump - 最佳单层决策树信息
        minError - 最小误差
        bestClasEst - 最佳的分类结果
    """
    dataMatrix = np.mat(dataArr); labelMat = np.mat(classLabels).T
    m,n = np.shape(dataMatrix)
    numSteps = 10.0; bestStump = {}; bestClasEst = np.mat(np.zeros((m,1)))
    minError = float('inf')     #最小误差初始化为正无穷大
    for i in range(n):                                                            #遍历所有特征
        rangeMin = dataMatrix[:,i].min(); rangeMax = dataMatrix[:,i].max()        #找到特征中最小的值和最大值
        stepSize = (rangeMax - rangeMin) / numSteps                                #计算步长
        for j in range(-1, int(numSteps) + 1):                                     
            for inequal in ['lt', 'gt']: #大于和小于的情况，均遍历。lt:less than，gt:greater than
                threshVal = (rangeMin + float(j) * stepSize) #计算阈值
                predictedVals = stumpClassify(dataMatrix, i, threshVal, inequal)#计算分类结果
                errArr = np.mat(np.ones((m,1)))             #初始化误差矩阵
                errArr[predictedVals == labelMat] = 0           #分类正确的,赋值为0
                weightedError = D.T * errArr    #计算误差
                # print("split: dim %d, thresh %.2f, thresh ineqal: %s, the weighted error is %.3f" % (i, threshVal, inequal, weightedError))
                if weightedError < minError:                #找到误差最小的分类方式
                    minError = weightedError
                    bestClasEst = predictedVals.copy()
                    bestStump['dim'] = i
                    bestStump['thresh'] = threshVal
                    bestStump['ineq'] = inequal
    return bestStump, minError, bestClasEst

• 使用迭代的方式训练弱分类器

ef adaBoostTrainDS(dataArr, classLabels, numIt = 40):
    """
    使用AdaBoost算法提升弱分类器性能
    Parameters:
        dataArr - 数据矩阵
        classLabels - 数据标签
        numIt - 最大迭代次数
    Returns:
        weakClassArr - 训练好的分类器
        aggClassEst - 类别估计累计值
    """
    weakClassArr = []
    m = np.shape(dataArr)[0]
    D = np.mat(np.ones((m, 1)) / m)                                            #初始化权重
    aggClassEst = np.mat(np.zeros((m,1)))
    for i in range(numIt):
        bestStump, error, classEst = buildStump(dataArr, classLabels, D)     #构建单层决策树
        # print("D:",D.T)
        alpha = float(0.5 * np.log((1.0 - error) / max(error, 1e-16)))         #计算弱学习算法权重alpha,使error不等于0,因为分母不能为0
        bestStump['alpha'] = alpha      #存储弱学习算法权重
        weakClassArr.append(bestStump)  #存储单层决策树
        # print("classEst: ", classEst.T)
        expon = np.multiply(-1 * alpha * np.mat(classLabels).T, classEst)     #计算e的指数项
        D = np.multiply(D, np.exp(expon))                                      
        D = D / D.sum()     #根据样本权重公式，更新样本权重
        #计算AdaBoost误差，当误差为0的时候，退出循环
        aggClassEst += alpha * classEst     #计算类别估计累计值                                
        # print("aggClassEst: ", aggClassEst.T)
        aggErrors = np.multiply(np.sign(aggClassEst) != np.mat(classLabels).T, np.ones((m,1)))     #计算误差
        errorRate = aggErrors.sum() / m
        # print("total error: ", errorRate)
        if errorRate == 0.0: break      #误差为0，退出循环
    return weakClassArr, aggClassEst

• 执行分类预测

def adaClassify(datToClass,classifierArr):
    """
    AdaBoost分类函数
    Parameters:
        datToClass - 待分类样例
        classifierArr - 训练好的分类器
    Returns:
        分类结果
    """
    dataMatrix = np.mat(datToClass)
    m = np.shape(dataMatrix)[0]
    aggClassEst = np.mat(np.zeros((m,1)))
    for i in range(len(classifierArr)): #遍历所有分类器，进行分类
        classEst = stumpClassify(dataMatrix, classifierArr[i]['dim'], classifierArr[i]['thresh'], classifierArr[i]['ineq'])            
        aggClassEst += classifierArr[i]['alpha'] * classEst
        print(aggClassEst)
    return np.sign(aggClassEst)

• 调度函数

if __name__ == '__main__':
    dataArr,classLabels = loadSimpData()
    weakClassArr, aggClassEst = adaBoostTrainDS(dataArr, classLabels)
    print(adaClassify([[0,0],[5,5]], weakClassArr))

预测结果

八、线性回归算法

8.1 算法简述

在统计学中，线性回归（Linear Regression）是利用称为线性回归方程的最小平方函数对一个或多个自变量和因变量之间关系进行建模的一种回归分析。这种函数是一个或多个称为回归系数的模型参数的线性组合（自变量都是一次方）。只有一个自变量的情况称为简单回归，大于一个自变量情况的叫做多元回归

8.2 算法的数学模型

假定输入数据存放在矩阵X中，结果存放在向量y中：

而回归系数存放在向量w中：

所以线性方程可表示为：y=wx+e,其中e为误差。

方程的求解

• 最小二乘法

  

  用矩阵表示还可以写做：

  

  我们要先明确一个目的：找到w，使平方误差和最小。因为我们认为平方误差和越小，说明线性回归拟合效果越好。
  现在，我们用矩阵表示的平方误差和对w进行求导：

  

  令上述公式等于0，得到：

  

  w上方的小标记表示，这是当前可以估计出的w的最优解。

  上面用到了矩阵的逆，所以我们假设所给的数据矩阵是可逆的。

• 局部加权法

  线性回归的一个问题是有可能出现欠拟合现象，因为它求的是具有小均方误差的无偏估 计。显而易见，如果模型欠拟合将不能取得好的预测效果。所以有些方法允许在估计中引入一 些偏差，从而降低预测的均方误差。
  其中的一个方法是局部加权线性回归（Locally Weighted Linear Regression，LWLR）。在该方法中，我们给待预测点附近的每个点赋予一定的权重。

  该算法对回归系数W的估计方程如下：

  

  其中W是一个矩阵，这个公式跟我们上面推导的公式的区别就在于W，它用来给每个点赋予权重。

  LWLR使用"核"（与支持向量机中的核类似）来对附近的点赋予更高的权重。核的类型可以自由选择，最常用的核就是高斯核，高斯核对应的权重如下：

  

  这样我们就可以通过改变k的值，调节回归效果。

• 梯度下降法

• 岭回归（ridge regression)和缩减方法

8.4 算法的运用实例

训练一下数据，得到其线性回归模型；

1.000000	0.067732	3.176513
1.000000	0.427810	3.816464
1.000000	0.995731	4.550095
1.000000	0.738336	4.256571
1.000000	0.981083	4.560815
1.000000	0.526171	3.929515
1.000000	0.378887	3.526170
1.000000	0.033859	3.156393
1.000000	0.132791	3.110301
1.000000	0.138306	3.149813
...

第一列都为1.0，即x0。第二列为x1，即x轴数据。第三列为x2，即y轴数据。

算法的实现

• 首先绘制下数据，看下数据分布

  import matplotlib.pyplot as plt
  import numpy as np
   
  def loadDataSet(fileName):
      """
      函数说明:加载数据
      Parameters:
          fileName - 文件名
      Returns:
          xArr - x数据集
          yArr - y数据集
      """
   
      numFeat = len(open(fileName).readline().split('\t')) - 1
      xArr = []; yArr = []
      fr = open(fileName)
      for line in fr.readlines():
          lineArr =[]
          curLine = line.strip().split('\t')
          for i in range(numFeat):
              lineArr.append(float(curLine[i]))
          xArr.append(lineArr)
          yArr.append(float(curLine[-1]))
      return xArr, yArr
   
  def plotDataSet():
      """
      函数说明:绘制数据集
      Parameters:
          无
      Returns:
          无
      """
      xArr, yArr = loadDataSet('ex0.txt')     #加载数据集
      n = len(xArr)                           #数据个数
      xcord = []; ycord = []                  #样本点
      for i in range(n):                                      
          xcord.append(xArr[i][1]); ycord.append(yArr[i])  #样本点
      fig = plt.figure()
      ax = fig.add_subplot(111)           #添加subplot
      ax.scatter(xcord, ycord, s = 20, c = 'blue',alpha = .5) #绘制样本点
      plt.title('DataSet')                #绘制title
      plt.xlabel('X')
      plt.show()
   
  if __name__ == '__main__':
      plotDataSet()
  

运行结果：

• 最小二乘法实现线性回归

  def standRegres(xArr,yArr):
      """
      函数说明:计算回归系数w
      Parameters:
          xArr - x数据集
          yArr - y数据集
      Returns:
          ws - 回归系数
      """
      xMat = np.mat(xArr); yMat = np.mat(yArr).T
      xTx = xMat.T * xMat                            #根据文中推导的公示计算回归系数
      if np.linalg.det(xTx) == 0.0:
          print("矩阵为奇异矩阵,不能求逆")
          return
      ws = xTx.I * (xMat.T*yMat)
      return ws
  

• 判断拟合效果

  线性拟合时我们一般可以采用相关系数来判断拟合的效果，数学模型如下:

  

  # 直接调用np.corrcoef()函数
  corr = np.corrcoef(yHat.T, yMat)
  

• 局部加权线性回归

  数学模型上面已经给出，这里是代码实现：

  ef lwlr(testPoint, xArr, yArr, k = 1.0):
  """
  函数说明:使用局部加权线性回归计算回归系数w
  Parameters:
      testPoint - 测试样本点
      xArr - x数据集
      yArr - y数据集
      k - 高斯核的k,自定义参数
  Returns:
      ws - 回归系数
  """
  xMat = np.mat(xArr); yMat = np.mat(yArr).T
  m = np.shape(xMat)[0]
  weights = np.mat(np.eye((m)))                                        #创建权重对角矩阵
  for j in range(m):                                                  #遍历数据集计算每个样本的权重
      diffMat = testPoint - xMat[j, :]                                 
      weights[j, j] = np.exp(diffMat * diffMat.T/(-2.0 * k**2))
  xTx = xMat.T * (weights * xMat)                                        
  if np.linalg.det(xTx) == 0.0:
      print("矩阵为奇异矩阵,不能求逆")
      return
  ws = xTx.I * (xMat.T * (weights * yMat))                            #计算回归系数
  return ws
  

• 局部加权线性回归的测试

  def lwlrTest(testArr, xArr, yArr, k=1.0):  
  """
  函数说明:局部加权线性回归测试
  Parameters:
      testArr - 测试数据集
      xArr - x数据集
      yArr - y数据集
      k - 高斯核的k,自定义参数
  Returns:
      ws - 回归系数
  """
  m = np.shape(testArr)[0]                                            #计算测试数据集大小
  yHat = np.zeros(m)    
  for i in range(m):                                                    #对每个样本点进行预测
      yHat[i] = testArr[i]*lwlr(testArr[i],xArr,yArr,k)
  return yHat
  

• 回归结果可视化

  def plotlwlrRegression():
  """
  函数说明:绘制多条局部加权回归曲线
  Parameters:
      无
  Returns:
      无
  """
  font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
  xArr, yArr = loadDataSet('ex0.txt')                                    #加载数据集
  yHat_1 = lwlrTest(xArr, xArr, yArr, 1.0)                            #根据局部加权线性回归计算yHat
  yHat_2 = lwlrTest(xArr, xArr, yArr, 0.01)                            #根据局部加权线性回归计算yHat
  yHat_3 = lwlrTest(xArr, xArr, yArr, 0.003)                            #根据局部加权线性回归计算yHat
  xMat = np.mat(xArr)                                                    #创建xMat矩阵
  yMat = np.mat(yArr)                                                    #创建yMat矩阵
  srtInd = xMat[:, 1].argsort(0)                                        #排序，返回索引值
  xSort = xMat[srtInd][:,0,:]
  fig, axs = plt.subplots(nrows=3, ncols=1,sharex=False, sharey=False, figsize=(10,8))                                        
  axs[0].plot(xSort[:, 1], yHat_1[srtInd], c = 'red')                        #绘制回归曲线
  axs[1].plot(xSort[:, 1], yHat_2[srtInd], c = 'red')                        #绘制回归曲线
  axs[2].plot(xSort[:, 1], yHat_3[srtInd], c = 'red')                        #绘制回归曲线
  axs[0].scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = 'blue', alpha = .5)                #绘制样本点
  axs[1].scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = 'blue', alpha = .5)                #绘制样本点
  axs[2].scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = 'blue', alpha = .5)                #绘制样本点
  #设置标题,x轴label,y轴label
  axs0_title_text = axs[0].set_title(u'局部加权回归曲线,k=1.0',FontProperties=font)
  axs1_title_text = axs[1].set_title(u'局部加权回归曲线,k=0.01',FontProperties=font)
  axs2_title_text = axs[2].set_title(u'局部加权回归曲线,k=0.003',FontProperties=font)
  plt.setp(axs0_title_text, size=8, weight='bold', color='red')  
  plt.setp(axs1_title_text, size=8, weight='bold', color='red')  
  plt.setp(axs2_title_text, size=8, weight='bold', color='red')  
  plt.xlabel('X')
  plt.show()
  

  可视化结果：

  

  可以看到，当k越小，拟合效果越好。但是当k过小，会出现过拟合的情况，例如k等于0.003的时候。

九、随机森林算法

9.1 算法简述

在决策树的基础上，结合多个分类器组合的思想，由多个决策树生成随机森林。随机森林的核心思想就是，对训练集进行重采样，组成多个训练子集，每个子集生成一个决策树，所有决策树通过投票的方式进行决策，组成随机森林。

9.2 算法的数学模型

bootstrap sampling（自助采样法）

随机有放回采样获取训练样本每一轮采样中，样本x被抽到的概率为1/m，因此，在m轮抽样后，该样本仍未被抽取到的概率为：

因此，原始数据集中有大约36.8%(1/3)的样本没有被抽到，这些样本就可以用于由新数据集D’训练得到的学习器的“包外估计”

Bagging（套袋法）

Bagging就是通过在每轮迭代中，首先通过自助采样法生成不同的训练集，之后基于这个“新”训练集对基学习器进行训练。最后，集成学习器的预测结果将由T轮迭代后生成的T个基学习器，进行简单的投票法

bagging的算法过程如下：

• 从原始样本集中使用Bootstraping方法随机抽取n个训练样本，共进行k轮抽取，得到k个训练集。（k个训练集之间相互独立，元素可以有重复）
• 对于k个训练集，我们训练k个模型（这k个模型可以根据具体问题而定，比如决策树，knn等）
• 对于分类问题：由投票表决产生分类结果；对于回归问题：由k个模型预测结果的均值作为最后预测结果。（所有模型的重要性相同）

决策树

常用的决策树算法有ID3，C4.5，CART三种。3种算法的模型构建思想都十分类似，只是采用了不同的指标。决策树模型的构建过程大致如下：

• ID3，C4.5决策树的生成

  输入：训练集D，特征集A，阈值eps 输出：决策树T

  • 若D中所有样本属于同一类Ck，则T为单节点树，将类Ck作为该结点的类标记，返回T
  • 若A为空集，即没有特征作为划分依据，则T为单节点树，并将D中实例数最大的类Ck作为该结点的类标记，返回T
  • 否则，计算A中各特征对D的信息增益(ID3)或者信息增益比(C4.5)，选择信息增益最大的特征Ag
  • 若Ag的信息增益（比）小于阈值eps，则置T为单节点树，并将D中实例数最大的类Ck作为该结点的类标记，返回T
  • 否则，依照特征Ag将D划分为若干非空子集Di，将Di中实例数最大的类作为标记，构建子节点，由结点及其子节点构成树T，返回T
  • 对第i个子节点，以Di为训练集，以A-{Ag}为特征集，递归地调用1~5，得到子树Ti，返回Ti

• CART决策树的生成

• 这里只简单介绍下CART与ID3和C4.5的区别。

  • CART树是二叉树，而ID3和C4.5可以是多叉树
  • CART在生成子树时，是选择一个特征一个取值作为切分点，生成两个子树
  • 选择特征和切分点的依据是基尼指数，选择基尼指数最小的特征及切分点生成子树

算法的运用实例

算法实现的一般流程

• 从原始训练集中使用Bootstraping方法随机有放回采样选出m个样本，共进行n_tree次采样，生成n_tree个训练集
• 对于n_tree个训练集，我们分别训练n_tree个决策树模型
• 对于单个决策树模型，假设训练样本特征的个数为n，那么每次分裂时根据信息增益/信息增益比/基尼指数选择最好的特征进行分裂
• 每棵树都一直这样分裂下去，直到该节点的所有训练样例都属于同一类。在决策树的分裂过程中不需要剪枝
• 将生成的多棵决策树组成随机森林。对于分类问题，按多棵树分类器投票决定最终分类结果；对于回归问题，由多棵树预测值的均值决定最终预测结果

算法的实现（以CART分类树作为弱分类器）

• 样本的选取和特征的选取

input:  data(list):原始数据集
        k(int):选择特征的个数
output: data_samples(list):被选择出来的样本
        feature(list):被选择的特征index
"""
def choose_sample(data,k):
    m,n = np.shape(data) # 样本的个数和样本特征的个数

    # 1、选择出k个特征的index
    feature = []
    for j in range(k):
        feature.append(rd.randint(0,n-2))  # n-1列是标签
    # 2、选择出m个样本的index
    index = []
    for i in range(m):
        index.append(rd.randint(0,m-1))
    # 3、从data中选择出m个样本的k个特征，组成数据集data_samples
    data_samples = []
    for i in range(m):
        data_tmp = []
        for fea in feature:
            data_tmp.append(data[index[i]][fea])  # 添加i行fea列
        data_tmp.append(data[index[i]][-1]) # 添加标签值
        data_samples.append(data_tmp)
    return data_samples,feature

• 构建CART分类树

from math import pow

"""
树的节点类
"""
class node:
    def __init__(self,fea=-1,value=None,results=None,right=None,left=None):
        self.fea = fea
        self.value =value
        self.results = results
        self.right = right
        self.left = left

"""
根据特征fea中的值value将数据集data划分成左右子树
input:  data(list):数据集
        fea(int):待分割特征的索引
        value(float):待分割的特征的具体值
output: (set1,set2)(tuple):分割后的左右子树
"""
def split_tree(data,fea,value):
    set_1 = []
    set_2 = []
    for x in data:
        if x[fea] >= value:  # 大于等于一般是在左子树
            set_1.append(x)
        else:
            set_2.append(x)
    return (set_1,set_2)

"""
统计数据集中不同的类标签label的个数
input:  data(list):原始数据集
output: label_uniq_cnt(int):样本中的标签的个数
"""
def label_uniq_cnt(data):
    label_uniq_cnt = {}  # 字典

    for x in data:
        label = x[len(x) -1]  
        if label not in label_uniq_cnt:  # 取得每一个样本的类标签label
            label_uniq_cnt[label] = 0
        label_uniq_cnt[label] += 1
    return label_uniq_cnt

"""
计算给定数据的gini指数
    input:  data(list):树中
    output: gini(float):Gini指数
"""
def cal_gini_index(data):
    total_sample = len(data)  # 样本的总个数 
    if len(data) == 0:
        return 0
    label_counts = label_uniq_cnt(data)  # 统计数据集中不同标签的个数

    # 计算数据集的Gini指数
    gini = 0
    for label in label_counts:
        gini = gini + pow(label_counts[label],2)

    gini = 1 - float(gini)/pow(total_sample,2)
    return gini

"""
构建树
    input:  data(list):训练样本
    output: node:树的根结点
"""
def build_tree(data):
    # 构建决策树，函数返回该决策树的根节点
    if len(data) == 0:
        return node()

    # 1、计算当前的Gini指数
    currentGini = cal_gini_index(data)

    bestGain = 0.0
    bestCriteria = None # 存储最佳切分属性以及最佳切分点
    bestSets = None

    feature_num = len(data[0]) -1  # 样本中特征的个数
    # 2、找到最好的划分
    for fea in range(0,feature_num):
        # 2.1、取得fea特征处所有可能的取值
        feature_values = {}  # 在fea位置处可能的取值
        for sample in data:  # 对每一个样本
            feature_values[sample[fea]] = 1  # 存储特征fea处所有可能的取值

        # 2.2、针对每一个可能的取值，尝试将数据集划分，并计算Gini指数  
        for value in list(feature_values.keys()):  # 遍历该属性的所有切分点
            # 2.2.1、 根据fea特征中的值value将数据集划分成左右子树
            (set_1,set_2) = split_tree(data,fea,value)
            # 2.2.2、计算当前的Gini指数
            nowGini = float(len(set_1)*cal_gini_index(set_1)+len(set_2)*cal_gini_index(set_2))/len(data)
            # 2.2.3、计算Gini指数的增加量
            gain = currentGini - nowGini
            # 2.2.4、判断此划分是否比当前的划分更好
            if gain > bestGain and len(set_1) > 0 and len(set_2) > 0:
                bestGain = gain
                bestCriteria = (fea,value)
                bestSets = (set_1,set_2)

    # 3、判断划分是否结束
    if bestGain > 0:
        right = build_tree(bestSets[0])
        left = build_tree(bestSets[1])
        return node(fea=bestCriteria[0],value=bestCriteria[1],right=right,left=left)
    else:
        return node(results=label_uniq_cnt(data))  # 返回当前的类别标签作为最终的类别标签

• 用多个决策树构建随机深林

"""
构建随机森林
input:  data_train(list):训练数据
        trees_num(int):分类树的个数
output: trees_result(list):每一棵树的最好划分
        trees_feature(list):每一棵树中对原始特征的选择
"""
def random_forest_training(data_train,tress_num):
    trees_result = []  # 构建好每一棵树的最好划分
    trees_feature = []
    n = np.shape(data_train)[1] # 样本的维数

    if n > 2:
        k = int(log(n-1,2)) + 1   # 设置特征的个数
    else:
        k = 1
    # 开始构建每一棵树
    for i in range(tress_num):
        # 1、随机选择m个样本, k个特征
        data_samples,feature = choose_sample(data_train,k)
        # 2、构建每一棵分类树
        tree = build_tree(data_samples)
        # 3、保存训练好的分类树
        trees_result.append(tree)
        # 4、保存好该分类树使用到的特征
        trees_feature.append(feature)

    return trees_result,trees_feature

• 对每一个样本的进行预测（测试）

"""
对每一个样本sample进行预测
    input:  sample(list):需要预测的样本
            tree(类):构建好的分类树
    output: tree.results:所属的类别
"""
def predict(sample,tree):
    # 1、只是树根
    if tree.results != None:
        return tree.results
    else:
        # 2、有左右子树
        val_sample = sample[tree.fea]
        branch = None
        if val_sample >= tree.value:
            branch = tree.right
        else:
            branch = tree.left
        return predict(sample, branch)
        
"""
选择特征
input:  data_train(list):训练数据集
        feature(list):要选择的特征
output: data(list):选择出来的数据集
"""
def split_data(data_train,feature):
    m = np.shape(data_train)[0]
    data = []

    for i in range(m):
        data_x_tmp = []
        for x in feature:
            data_x_tmp.append(data_train[i][x])
        data_x_tmp.append(data_train[i][-1])
        data.append(data_x_tmp)
    return data

"""
trees_fiture:保存了每一颗分类树中随机选择的特征
"""
def get_predict(trees_result,trees_fiture,data_train):
    m_tree = len(trees_result) # 有多少棵树
    m = np.shape(data_train)[0]  # 训练集的数目

    result = []
    for i in range(m_tree):
        clf = trees_result[i]
        feature = trees_fiture[i]
        data = split_data(data_train,feature)
        result_i = []
        for i in range(m):
            result_i.append((list(predict(data[i][0:-1],clf).keys()))[0])
        result.append(result_i)
    final_predict = np.sum(result,axis=0)
    return final_predict

• 计算错误率

def cal_correct_rate(data_train,final_predict):
    m = len(final_predict)
    corr = 0.0
    for i in range(m):
        if data_train[i][-1] * final_predict[i] > 0:
            corr += 1
    return corr/m

• 入口函数

if __name__ == "__Main__":
    # 训练数据集
    file_path_train = 'filepath1'
    
    # 需要预测的数据
    file_path_predict = 'filepath2'
    # 1、导入数据
    print("----------- 1、load data -----------")
    data_train = load_data(file_path_train)
    # 2、训练random_forest模型
    print("----------- 2、random forest training ------------")
    trees_result,trees_feature = random_forest_training(data_train,1)
    # 3、得到训练的准确性
    print("------------ 3、get prediction correct rate ------------")
    result = get_predict(trees_result,trees_feature,data_train)
    corr_rate = cal_correct_rate(data_train,result)
    print("\t------correct rate: ", corr_rate) 
    
    # 4、加载测试数据
    print("----------- 1、load data -----------")
    data_test = load_data(file_path_predict)
    
    # 5、进行预测
    print("------------ 2、get prediction correct rate ------------")
    result = get_predict(trees_result,trees_feature,data_test)

十、神经网络/人工神经网络

10.1 算法简述

10.2 算法的数学模型

10.4 算法的运用实例

算法实现的一般流程

算法的实现