机器学习之Logistic回归解析及实例应用

参考博客：https://techlog.cn/article/list/10183274https://blog.csdn.net/u012762410/article/details/78912667前言：学习《机器学习实战》这本书到现在，这一章节算是数学理论较多的，也很高兴自己能通过搜索资料和学习他人博客推导出运算公式并了解代码含义，对自己而言也是一个小的突破，继续写下机器学习博客记录。Lo

疯狂奔跑的少年

2264人浏览 · 2020-08-29 14:27:03

疯狂奔跑的少年 · 2020-08-29 14:27:03 发布

前言：学习《机器学习实战》这本书到现在，这一章节算是数学理论较多的，也很高兴自己能通过搜索资料和学习他人博客推导出运算公式并了解代码含义，对自己而言也是一个小的突破，继续写下机器学习博客记录。

Logistic回归

回归：假设有一些数据点，我们用一条直线对这些点进行拟合（该线称为最佳拟合直线），这个拟合的过程就称为回归。

拟合：形象上说就是把平面上一系列的点，用一条光滑的曲线连接起来，因为这条曲线有无数种可能，从而有各种拟合法，拟合的曲线一般可以用函数表示，根据这个函数的不同有不同的拟合名字。

利用Logistic回归进行分类的主要思想：根据现有的数据对分类边界线建立回归公式，以次进行分类。

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

Logistic回归的一般过程：

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

基于Logistic回归和Sigmoid函数的分类：
Logistic回归要进行分类，就要构建分类器，找到最佳的拟合，即要找到最佳的拟合参数，我们所想到的函数应该是能接受所有的输入然后预测出其类别。比如在两个类的情况下，上述函数输出0或1。

Sigmoid函数：
在上述情况下，我们可以使用Simoid函数来解决，同样，它对于Logistic回归分类也是十分重要的。它的计算公式如下： $\sigma(z)=\frac{1}{1+e^{-z}}$ 我们可以看看它的函数图像在这里插入图片描述
这是Sigmoid函数在不同坐标尺度下的两条曲线图，我们可以看到，当x为0时，Sigmoid函数值为0.5，随着x的增大，其值也不断趋近于1，相反，当x减小，其值不断趋近于0，如果横坐标足够大，比如下图，函数看起来就很像一个阶跃函数。

因此实现Logistic回归分类器，我们可以每个特征都乘上一个回归系数，然后所有值相加，将这个总和作为z代入Sigmoid函数中，进而得到一个范围在0-1的数值，将大于0.5的数据分入1类，小于0.5的即分为0类，所以，Logistic回归也可以看作是一个概率估计。

但是，了解了这个函数还不够我们构架Logistic回归分类器，下面我们就要进行一些数学公式推导。

几率函数和Logit函数：
一个事件发生的概率是p，那么其不发生的概率就是1-p，那么p/(1-p)就是事件发生的几率，而 $f (P) = P / (1 - P)$ 就是几率函数。比如投篮，投中的概率是90%，不中的概率是10%，f(0.9)=9，即投10球的情况下，其可以中9个。那么对几率函数取对数就是Logit函数 $Logit(P)=log\frac{P}{1-P}$ 如果我们将上诉例子看成是条件是投10个球，投中概率是90%时得到中9个的特征的结果，但我们需要解决的是符合一定特征的样本其属于某个类别的概率，就需要对上面的函数求反，得到Sigmoid函数，也称为Logistic函数 $\sigma(z)=\frac{1}{1+e^{-z}}$
代价函数：
假设我们现在有个样本数据x，其有n个特征值x₀,x₁,…x_n，我们将其与系数θ相乘作为Sigmoid函数的输入值 $Z=\left[ \begin{matrix} \theta_0 & \theta_2 &... & \theta_n \\ \end{matrix} \right] \left[ \begin{matrix} x_0 \\ x_2 \\ ...\\ x_n \\ \end{matrix} \right] =\theta^Tx$
进而我们可以得到 $h_\theta(x)=Sigmoid(\theta^Tx)=\frac{1}{1+e^{-\theta^Tx}}$
我们已经知道Sigmoid函数的阈值就是在0-1之间，我们只要有合适的系数的列向量和样本数据集的列向量，就可以计算出该样本对应的概率值，从而对其进行分类（上述提到以0.5作为分界）。那么我们所要做的就是找出最佳的系数。

解决二分类的问题，我们可以这样写出假设： $P(Y=1|x;\theta)=h_\theta(x)\\P(Y=0|x;\theta)=1-h_\theta(x)$
这两个公式的含义也很好理解，我们要计算最终的分类是属于1还是0，那么在样本数据和系数已知的条件下，样本集属于1类（即y=1)的概率和属于0类的概率，因为y仅有两种取值，所以我们将两个公式合并 $P(Y|x;\theta)=(h_\theta(x))^Y(1-h_\theta(x))^{1-Y}$ 也就是 $Cost(h_\theta(x),Y)=(h_\theta(x))^Y(1-h_\theta(x))^{1-Y}$ 这也就是代价函数。

为了后面的简化问题，我们可以对其进行求对数，即 $Cost(h_\theta(x),Y)=Ylog(h_\theta(x))+(1-Y)log(1-h_\theta(x))$ 我们通过求这个代价函数即是求样本所属的类别概率，我们当然希望这个概率能更准确的反应其类别，所以自然希望它的该概率值越大，相应的，也就希望系数也就是最佳的，而找系数θ的过程就是极大似然估计。极大似然估计的目的也就是利用已知的样本结果，反推最大概率导致这样结果的参数值，假定样本间的特征相互独立，那么整个样本集生成的概率即为所有样本生成概率的乘积 $L(\theta|x,Y)=\prod_{i=1}^mY^{(i)}log(h_\theta(x))+(1-Y^{(i)})log(1-h_\theta(x^{(i)}))$ 为了和上面的函数保持一致同样换成对数似然函数： $j(\theta)=L(\theta|x,Y)=\sum_{i=1}^mY^{(i)}log(h_\theta(x))+(1-Y^{(i)})log(1-h_\theta(x^{(i)}))$ 其中m是样本总数，Y⁽ⁱ⁾是第i个样本的类别，x⁽ⁱ⁾是第i个样本，其中θ和x⁽ⁱ⁾是多维向量。

梯度上升优化算法：
我们要求取J(θ)最大的θ值，所以我们后面采用的是梯度上升法。它的思想是：要找到某函数的最大值，最好的方法就是沿着该函数的梯度方向去探寻。

这就好比是爬山，当不断向上快到山顶的时候，我们也就越接近极值。在这里插入图片描述
就如图函数，每次上升，就与之前的值比较，若是大于前者说明处于上升的阶段，若是小于上次的值，说明处于下降的阶段，那么前者就是极值点了。与此同时，我们要控制步长，步子太小，耗费的时间会很多，若是步子太大，很可能会错过极值。由于函数的导数随着接近极值点是不断减小的（在极值点为0），所以只要步长正比于函数的导数，步子就可以在接近极值点的时候越来越小，定义上升函数： $\theta_{j+1}=\theta_j +\alpha\frac{\partial j_{(\theta)}}{\theta_j}$ 其中α是步长系数。该公式将一直被迭代执行，直到达到某个条件为止，比如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。

我们知道上升公式： $\theta_{j+1}=\theta_j +\alpha\frac{\partial j_{(\theta)}}{\theta_j}$
所以我们可以求j_(θ)的偏导，继而求得极大值。我们前面已经求得： $j(\theta)=L(\theta|x,Y)=\sum_{i=1}^mY^{(i)}log(h_\theta(x))+(1-Y^{(i)})log(1-h_\theta(x^{(i)}))$ ,并且 $h_\theta(x)=Sigmoid(\theta^Tx)=\frac{1}{1+e^{-\theta^Tx}}$ 所以对于 $\alpha\frac{\partial j_{(\theta)}}{\theta_j}$ 的导数，我们可以分为三部分来求解，即 $\frac{\partial j_{(\theta)}}{\theta_j}=\frac{\partial j_{(\theta)}}{\partial Sigmoid(\theta^Tx)}*\frac{\partial Sigmoid(\theta^Tx)}{\partial \theta^Tx}*\frac{\partial \theta^Tx}{\partial \theta_j}$
其中第一部分的求导 $\frac{\partial j_{(\theta)}}{\partial Sigmoid(\theta^Tx)}=Y*\frac{1}{Sigmoid(\theta^Tx)}+(Y-1)*\frac{1}{1-Sigmoid(\theta^Tx)}$ 第二部分的求导是 $\frac{\partial Sigmoid(\theta^Tx)}{\partial \theta^Tx}=\frac{1}{1+e^{-\theta^Tx}}*(1-\frac{1}{1+e^{-\theta^Tx}})=Sigmoid(\theta^Tx)*(1-Sigmoid(\theta^Tx))$ 而第三部分的求导是 $\frac{\partial \theta^Tx}{\partial \theta_j}=\frac{\partial (\theta_1x_1+\theta_2x_2+...+\theta_nx_n)}{\partial \theta_j}=x_j$ 由此三部分的导数已经求解完成，我们可以整合算式得出最后的结果就是 $\frac{\partial j_{(\theta)}}{\theta_j}=(Y-h_{(\theta)}x)x_j$ 所以最后的梯度上升迭代公式就是 $\theta_j:=\theta_j+\alpha*\sum_{i=1}^n(Y^{(i)}-h_{(\theta)}x^{(i)})x_j^{(i)}$

准备数据：
这次的数据同样与之前的一样下载，在testSet.txt文本中共有100个数据，每个数据包含两个数值型的特征X1和X2，与此同时还有其所属的类别，我们通过这个数据集来找出其最佳回归系数。所以我们可以编写代码先获取准备数据：

from numpy import *
def loadDataSet():
    dataMat = []
    labelMat = []
    fr = open('testSet.txt')
    #读取文本文件
    for line in fr.readlines():
        #读取全部行内容并且逐行遍历
        lineArr = line.strip().split()
        #对单行进行处理，伤处空格换行符并且切分
        dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])
        #将文本中的数据的前两个x1,x2和1.0作为x0一起放入列表中
        labelMat.append(int(lineArr[2]))
        #将标签放入标签列表
    return dataMat,labelMat
if __name__=='__main__':
	dataArr,labelMat = loadDataSet()
	print(dataArr)
    print(labelMat)

在这里插入图片描述
这样我们就获取了文本中的数据，并且将其切分成我们所需要的形式。

对于上述代码理解自然没有问题，但是其中数据列表中添加元素的时候，为什么除了数据本身具有的两个特征值，还加入了一个1.0作为X0，这个我们后面解释。

训练算法：使用梯度上升找到最佳参数
我们前面提到Sigmoid函数的输入我们定义为Z，且Z为最佳系数与特征值的向量积，再结合我们之前推导出的梯度上升迭代公式，将迭代公式矢量化 $\theta_j:=\theta_j+\alpha*(Y-Sigmoid(\theta x))x^T$ 其中Y是向量。我们就可以编写代码：

from numpy import *
def loadDataSet():
    dataMat = []
    labelMat = []
    fr = open('testSet.txt')
    #读取文本文件
    for line in fr.readlines():
        #读取全部行内容并且逐行遍历
        lineArr = line.strip().split()
        #对单行进行处理，伤处空格换行符并且切分
        dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])
        #将文本中的数据的前两个x1,x2和1.0作为x0一起放入列表中
        labelMat.append(int(lineArr[2]))
        #将标签放入标签列表
    return dataMat,labelMat

def sigmoid(inX):#sigmoid函数计算
    return 1.0/(1+exp(-inX))

def gradAscent(dataMatIn,classLabels):#梯度上升优化算法求最佳回归系数
    dataMatrix = mat(dataMatIn)
    #将训练特征数据转为numpy矩阵
    labelMat = mat(classLabels).transpose()
    #将标签列表转为numpy矩阵并将其转置
    m,n = shape(dataMatrix)
    #获得矩阵的行数和列数
    alpha = 0.001
    #向目标移动的步长，控制速率
    maxCycles = 500
    #迭代次数
    weights = ones((n,1))
    #初始化回归系数为1(以列数即特征数量为基准)
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights)
        #梯度上升矢量化公式
        error = (labelMat - h)
        #公式计算
        weights = weights + alpha * dataMatrix.transpose()*error
    return weights.getA()
if __name__=='__main__':
    dataArr,labelMat = loadDataSet()
    weights = gradAscent(dataArr,labelMat)
    print(weights)

此处有所改动，为了方便后续的操作，所以gradAscent函数的返回值先将其从矩阵类型转为数组类型。在这里插入图片描述
由此我们便获得了其最佳回归系数。我们来解析代码，为什么要转置，我们清楚dataMatrix是mn的矩阵，而LabelMat是m1的矩阵，在参与Sigmoid函数运算时，dataMatrix与weights进行矩阵相乘，而weights是n*1的矩阵。我们知道矩阵相乘要满足前者列数要与后者的行数相等才能进行乘运算，所以要进行转置，所以在数据集列表前加入X0=1的一个原因就是方便计算。对于这个公式计算不是很理解的，大家可以看看这篇博客，梯度上升算法直观理解。

分析数据：画出决策边界
解出了回归系数，我们就确定了不同类别的数据之间的分隔线，那我们就可以将其画出来更深层次感受回归系数。用Python画函数图我们一样仍然用matplotlib。

from numpy import *
import matplotlib.pyplot as plt

def loadDataSet():
    dataMat = []
    labelMat = []
    fr = open('testSet.txt')
    #读取文本文件
    for line in fr.readlines():
        #读取全部行内容并且逐行遍历
        lineArr = line.strip().split()
        #对单行进行处理，伤处空格换行符并且切分
        dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])
        #将文本中的数据的前两个x1,x2和1.0作为x0一起放入列表中
        labelMat.append(int(lineArr[2]))
        #将标签放入标签列表
    return dataMat,labelMat

def sigmoid(inX):#sigmoid函数计算
    return 1.0/(1+exp(-inX))

def gradAscent(dataMatIn,classLabels):#梯度上升优化算法求最佳回归系数
    dataMatrix = mat(dataMatIn)
    #将训练特征数据转为numpy矩阵
    labelMat = mat(classLabels).transpose()
    #将标签列表转为numpy矩阵并将其转置
    m,n = shape(dataMatrix)
    #获得矩阵的行数和列数
    alpha = 0.001
    #向目标移动的步长，控制速率
    maxCycles = 500
    #迭代次数
    weights = ones((n,1))
    #初始化回归系数为1(以列数即特征数量为基准)
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights)
        #梯度上升矢量化公式
        error = (labelMat - h)
        #公式计算
        weights = weights + alpha * dataMatrix.transpose()*error
    return weights.getA()

def plotBestFit(weights):#绘制数据集和Logistic回归最佳拟合直线
    dataMat,labelMat = loadDataSet()
    #获取数据集和标签集
    dataArr = array(dataMat)
    #将数据集转为array数组
    n = shape(dataArr)[0]
    #获取数据集的行数，即数据的个数
    xcord1 = []
    ycord1 = []
    xcord2 = []
    ycord2 = []
    for i in range(n):#遍历数据
        if int(labelMat[i]) == 1:#若分类标签为1
            xcord1.append(dataArr[i,1])
            ycord1.append(dataArr[i,2])
            #添加正样本
        else:#类别为0
            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',alpha=.5)
    ax.scatter(xcord2,ycord2,s=30,c="green",alpha=.5)
    # 绘制散点图，前两个是数组，s是标量，c是颜色设置，marker是点的形状设置(
    #默认是圆，即o，s是正方形),alpha是透明度，在0-1之间
    x = arange(-3.0,3.0,0.1)
    #与range区别是可以为小数，指从-3到3以间隔为0.1递增
    y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x,y)
    #绘制函数图线
    plt.xlabel('X1')
    plt.ylabel('X2')
    #绘制label
    plt.show()
    #展示
if __name__='__main__':
	dataArr,labelMat = loadDataSet()
    weights = gradAscent(array(dataArr),labelMat)
	plotBestFit(weights)

在这里插入图片描述
代码中，将标签为1的数据和标签为0的数据区分画出，并且以X1的特征值作为X轴的坐标，X2的特征值作为Y轴的坐标，并且确定拟合直线的长度范围及增长趋势。对于y = (-weights[0]-weights[1]*x)/weights[2]其实也很好理解。

我们之前提到过，Sigmoid函数是以0.5作为分界线的，概率小于0.5的我们划分为0类别，大于0.5的我们划分为1类别。而我们所要画的拟合直线就是区分两者类别的直线，而当Sigmoid函数为0.5时,x的取值就是为0，也就代表着Z=0，所以我们可以写出公式0=w₀x₀+w₁x₁+w₂x₂，w是最佳系数，x是特征值，由于X₀都是1.0，所以公式可以写成0=w₀+w₁x₁+w₂x₂，上述说到X1作为x轴，X2作为y轴，所以转换之后可以写成0=w₀+w₁x+w₂y,再带入相应的回归参数，就是代码中的公式，在-3-3之间算出y的值，即X2。那我们也可以发现，W₀最后就是区分列别划分线的截距，即y=ax+b中的b，所以之前加入X₀=1.0就是匹配这个方程式，W₀*1.0不会改变其值。

虽说上面的分类效果不错，但是这个方法却需要大量的计算，每次的迭代都要进行300次的乘法，复杂程度有点高，所以我们对其进行改进。

训练算法：随机梯度上升
梯度上升算法每次更新回归系数的时候都需要遍历整个数据集，处理少量的数据尚可，但是如果有数十亿样本和成千上万的特征，该方法的复杂程度就很高。一种改进的方法就是一次仅用一个样本来更新回归系数，该方法称为随机梯度上升算法。由于可以在新样本到来的时对分类器进行增量式更新，因而随机梯度上升算法是一个在线学习算法。与”在线学习“相对应，一次处理所有数据被称作是批处理。我们写代码实现：

from numpy import *
def stocGradAscent0(dataMatrix,classLabels):#随机梯度上升算法
    m,n = shape(dataMatrix)
    #获取数据集的行数与列数
    alpha = 0.01
    #步长系数
    weights = ones(n)
    for i in range(m):
        h = sigmoid(sum(dataMatrix[i]*weights))
        error = classLabels[i]-h
        weights = weights + alpha * error * dataMatrix[i]
    return weights
if __name__=='__main__':
    dataArr,labelMat = loadDataSet()
    weights = stocGradAscent0(array(dataArr),labelMat)
    plotBestFit(weights)

在这里插入图片描述
这个分类器错分了三分之一的样本，大体上的效果还不错，而且与之前的梯度上升算法也有些区别，前者的变量h和误差error都是向量，后者全是数值，并且后者没有矩阵的转换的过程，所有变量的数据类型都是NumPy数组。我们来看看之前的梯度上升算法其最佳的三个回归系数与迭代数的变化关系，给出代码：

from numpy import *
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
def loadDataSet():
    dataMat = []
    labelMat = []
    fr = open('testSet.txt')
    #读取文本文件
    for line in fr.readlines():
        #读取全部行内容并且逐行遍历
        lineArr = line.strip().split()
        #对单行进行处理，伤处空格换行符并且切分
        dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])
        #将文本中的数据的前两个x1,x2和1.0作为x0一起放入列表中
        labelMat.append(int(lineArr[2]))
        #将标签放入标签列表
    return dataMat,labelMat
    
def sigmoid(inX):#sigmoid函数计算
    return 1.0/(1+exp(-inX))
    # return .5 * (1 + tanh(.5 * inX))

def gradAscent(dataMatIn,classLabels):#梯度上升优化算法求最佳回归系数
    dataMatrix = mat(dataMatIn)
    #将训练特征数据转为numpy矩阵
    labelMat = mat(classLabels).transpose()
    #将标签列表转为numpy矩阵并将其转置
    m,n = shape(dataMatrix)
    #获得矩阵的行数和列数
    alpha = 0.001
    #向目标移动的步长，控制速率
    maxCycles = 500
    #迭代次数
    weights = ones((n,1))
    weights_array = array([])
    #初始化回归系数为1(以列数即特征数量为基准)
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights)
        #梯度上升矢量化公式
        error = (labelMat - h)
        #公式计算
        weights = weights + alpha * dataMatrix.transpose()*error
        weights_array = append(weights_array,weights)
        #添加到数组中
    weights_array = weights_array.reshape(maxCycles,n)
    #改变维度
    return weights.getA(),weights_array
    
def plotWeights(weights_array):
    fig = plt.figure()
    font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc",size=14)
    x = arange(0,len(weights_array),1)
    ax1 = fig.add_subplot(3,1,1)
    ax1.plot(x,weights_array[:,0])
    plt.ylabel('W0')
    ax2 = fig.add_subplot(3,1,2)
    ax2.plot(x,weights_array[:,1])
    plt.ylabel('W1')
    ax3 = fig.add_subplot(3,1,3)
    ax3.plot(x,weights_array[:,2])
    plt.xlabel('迭代次数',fontproperties=font)
    plt.ylabel('W3')
    plt.show()
    
if __name__=='__main__':
    dataArr,labelMat = loadDataSet()
    weights,weights_array = gradAscent(dataArr,labelMat)
    plotWeights(weights_array)

在这里插入图片描述

图中的波动情况显示不是很清晰，收敛情况也不明显，所以我们可以改变步长系数来让结果更加直观些，alpha = 0.01，将步长系数增大点，运行结果：
在这里插入图片描述

这样我们可以大致看出它的系数与迭代次数的波动关系，算法中总共迭代了500次，它在300多次才差不多开始收敛，我们要知道，一个判断优化算法优劣的可靠方法就是看它是否收敛，也就是说参数是否达到了稳定值，是否还会不断变化。所以我们可以再来看看随机梯度上升算法的迭代关系图，由于它迭代的次数过少，我们根本无法看出变化情况，所以我们对其在外层增加200次迭代：

from numpy import *
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
def loadDataSet():
    dataMat = []
    labelMat = []
    fr = open('testSet.txt')
    #读取文本文件
    for line in fr.readlines():
        #读取全部行内容并且逐行遍历
        lineArr = line.strip().split()
        #对单行进行处理，伤处空格换行符并且切分
        dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])
        #将文本中的数据的前两个x1,x2和1.0作为x0一起放入列表中
        labelMat.append(int(lineArr[2]))
        #将标签放入标签列表
    return dataMat,labelMat
    
def sigmoid(inX):#sigmoid函数计算
    return 1.0/(1+exp(-inX))
    # return .5 * (1 + tanh(.5 * inX))

def stocGradAscent0(dataMatrix,classLabels):
    m,n = shape(dataMatrix)
    alpha = 0.01
    weights = ones(n)
    weights_array = array([])
    for j in range(200):
        for i in range(m):
            h = sigmoid(sum(dataMatrix[i] * weights))
            error = classLabels[i] - h
            weights = weights + alpha * error * dataMatrix[i]
            weights_array = append(weights_array, weights)
    weights_array = weights_array.reshape(200*m,n)
    return weights,weights_array
    
def plotWeights(weights_array):
    fig = plt.figure()
    font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc",size=14)
    x = arange(0,len(weights_array),1)
    ax1 = fig.add_subplot(3,1,1)
    ax1.plot(x,weights_array[:,0])
    plt.ylabel('W0')
    ax2 = fig.add_subplot(3,1,2)
    ax2.plot(x,weights_array[:,1])
    plt.ylabel('W1')
    ax3 = fig.add_subplot(3,1,3)
    ax3.plot(x,weights_array[:,2])
    plt.xlabel('迭代次数',fontproperties=font)
    plt.ylabel('W3')
    plt.show()
    
if __name__=='__main__':
    dataArr,labelMat = loadDataSet()
    weights,weights_array = stocGradAscent0(array(dataArr),labelMat)
    plotWeights(weights_array)

在这里插入图片描述
虽然迭代次数改为200，但它对100个样本的内部遍历，所以回归系数的更新为100*200次，系数1比0和2都更快收敛，但他们还需要更多次的迭代。我们也可以注意到，在大的波动停止后，还有一些小的波动，产生的原因就是存在一些不能正确分类的样本点（数据集并非线性可分），每次的迭代时会引发系数的剧烈改变。所以我们希望算法可以避免来回波动，从而收敛到某个值，收敛速度也可以加快，所以我们对随机梯度上升算法进行改进。

改进的随机梯度上升算法：

from numpy import *
def stocGradscent1(dataMatrix,classLabels,numIter=150):#改进的随机梯度上升算法
    #定义了迭代次数，可更改
    m,n = shape(dataMatrix)
    #获得数据集的行列数
    weights = ones(n)
    #初始化回归系数为1
    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))
            #计算函数值
            error = classLabels[randIndex] - h
            weights = weights + alpha*error*dataMatrix[randIndex]
            #公式计算最佳回归系数
            del(dataIndex[randIndex])
            #删除适用过的数据的索引
    return weights#返回系数

我们同样看其的图的效果：在这里插入图片描述
对比之前的随机梯度上升算法，改动了两个地方，一个是alpha的取值，alpha = 4/(1.0+j+i)+0.01，每次迭代的时候都会调整，可以缓解数据上的波动或者高频波动，虽然随着迭代alpha的值会不断减小，但是永远不会减小到0，因为它还有0.01这个常数项，这样做的原因是为了保证在多次迭代后新的数据仍然具有一定的影响。如果要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获得更大的回归系数。同时降低alpha的函数中，alpha每次减少1/(j+i)，其中j是迭代次数，i是样本的下标。另外一处改动在之前也提到过，随机选取样本，从列表中选取一个值，在从列表中删除该值，这样可以减少周期性波动。

同样看看其波动情况：

from numpy import *
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
def loadDataSet():
    dataMat = []
    labelMat = []
    fr = open('testSet.txt')
    #读取文本文件
    for line in fr.readlines():
        #读取全部行内容并且逐行遍历
        lineArr = line.strip().split()
        #对单行进行处理，伤处空格换行符并且切分
        dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])
        #将文本中的数据的前两个x1,x2和1.0作为x0一起放入列表中
        labelMat.append(int(lineArr[2]))
        #将标签放入标签列表
    return dataMat,labelMat
    
def sigmoid(inX):#sigmoid函数计算
    return 1.0/(1+exp(-inX))
    # return .5 * (1 + tanh(.5 * inX))

def stocGradAscent1(dataMatrix,classLabels,numIter=150):#随机梯度上升算法
    #定义了迭代次数，可更改
    m,n = shape(dataMatrix)
    #获得数据集的行列数
    weights = ones(n)
    #初始化回归系数为1
    weights_arry = array([])
    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))
            #计算函数值
            error = classLabels[randIndex] - h
            weights = weights + alpha*error*dataMatrix[randIndex]
            #公式计算最佳回归系数
            weights_arry = append(weights_arry,weights,axis=0)
            del(dataIndex[randIndex])
            #删除适用过的数据的索引
    weights_arry = weights_arry.reshape(numIter*m,n)
    return weights,weights_arry
    
def plotWeights(weights_array):
    fig = plt.figure()
    font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc",size=14)
    x = arange(0,len(weights_array),1)
    ax1 = fig.add_subplot(3,1,1)
    ax1.plot(x,weights_array[:,0])
    plt.ylabel('W0')
    ax2 = fig.add_subplot(3,1,2)
    ax2.plot(x,weights_array[:,1])
    plt.ylabel('W1')
    ax3 = fig.add_subplot(3,1,3)
    ax3.plot(x,weights_array[:,2])
    plt.xlabel('迭代次数',fontproperties=font)
    plt.ylabel('W3')
    plt.show()
    
if __name__=='__main__':
    dataArr,labelMat = loadDataSet()
    weights,weights_array = stocGradAscent1(array(dataArr),labelMat)
    plotWeights(weights_array)

在这里插入图片描述

我们能看的出来改进的随机梯度算法收敛效果更好，它迭代了150次，我们总共有100个样本，所以我们总共更新了150*100次回归系数，也就是上图中的横坐标。它在差不多2000次的时候收敛完成，也就是迭代20次的时候收敛基本完成。对比之前的两种方法，它的效率就更加突出了。

应用实例：从疝气病症预测病马的死马率
我们已经了解了Logistic回归，也实现了梯度上升算法，下面我们用Logistic回归来预测患有疝病的马的存活问题。
这里给出的数据包含368个样本和28个特征，疝病是描述马胃肠痛的术语，然而这种病不一定源自马的胃肠问题，其它问题也可能引发马疝病。该数据集中包含了医院检测马疝病的一些指标，有的指标比较主观，有的指标难以测量，例如马的疼痛的级别。

从疝气病症预测病马的死亡率流程：

1.收集数据：给定数据文件
2.准备数据：用python解析文本文件并填充缺失值
3.分析数据：可视化并观察数据
4.训练算法：使用优化算法，找到最佳的系数
5.测试算法：为了量化回归的效果，需要观察错误率。根据错误率决定是否退回到训练阶段，通过改变迭代的次数和步长等参数来得到更好的回归系数。
6.使用算法：实现一个简单的命令行程序来收集马的疝状并输出预测结果。

这里我们给出下载地址数据下载地址，点击Source code即可下载，但是需要说明的是，这些原始数据有30%数据是缺失的（当然我们下载的是处理好的数据），所以我们可以采用一些办法来处理数据集中的数据缺失的问题。

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

可选的做法：

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

所以我们先对原始数据进行预处理，其一，将所有的缺失值必须用一个实数值来替换，因为我们使用的NumPy数据类型不允许包含缺失值。这里选择0来替换所有的缺失值，回归系数的更新公式weights = weights + alpha*error*dataMatrix[randIndex]，如果dataMatrix的某特征对应值为0，那么特征的系数将不作更新，weights=weights；对于h = sigmoid(sum(dataMatrix[randIndex]*weights))，因为Sigmoid(0)=0.5，即它对结果的预测不具有任何的倾向性。其二，如果在测试数据集中发现了一条数据的类别标签已经缺失，简单的做法就是将该条数据丢弃，这是因为类别标签与特征不同，很难确定采用某个合适的值来替换。

测试算法：用Logistic回归进行分类
下面我们应用这个例子用Logistic回归在分类上做尝试。其实分类的工作很简单，就是把测试集上的每个特征向量乘以最优化方法得来的回归系数，再将该乘积结果求和，最后输入到Sigmoid函数中即可。如果对应的Sigmoid值大于0.5就预测类别为1，否则为0。
废话不多说，直接上代码：

def classifyVector(inX,weights):#分类函数
    prob = sigmoid(sum(inX*weights))
    if prob > 0.5:
        return 1.0
    else:
        return 0.0
    #计算结果大于0.5返回类别1，否则返回0

def colicTest():
    frTrain = open('horseColicTraining.txt')
    frTest = open('horseColicTest.txt')
    #打开训练集和测试集文本
    trainingSet = []
    trainingLabels = []
    for line in frTrain.readlines():
        currLine = line.strip().split('\t')
        lineArr = []
        #将训练文本划分为列表
        for i in range(21):
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[21]))
        #将数据集和类标签加入列表
    trainWeights = stocGradscent1(array(trainingSet),trainingLabels,500)
    #随机梯度上升算法求出最佳回归系数
    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

def multiTest():#计算10次平均值
    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)))

运行的结果来看：
在这里插入图片描述
从结果上看10次的迭代之后的平均的错误率为35%，事实上这个结果并不差，因为有30%的数据缺失。当然如果调整colicTest()中的迭代次数和stochGradAscent1()中的步长，平均错误率可以降到20%左右。

我们也发现运行的过程中回出现这样的报错：RuntimeWarning: overflow encountered in exp return 1.0/(1+exp(-inX))这是因为参数值inx很大时，exp(inx)可能会发生溢出，我网上搜索别人的博客论坛得到这样的解决办法：return .5 * (1 + tanh(.5 * inX))用这个代替sigmoid()函数原来内容，虽然这段代码我目前还不是很懂，可以参考这篇博客，https://www.cnpython.com/qa/84507，运行之后在这里插入图片描述
报错解决了，但是结果与原来的也有一些偏差。