《机器学习实战》之逻辑回归--基于Python3--01

Posted 2021-04-06 华少的知识宝典

分类和回归是机器学习可以解决两大主要问题，从预测值的类型上看，连续变量预测的定量输出称为回归；离散变量预测的定性输出称为分类。一般来说，回归不会用在分类问题上，但逻辑回归可以处理分类问题。其通常用于研究某些因素条件下某个结果是否发生，比如医学中根据病人的一些症状来判断他是否患有某种疾病。本文将会介绍逻辑回归的原理，并且使用逻辑回归实现从疝气病症预测病马的死亡率。

逻  辑  回  归

1、逻辑回归简介     

逻辑回归（Logistic Regression）也被称为“对数几率回归”，其用以对事物发生的可能性进行估计，或者研究某些因素条件下某个结果是否发生。例如在广告预测中，使用逻辑回归预测某广告被用户点击的可能性，然后就可以将最可能被用户点击的广告摆放在其所能看到的地方，以实现广告的精准投放。其也被用于疾病的诊断，根据病人的一些症状来判断他是否患有某种疾病。

使用线性回归，可以根据输入值x来预测其输出值。而线性回归的关键也就是要去找出模型的参数 θ 的值，其可以通过最小二乘法或梯度下降法得到。

线性回归完成的是回归拟合任务，而对于分类任务，我们同样需要一条线，但不是去拟合每个数据点，而是把不同类别的样本区分开来。也就是说，在分类问题中，我们要求其输出的y值应在（0,1）之间，但是线性回归的输出值很明显会超出这个范围，可能会大于1，或者小于0。所以，线性回归的假设函数显然不适用于分类问题，需要重新找一个新的假设函数用来拟合分类问题的期望输出。因而，我们可以通过引入联系函数，将线性回归方程z变换为g(z)，并且令g(z)的值分布在（0,1）之间，且当g(z)接近0时样本的标签为类别0，即样本属于负例，当g(z)接近1时样本的标签为类别1，即样本属于正例。这样就得到了一个分类模型。

2、sigmoid函数      

对于要满足上诉要求的的函数，你可能第一反应是阶跃函数。但是，阶跃函数在跳跃点上从0瞬间跳跃到1，这个瞬间跳跃过程有时很难处理。不过，另一个函数与阶跃函数具有相似的性质，并且在数学上更容易处理，这就是sigmoid函数。

其函数图像如下图。可以看到当坐标轴的横坐标的尺度足够大时，在x=0处sigmoid函数看起来很像阶跃函数。并且由函数图像可知，sigmod函数能够满足逻辑回归的基本要求，并且当g(z)>0.5时可将其近似认为g(z)=1，而当g(z)≤0.5时，认为g(z)=0。g(z)的输出，在某种程度上表示一个分类问题在给定x的条件下等于0或者1的概率。

3、梯度上升法   

现在已经确定了分类器的函数形式，接下来的问题就是要找到参数 θ 的值，也就是最佳回归系数。而要找到他，就可以使用梯度上升法。

梯度的本意是一个向量，表示某一函数在该点处的方向导数沿着该方向取得最大值，即函数在该点处沿着该方向（此梯度的方向）变化最快，变化率最大（为该梯度的模）。 梯度上升法基于的思想为：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。

其就好像大雾天去爬山，我们只能看清周围的路，远处都看不见，那么怎么才能到达山顶呢？方法就是先找到我们周围最陡峭的地方，也就是该点变化最快的方向，然后沿着该方向走几步，接着在看看周围最陡峭的地方，然后沿着这个方向再走几步，不断重复，就可以到达山顶了。

面对这些专业的文字描述，你可能还是不太懂，那我们来做个高中简单的数学题来理解吧。这是一个简单的求极大值的例子。首先给大家一个简单的函数：

对于这个函数，极值怎么求？显然这个函数开口向下，存在极大值，它的函数图像为：

求极值，先求函数的导数：

令导数为0，可求出x=2即取得函数f(x)的极大值。极大值等于f(2)=4。

但是真实环境中的函数不会像上面这么简单，就算求出了函数的导数，也很难精确计算出函数的极值。此时我们就可以用迭代的方法来做。就像爬坡一样，一点一点逼近极值。这种寻找最佳拟合参数的方法，就是最优化算法。爬坡这个动作用数学公式表达即为：

其中，α为步长，也就是学习速率，控制更新的幅度。效果如下图所示：

比如从(0,0)开始，迭代路径就是1->2->3->4->…->n，直到求出的x为函数极大值的近似值，停止迭代。

3、代码实现     

3.1 数据处理

这里使用书中提供的数据。

import numpy as np
"""函数说明: 加载并解析数据
Parameters: 无Returns: dataMat - 数据列表。注意：这里为 list类型 labelMat - 标签列表。注意：这里为 list类型"""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])]) # 这里将 X0 的值设为了 1.0，读取到的文件的每行前两个值分别为 X1 和 X2 labelMat.append(int(lineArr[2])) # 读取到的文件的每行的第三个值为 类别标签 return dataMat, labelMat
"""函数说明: sigmoid函数
Parameters: inX - 输入数据Returns: 计算得到的sigmoid函数值"""def sigmoid(inX):    return 1.0 / (1 + np.exp(-inX))
"""函数说明: 梯度上升法
Parameters: dataMatIn - 数据集，是一个2维NumPy数组, 每列分别代表每个不同的特征, 每行则代表每个训练样本。 classLabels - 数据类别标签。它是一个 1*100 的行向量. 为了便于矩阵计算, 需要将该 行向量 转换为 列向量, 做法是将原向量转置, 再将它赋值给labelMatReturns: weights - 返回权重 数组，即最优的参数"""def gradAscent(dataMatIn, classLabels): dataMatrix = np.mat(dataMatIn) # 将 NumPy数组 转化为 矩阵 labelMat = np.mat(classLabels).transpose() # 将数据转化为 矩阵，然后再 转置 m, n = np.shape(dataMatrix) # 得到 dataMatrix 矩阵的 大小（m行，n列）。即：m个样本数, n个特征 alpha = 0.001 # 步长 maxCycles = 500 # 迭代次数 weights = np.ones((n, 1)) # 回归系数。这里生成了一个长度与 特征数 相同（n行，1列）的 矩阵，且 里面的 所有值都为 1 for k in range(maxCycles): # 每一行 特征值 乘以 一个回归系数，然后把所有的值相加，最后将各自的总和带入Sigmoid函数中，得到一个矩阵 h = sigmoid(dataMatrix*weights) # 注意：m*n矩阵 乘以 n*1 矩阵 = m*1矩阵。所以这里的 h 不是 一个数，而是 一个列向量，其元素的个数等于样本数 error = (labelMat - h) # 计算 误差。这里得到的结果还是一个矩阵 weights = weights + alpha * dataMatrix.transpose()* error # 根据误差 更新回归系数的向量 return weights # 返回回归系数 # 这里关于 for 循环里的 详细解说，可以参考这篇博客（主要在其3.4部分的公式演绎）： # https://blog.csdn.net/achuo/article/details/51160101
if __name__ == "__main__": dataMat, labelMat = loadDataSet() weights = gradAscent(dataMat, labelMat) print(weights)

3.2 梯度上升算法

"""函数说明: 梯度上升法
Parameters: dataMatIn - 数据集，是一个2维NumPy数组, 每列分别代表每个不同的特征, 每行则代表每个训练样本。 classLabels - 数据类别标签。它是一个 1*100 的行向量. 为了便于矩阵计算, 需要将该 行向量 转换为 列向量, 做法是将原向量转置, 再将它赋值给labelMatReturns: weights - 返回权重 数组，即最优的参数"""def gradAscent(dataMatIn, classLabels): dataMatrix = np.mat(dataMatIn) # 将 NumPy数组 转化为 矩阵 labelMat = np.mat(classLabels).transpose() # 将数据转化为 矩阵，然后再 转置 m, n = np.shape(dataMatrix) # 得到 dataMatrix 矩阵的 大小（m行，n列）。即：m个样本数, n个特征 alpha = 0.001 # 步长 maxCycles = 500 # 迭代次数 weights = np.ones((n, 1)) # 回归系数。这里生成了一个长度与 特征数 相同（n行，1列）的 矩阵，且 里面的 所有值都为 1 for k in range(maxCycles): # 每一行 特征值 乘以 一个回归系数，然后把所有的值相加，最后将各自的总和带入Sigmoid函数中，得到一个矩阵 h = sigmoid(dataMatrix*weights) # 注意：m*n矩阵 乘以 n*1 矩阵 = m*1矩阵。所以这里的 h 不是 一个数，而是 一个列向量，其元素的个数等于样本数 error = (labelMat - h) # 计算 误差。这里得到的结果还是一个矩阵 weights = weights + alpha * dataMatrix.transpose()* error # 根据误差 更新回归系数的向量    return weights  # 返回回归系数

 这里关于 for 循环里的 详细解说，可以参考这篇博客（主要在其3.4部分的公式演绎）：https://blog.csdn.net/achuo/article/details/51160101

到这里，可以说基本完成了一个简单的逻辑回归。我们可以运行一下，来看看使用这个梯度上升算法最后得到的回归系数。

if __name__ == "__main__": dataMat, labelMat = loadDataSet() weights = gradAscent(dataMat, labelMat) print(weights)

运行结果：

3.3 画出决策边界

"""函数说明: 画出 数据集 和 Logistic回归最佳拟合直线
Parameters: dataMat - 数据 labelMat - 样本的类别标签 weights - 回归系数Returns: 无"""def plotBestFit(dataMat, labelMat, weights): import matplotlib.pyplot as plt
 dataArr = np.array(dataMat) # 将 dataMat 转换为 array 数组 n = np.shape(dataArr)[0] # 返回数据的个数。这里取里面标签为0处的数据，即 行数 # 正样本(分类为1）和负样本（分类为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]) # 取当前遍历到的第i行数据，取其前两个数据（也就是x和y，第三个数据为类别标签，在loadDataSet函数里就已经将其抽取到labelMat里了） else: # 数据的 类别标签 不为 1，那么就将其 放到 负样本 里面 xcord2.append(dataArr[i, 1]); ycord2.append(dataArr[i, 2])
 fig = plt.figure() # 创建 绘图 的 对象 ax = fig.add_subplot(111) # 添加子图，1行1列，在第1个位置 # 绘制 数据集 ax.scatter(xcord1, ycord1, s=30, c='red', marker='s') ax.scatter(xcord2, ycord2, s=30, c='green') # 绘制 Logistic回归最佳拟合直线 x = np.arange(-3.0, 3.0, 0.1) # 为 x 取点。-3.0到3.0之间，步长为0.1 y = (-weights[0]-weights[1]*x)/weights[2] # 计算 y。w0+w1*x+w2*y=0 => y = (-w0-w1*x)/w2 # 如果这里出现报错： # x and y must have same first dimension, but have shapes (60,) and (1, 60) # 这是说：x和y的第一维度必须相同，但一个是60*1，一个是1*60。所以将y转置一下就好了。 y = y.transpose() # 将y转置 ax.plot(x, y) # 绘制出 拟合直线 plt.xlabel('X1'); plt.ylabel('X2') # 给坐标轴添加标签 plt.show()

运行一下：

if __name__ == "__main__": dataMat, labelMat = loadDataSet() weights = gradAscent(dataMat, labelMat) plotBestFit(dataArr,labelMat,weights)

结果：

可以看到，分类结果还是不错的，只分错了2到4个点。但虽然这里的数据集很小，却需要大量的计算，所以算法还有待改进。

 4、梯度上升算法改进   

4.1 随机梯度上升

梯度上升算法每次更新回归系数时都要遍历整个数据集，所以当数据集很大时，其计算复杂度就会变得很高。这里可以用一种 改进方法，那就是一次仅用一个样本点来更新回归系数，也就是随机梯度上升算法。由于可以在样本到来时对分类器进行增量式更新，因而随机梯度上升算法是一个“在线学习”算法，相应的，像上面那样一次处理所有数据的被称为“批处理”。

"""函数说明: 随机梯度上升法
Parameters: dataMatIn - 数据集，是一个2维NumPy数组, 每列分别代表每个不同的特征, 每行则代表每个训练样本。 classLabels - 数据类别标签。它是一个 1*100 的行向量。Returns: weights - 返回权重 数组，即最优的参数"""def stocGradAscent0(dataMatrix, classLabels): m, n = np.shape(dataMatrix) # 得到 dataMatrix 矩阵的 大小（m行，n列）。即：m个样本数, n个特征 alpha = 0.01 # 步长 weights = np.ones(n) # 初始化长度为n的数组, 元素全部为 1 for i in range(m): # sum(dataMatrix[i]*weights)是求 f(x)的值，f(x) = a1*x1+b2*x2+..+nn*xn, # 此处求出的 h 是一个具体的数值，而不是一个矩阵 h = sigmoid(sum(dataMatrix[i]*weights)) error = classLabels[i] - h # 计算真实类别与预测类别之间的差值，即 误差 weights = weights + alpha * error * dataMatrix[i] # 根据误差调整回归系数    return weights

由代码可以看出，随机梯度上升算法和梯度上升算法在代码上很相似，其主要区别在于：

• 梯度上升算法的变量h和误差error都是向量，而随机梯度上升都是数值；

• 随机梯度上升没有矩阵的转换过程，所有变量的数据类型都是NumPy数组。

我们来看看效果：

if __name__ == "__main__": dataMat, labelMat = loadDataSet() weights = stocGradAscent0(np.array(dataMat), labelMat) # 这里 注意 第一个参数的类型 plotBestFit(dataArr,labelMat,weights)

注：weights = stocGradAscent0(np.array(dataMat),labelMat)

这里的第一个参数为什么需要np.array(dataMat)?

在梯度上升函数gradAscend函数中，两个参数都是array_like的，也就是list, array都可以，而在随机梯度上升stocGradAscent0函数中, classLabels是array_like的, 而dataMatrix则必须是array. 

如果dataMatrix传入和gradAscend一样的list则会报错：

numpy.float64不能解释为整数.

其实错误就出现在error*datamatrix[i]这里，error是numpy.float64，而dataMatrix[i]则是一个含有三个元素的列表。

呃，好像误差有点大！！

4.2 再次改进

判断优化算法优劣的可靠方法是看它是否收敛,  也就是说参数是否达到了稳定值，是否还会不断地变化？

下图展示了随机梯度上升算法在200次迭代过程中回归系数的变化情况。可以看到，其中的系数2（也就是 X2 ）只经过了 50 次迭代就达到了稳定值,  但系数 1 和 0 则需要更多次的迭代。

产生这种现象的原因在于存在一些不能正确分类的样本点（数据集并非线性可分），在每次迭代时会引发系数的剧烈变化。我们期望的是算法能够避免来回避波动，从而收敛到某个值，另外，收敛速度也要加快。

对此，我们对随机梯度上升算法进行改进。

"""函数说明: 改进的 随机梯度上升法
Parameters:Parameters: dataMatIn - 数据集，是一个2维NumPy数组, 每列分别代表每个不同的特征, 每行则代表每个训练样本。 classLabels - 数据类别标签。它是一个 1*100 的行向量。 numIter - 迭代次数，整数Returns: weights - 返回权重 数组，即最优的参数"""def stocGradAscent1(dataMatrix, classLabels, numIter=150): m, n = np.shape(dataMatrix) weights = np.ones(n) for j in range(numIter): dataIndex = list(range(m)) # 产生一个0至m的列表 for i in range(m): alpha = 4/(1.0+j+i)+0.0001 # 每次迭代时都调整 步长。随着i和j的不断增大，alpha的值会不断减少，但是不为0，因为里面有一个常数项 randIndex = int(np.random.uniform(0, len(dataIndex))) # 随机产生一个 0～len()之间的一个 整数值，用于后面 通过随机选取样本来更新回归系数 # sum(dataMatrix[i]*weights)是求 f(x)的值，f(x) = a1*x1+b2*x2+..+nn*xn, # 此处求出的 h 是一个具体的数值，而不是一个矩阵 h = sigmoid(sum(dataMatrix[randIndex]*weights)) error = classLabels[randIndex] - h weights = weights + alpha * error * dataMatrix[randIndex] del(dataIndex[randIndex]) # 使用完选出的值后，将其从列表中删掉。因为dataIndex列表是从0到m-1的数，所以每个数刚好就等于其索引值 return weights

运行看一下：

if __name__ == "__main__": dataMat, labelMat = loadDataSet() weights = stocGradAscent1(np.array(dataMat), labelMat) # 这里 注意 第一个参数的类型 plotBestFit(dataMat,labelMat,weights)

效果不错！

这里改进了三个地方：

首先是， alpha在每次迭代的时候都会调整，这会缓解数据波动或者高频波动。另外，虽然alpha会随着迭代次数不断减小但永远不会减小到0，这是因为其中还存在一个常数项，必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。如果要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获得更大的回归系数。另一点值得注意的是，在降低alpha的函数中，alpha每次减少 1/(j+i)，其中 j 是迭代次数，i 是样本点的下标。这样当 j<<max(i) 时，alpha就不是严格下降的。避免参数的严格下降也常见于模拟退火算法等其他优化算法中。

其次是，算法通过随机选取样本来更新回归系数。这种方法将减少周期性的波动，其每次随机从列表中选出一个值，然后从列表中删掉该值(再进行下一次迭代)。

最后，改进的算法还增加了一-个迭代次数作为第3个参数。如果该参数没有给定的话，算法将默认送代150次。如果给定，那么算法将按照新的参数值进行迭代。

5、杂七杂八          

梯度下降中 迭代更新的方式： 批量梯度下降（batch gradient descent）： 也就是是梯度下降法最原始的形式，对全部的训练数据求得误差后再对θ进行更新。 优点是每步都趋向全局最优解； 缺点是对于大量数据，由于每步要计算整体数据，训练过程慢； 随机梯度下降（stochastic gradient descent）： 每一步随机选择一个样本对θ进行更新。 优点是训练速度快； 缺点是每次的前进方向不好确定，容易陷入局部最优； 微型批量梯度下降（mini-batch gradient descent）： 每步选择一小批数据进行批量梯度下降更新θ， 属于批量梯度下降和随机梯度下降的一种折中。非常适合并行处理。

6、总结            

逻辑回归算法的计算复杂度不算高，也容易实现，主要可用于解决线性分类问题。但其容易欠拟合，且分类的精度可能不高。

逻辑回归是很多分类算法的基础组件，它的好处是输出值自然地落在0到1之间，并且有概率意义。但因为它本质上是一个线性的分类器，所以处理不好特征之间相关的情况。其拟合出来的参数就代表了每一个特征对结果的影响。

使用逻辑回归实现从疝气病症预测病马的死亡率的代码放到下一篇文章吧，这里写的太长了。