diff --git a/docs/5.Logistic回归.md b/docs/5.Logistic回归.md index 5e14a9cd..700946ed 100644 --- a/docs/5.Logistic回归.md +++ b/docs/5.Logistic回归.md @@ -41,6 +41,15 @@ Sigmoid 函数的输入记为 z ，由下面公式得到: 上图展示的，梯度上升算法到达每个点后都会重新估计移动的方向。从 P0 开始，计算完该点的梯度，函数就根据梯度移动到下一点 P1。在 P1 点，梯度再次被重新计算，并沿着新的梯度方向移动到 P2 。如此循环迭代，直到满足停止条件。迭代过程中，梯度算子总是保证我们能选取到最佳的移动方向。 +介绍一下几个相关的概念： +``` +例如：y = w1x1 + w2x2 + ... + wnxn +梯度：参考上图的例子，二维图像，x方向是代表第一个系数，也就是w1，y方向代表第二个系数也就是w2，这样的向量就是梯度。 +α：上面的梯度算法的迭代公式中的阿尔法，这个代表的是移动步长。移动步长会影响最终结果的拟合程度，最好的方法就是随着迭代次数更改移动步长。 +步长通俗的理解，100米，如果我一步走10米，我需要走10步；如果一步走20米，我只需要走5步。这里的一步走多少米就是步长的意思。 +▽f(w)：代表沿着梯度变化的方向。 +``` + 上图中的梯度上升算法沿梯度方向移动了一步。可以看到，梯度算子总是指向函数值增长最快的方向。这里所说的是移动方向，而未提到移动量的大小。该量值称为步长，记作 α 。用向量来表示的话，梯度上升算法的迭代公式如下: ![梯度上升迭代公式](../images/5.Logistic/LR_9.png) @@ -268,124 +277,250 @@ def stocGradAscent0(dataMatrix, classLabels): return weights ``` -![logistic回归梯度上升法](/images/5.Logistic/梯度上升算法.png "梯度上升法") +可以看到，随机梯度上升算法与梯度上升算法在代码上很相似，但也有一些区别: 第一，后者的变量 h 和误差 error 都是向量，而前者则全是数值；第二，前者没有矩阵的转换过程，所有变量的数据类型都是 NumPy 数组。 -介绍一下几个相关的概念： -``` - 例如：y = w1x1 + w2x2 + ... + wnxn - 梯度：参考上图的例子，二维图像，x方向是代表第一个系数，也就是w1，y方向代表第二个系数也就是w2，这样的向量就是梯度。 - α：上面的梯度算法的迭代公式中的阿尔法，这个代表的是移动步长。移动步长会影响最终结果的拟合程度，最好的方法就是随着迭代次数更改移动步长。 - 步长通俗的理解，100米，如果我一步走10米，我需要走10步；如果一步走20米，我只需要走5步。这里的一步走多少米就是步长的意思。 - ▽f(w)：代表沿着梯度变化的方向。 +判断优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到了稳定值，是否还会不断地变化？下图展示了随机梯度上升算法在 200 次迭代过程中回归系数的变化情况。其中的系数2，也就是 X2 只经过了 50 次迭代就达到了稳定值，但系数 1 和 0 则需要更多次的迭代。如下图所示: + +![回归系数与迭代次数的关系图](../images/5.Logistic/LR_12.png) + +针对这个问题，我们改进了之前的随机梯度上升算法，如下: + +```python +# 随机梯度上升算法（随机化） +def stocGradAscent1(dataMatrix, classLabels, numIter=150): + m,n = shape(dataMatrix) + weights = ones(n) # 创建与列数相同的矩阵的系数矩阵，所有的元素都是1 + # 随机梯度, 循环150,观察是否收敛 + for j in range(numIter): + # [0, 1, 2 .. m-1] + dataIndex = range(m) + for i in range(m): + # i和j的不断增大，导致alpha的值不断减少，但是不为0 + alpha = 4/(1.0+j+i)+0.0001 # alpha 会随着迭代不断减小，但永远不会减小到0，因为后边还有一个常数项0.0001 + # 随机产生一个 0～len()之间的一个值 + # random.uniform(x, y) 方法将随机生成下一个实数，它在[x,y]范围内,x是这个范围内的最小值，y是这个范围内的最大值。 + randIndex = int(random.uniform(0,len(dataIndex))) + # sum(dataMatrix[i]*weights)为了求 f(x)的值， f(x)=a1*x1+b2*x2+..+nn*xn + h = sigmoid(sum(dataMatrix[randIndex]*weights)) + error = classLabels[randIndex] - h + # print weights, '__h=%s' % h, '__'*20, alpha, '__'*20, error, '__'*20, dataMatrix[randIndex] + weights = weights + alpha * error * dataMatrix[randIndex] + del(dataIndex[randIndex]) + return weights ``` -拟合程度简介 -参考： http://blog.csdn.net/willduan1/article/details/53070777 +上面的改进版随机梯度上升算法，我们修改了两处代码。 -下面是原始数据集： +第一处改进为 alpha 的值。alpha 在每次迭代的时候都会调整，这回缓解上面波动图的数据波动或者高频波动。另外，虽然 alpha 会随着迭代次数不断减少，但永远不会减小到 0，因为我们在计算公式中添加了一个常数项。 -![拟合程度图](/images/5.Logistic/拟合程度示例图.png "拟合程度示例图") +第二处修改为 randIndex 更新，这里通过随机选取样本拉来更新回归系数。这种方法将减少周期性的波动。这种方法每次随机从列表中选出一个值，然后从列表中删掉该值（再进行下一次迭代）。 -下面是拟合程度较好的： +程序运行之后能看到类似于下图的结果图。 -![拟合程度较好示例图](/images/5.Logistic/拟合程度较好示例图.png "拟合程度较好示例图") +![改进随机梯度下降结果图](../images/5.Logistic/LR_13.png) -欠拟合：模型没有很好地捕捉到数据特征，不能很好地拟合数据。 +[完整代码地址](https://github.com/apachecn/MachineLearning/blob/master/src/python/5.Logistic/logistic.py): -![欠拟合示例图](/images/5.Logistic/欠拟合示例图.png "欠拟合示例图") +### 项目案例3: 从疝气病症预测病马的死亡率 -过拟合：模型把数据学习的太彻底了，以至于把噪声数据的特征也学习到了，这样就会导致在后期测试的时候不能够很好地识别数据，即不能正确地分类，模型泛化能力太差。 +#### 项目概述 -![过拟合示例图](/images/5.Logistic/过拟合示例图.png "过拟合示例图") +使用 Logistic 回归来预测患有疝病的马的存活问题。疝病是描述马胃肠痛的术语。然而，这种病不一定源自马的胃肠问题，其他问题也可能引发马疝病。这个数据集中包含了医院检测马疝病的一些指标，有的指标比较主观，有的指标难以测量，例如马的疼痛级别。 + +#### 开发流程 -梯度上升法的伪代码如下： ``` - 某个回归系数初始化为1 - 重复R次： - 计算整个数据集的梯度 - 使用 alpha X grandient 更新回归系数的向量 - 返回回归系数 +收集数据: 给定数据文件 +准备数据: 用 Python 解析文本文件并填充缺失值 +分析数据: 可视化并观察数据 +训练算法: 使用优化算法，找到最佳的系数 +测试算法: 为了量化回归的效果，需要观察错误率。根据错误率决定是否回退到训练阶段， + 通过改变迭代的次数和步长的参数来得到更好的回归系数 +使用算法: 实现一个简单的命令行程序来手机马的症状并输出预测结果并非难事， + 这可以作为留给大家的一道习题 ``` +> 收集数据: 给定数据文件 -梯度上升算法在每次回归系数时都需要遍历整个数据集，该方法在处理100个左右的数据集时尚可，但如果有数十亿样本和成千上万的特征， -那么该方法的计算复杂度就太高了。一种改进方法是一次仅用一个样本点来更新回归系数，该方法称为随机梯度上升算法。 -由于可以在新样本到来时对分类器进行增量式更新，因而随机梯度上升算法是一个在线学习算法。与“在线学习”相对应， -一次处理所有数据被称作是“批处理”。 - 随机梯度上升算法可以写成如下的伪代码： +病马的训练数据已经给出来了，如下形式存储在文本文件中: - 所有回归系数初始化为1 - 对数据集中每个样本 - 计算该样本的梯度 - 使用 alpha X gradient 更新回归系数值 - 返回回归系数值 +``` +1.000000 1.000000 39.200000 88.000000 20.000000 0.000000 0.000000 4.000000 1.000000 3.000000 4.000000 2.000000 0.000000 0.000000 0.000000 4.000000 2.000000 50.000000 85.000000 2.000000 2.000000 0.000000 +2.000000 1.000000 38.300000 40.000000 24.000000 1.000000 1.000000 3.000000 1.000000 3.000000 3.000000 1.000000 0.000000 0.000000 0.000000 1.000000 1.000000 33.000000 6.700000 0.000000 0.000000 1.000000 +``` -## 梯度下降最优化算法 +> 准备数据: 用 Python 解析文本文件并填充缺失值 -你最经常听到的应该是梯度下降算法，它与这里的梯度上升算法是一样的，只是公式中的加法需要变成减法。因此，对应的公式可以写成： - w:=w-a▽f(w) -梯度上升算法用来求函数的最大值，而梯度下降算法用来求函数的最小值。 +处理数据中的缺失值 -## 数据中的缺失项处理 - -数据中的缺失值是个非常棘手的问题，有很多文献都致力于解决这个问题。这个问题没有标准答案，取决于实际应用中的需求。 -那么，数据缺失究竟带来了什么问题？假设有100个样本和20个特征， -这些数据都是机器手机回来的。若机器上的某个传感器损坏导致一个特征无效时该怎么办？此时是否要扔掉整个数据？这种情况下，另外19个特征怎么办？ +假设有100个样本和20个特征，这些数据都是机器收集回来的。若机器上的某个传感器损坏导致一个特征无效时该怎么办？此时是否要扔掉整个数据？这种情况下，另外19个特征怎么办？ 它们是否还可以用？答案是肯定的。因为有时候数据相当昂贵，扔掉和重新获取都是不可取的，所以必须采用一些方法来解决这个问题。 -* 下面给出了一些可选的做法： - * 使用可用特征的均值来填补缺失值； - * 使用特殊值来填补缺失值，如 -1； - * 忽略有缺失值的样本； - * 使用有相似样本的均值添补缺失值； - * 使用另外的机器学习算法预测缺失值。 -## Logistic 回归总结 +下面给出了一些可选的做法： +* 使用可用特征的均值来填补缺失值； +* 使用特殊值来填补缺失值，如 -1； +* 忽略有缺失值的样本； +* 使用有相似样本的均值添补缺失值； +* 使用另外的机器学习算法预测缺失值。 -* 逻辑回归(Logistic Regression) - * 5.1 分类问题 - * 在分类问题中，尝试预测的是结果是否属于某一个类（例如正确或错误）。 - * 分类问题的例子有： - * 判断一封电子邮件是否是垃圾邮件； - * 判断一次金融交易是否是欺诈等等。 - * 从二元的分类问题开始讨论: - 将因变量(dependant variable)可能属于的两个类分别称为负向类（negative class）和正向类（positive class），则因变量 - y属于{0，1} - 注：其中 0 表示负向类，1 表示正向类。 - * 5.2 假说表示 +现在，我们对下一节要用的数据集进行预处理，使其可以顺利地使用分类算法。在预处理需要做两件事: +* 所有的缺失值必须用一个实数值来替换，因为我们使用的 NumPy 数据类型不允许包含缺失值。我们这里选择实数 0 来替换所有缺失值，恰好能适用于 Logistic 回归。这样做的直觉在于，我们需要的是一个在更新时不会影响系数的值。回归系数的更新公式如下: - * 5.3 判定边界 - * 在逻辑回归中，我们预测： - - 当 hθ 大于等于 0.5 时，预测 y=1 - 当 hθ 小于 0.5 时，预测 y=0 - * 根据上面绘制出的 S 形函数图像，我们知道当 - z=0时 ，g(z)=0.5 - z>0时 ，g(z)>0.5 - z<0时 ，g(z)<0.5 - 又z=θ的T次方与X的积，即： - z大于等于0时，预测：y=1 - z小于0时，预测：y=0 - * 现在假设我们有一个模型：Hθ(x)=g(θ0+θ1*x1+θ2*x2) - 并且参数θ是向量[-3 1 1]。则当-3+x1+x2大于等于0，即x1+x2大于等于3时，模型将预测y=1。 - 我们可以绘制直线x1+x2=3，这条线便是我们模型的分界线，将预测为1的区域和预测为0的区域分隔开。 - * 假使我们的数据呈现这样的分布情况，怎样的模型才能适合呢？ - 因为需要用曲线才能分隔 y=0 的区域和 y=1 的区域，我们需要二次方特征： 假设参数是Hθ(x)=g(θ0+θ1*x1+θ2*x2+θ3*(x1^2)+θ4*(x2^2)+θ4*(x2^2)) - 是[-1 0 0 1 1]，则我们得到的判定边界恰好是圆点在原点且半径为 1 的圆形。可以用非常复杂的模型来适应非常复杂形状的判定边界。 - * 5.4 代价函数 - * 5.5 简化的成本函数和梯度下降 - * 5.6 高级优化 + weights = weights + alpha * error * dataMatrix[randIndex] -## Logistic 本章小结 + 如果 dataMatrix 的某个特征对应值为 0，那么该特征的系数将不做更新，即: + + weights = weights + + 另外，由于 Sigmoid(0) = 0.5 ，即它对结果的预测不具有任何倾向性，因此我们上述做法也不会对误差造成任何影响。基于上述原因，将缺失值用 0 代替既可以保留现有数据，也不需要对优化算法进行修改。此外，该数据集中的特征取值一般不为 0，因此在某种意义上说它也满足 “特殊值” 这个要求。 + +* 如果在测试数据集中发现了一条数据的类别标签已经缺失，那么我们的简单做法是将该条数据丢弃。这是因为类别标签与特征不同，很难确定采用某个合适的值来替换。采用 Logistic 回归进行分类时这种做法是合理的，而如果采用类似 kNN 的方法就可能不太可行。 + +原始的数据集经过预处理后，保存成两个文件: horseColicTest.txt 和 horseColicTraining.txt 。 + +> 分析数据: 可视化并观察数据 + +将数据使用 MatPlotlib 打印出来，观察数据是否是我们想要的格式 + +> 训练算法: 使用优化算法，找到最佳的系数 + +下面给出 原始的梯度上升算法，随机梯度上升算法，改进版随机梯度上升算法 的代码: + +```python +# 正常的处理方案 +# 两个参数：第一个参数==> dataMatIn 是一个2维NumPy数组，每列分别代表每个不同的特征，每行则代表每个训练样本。 +# 第二个参数==> classLabels 是类别标签，它是一个 1*100 的行向量。为了便于矩阵计算，需要将该行向量转换为列向量，做法是将原向量转置，再将它赋值给labelMat。 +def gradAscent(dataMatIn, classLabels): + # 转化为矩阵[[1,1,2],[1,1,2]....] + dataMatrix = mat(dataMatIn) # 转换为 NumPy 矩阵 + # 转化为矩阵[[0,1,0,1,0,1.....]]，并转制[[0],[1],[0].....] + # transpose() 行列转置函数 + # 将行向量转化为列向量 => 矩阵的转置 + labelMat = mat(classLabels).transpose() # 首先将数组转换为 NumPy 矩阵，然后再将行向量转置为列向量 + # m->数据量，样本数 n->特征数 + m,n = shape(dataMatrix) + # print m, n, '__'*10, shape(dataMatrix.transpose()), '__'*100 + # alpha代表向目标移动的步长 + alpha = 0.001 + # 迭代次数 + maxCycles = 500 + # 生成一个长度和特征数相同的矩阵，此处n为3 -> [[1],[1],[1]] + # weights 代表回归系数， 此处的 ones((n,1)) 创建一个长度和特征数相同的矩阵，其中的数全部都是 1 + weights = ones((n,1)) + for k in range(maxCycles): #heavy on matrix operations + # m*3 的矩阵 * 3*1 的单位矩阵 ＝ m*1的矩阵 + # 那么乘上单位矩阵的意义，就代表：通过公式得到的理论值 + # 参考地址： 矩阵乘法的本质是什么？ https://www.zhihu.com/question/21351965/answer/31050145 + # print 'dataMatrix====', dataMatrix + # print 'weights====', weights + # n*3 * 3*1 = n*1 + h = sigmoid(dataMatrix*weights) # 矩阵乘法 + # print 'hhhhhhh====', h + # labelMat是实际值 + error = (labelMat - h) # 向量相减 + # 0.001* (3*m)*(m*1) 表示在每一个列上的一个误差情况，最后得出 x1,x2,xn的系数的偏移量 + weights = weights + alpha * dataMatrix.transpose() * error # 矩阵乘法，最后得到回归系数 + return array(weights) + + +# 随机梯度上升 +# 梯度上升优化算法在每次更新数据集时都需要遍历整个数据集，计算复杂都较高 +# 随机梯度上升一次只用一个样本点来更新回归系数 +def stocGradAscent0(dataMatrix, classLabels): + m,n = shape(dataMatrix) + alpha = 0.01 + # n*1的矩阵 + # 函数ones创建一个全1的数组 + weights = ones(n) # 初始化长度为n的数组，元素全部为 1 + for i in range(m): + # sum(dataMatrix[i]*weights)为了求 f(x)的值， f(x)=a1*x1+b2*x2+..+nn*xn,此处求出的 h 是一个具体的数值，而不是一个矩阵 + h = sigmoid(sum(dataMatrix[i]*weights)) + # print 'dataMatrix[i]===', dataMatrix[i] + # 计算真实类别与预测类别之间的差值，然后按照该差值调整回归系数 + error = classLabels[i] - h + # 0.01*(1*1)*(1*n) + print weights, "*"*10 , dataMatrix[i], "*"*10 , error + weights = weights + alpha * error * dataMatrix[i] + return weights + + +# 随机梯度上升算法（随机化） +def stocGradAscent1(dataMatrix, classLabels, numIter=150): + m,n = shape(dataMatrix) + weights = ones(n) # 创建与列数相同的矩阵的系数矩阵，所有的元素都是1 + # 随机梯度, 循环150,观察是否收敛 + for j in range(numIter): + # [0, 1, 2 .. m-1] + dataIndex = range(m) + for i in range(m): + # i和j的不断增大，导致alpha的值不断减少，但是不为0 + alpha = 4/(1.0+j+i)+0.0001 # alpha 会随着迭代不断减小，但永远不会减小到0，因为后边还有一个常数项0.0001 + # 随机产生一个 0～len()之间的一个值 + # random.uniform(x, y) 方法将随机生成下一个实数，它在[x,y]范围内,x是这个范围内的最小值，y是这个范围内的最大值。 + randIndex = int(random.uniform(0,len(dataIndex))) + # sum(dataMatrix[i]*weights)为了求 f(x)的值， f(x)=a1*x1+b2*x2+..+nn*xn + h = sigmoid(sum(dataMatrix[randIndex]*weights)) + error = classLabels[randIndex] - h + # print weights, '__h=%s' % h, '__'*20, alpha, '__'*20, error, '__'*20, dataMatrix[randIndex] + weights = weights + alpha * error * dataMatrix[randIndex] + del(dataIndex[randIndex]) + return weights ``` - Logistic回归的目的是寻找一个非线性函数 Sigmoid 的最佳拟合参数，求解过程可以由最优化算法来完成。 - 在最优化算法中，最常用的就是梯度上升算法，而梯度上升算法又可以简化为随机梯度上升算法。 - 随机梯度上升算法与梯度上升算法的效果相当，但占用更少的计算资源。 - 此外，随机梯度上升是一个在线算法，它可以在新数据到来时就完成参数更新，而不需要重新读取整个数据集来进行批处理运算。 + +> 测试算法: 为了量化回归的效果，需要观察错误率。根据错误率决定是否回退到训练阶段，通过改变迭代的次数和步长的参数来得到更好的回归系数 + +Logistic 回归分类函数 + +```python +def classifyVector(inX, weights): + prob = sigmoid(sum(inX*weights)) + if prob > 0.5: return 1.0 + else: return 0.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) + trainLabels.append(float(currLine[21])) + trainWeights = stocGradAscent1(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(): + 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)) ``` +> 使用算法: 实现一个简单的命令行程序来手机马的症状并输出预测结果并非难事，这可以作为留给大家的一道习题 + +[完整代码地址](https://github.com/apachecn/MachineLearning/blob/master/src/python/5.Logistic/logistic.py): + * * * -* **作者：[羊三](http://www.apache.wiki/display/~xuxin) [小瑶](http://www.apache.wiki/display/~chenyao)** +* **作者：[羊三](http://cwiki.apachecn.org/display/~xuxin) [小瑶](http://cwiki.apachecn.org/display/~chenyao)** * [GitHub地址](https://github.com/apachecn/MachineLearning): * **版权声明：欢迎转载学习 => 请标注信息来源于 [ApacheCN](http://www.apachecn.org/)** diff --git a/images/5.Logistic/LR_12.png b/images/5.Logistic/LR_12.png new file mode 100644 index 00000000..82e44c2f Binary files /dev/null and b/images/5.Logistic/LR_12.png differ diff --git a/images/5.Logistic/LR_13.png b/images/5.Logistic/LR_13.png new file mode 100644 index 00000000..580774a5 Binary files /dev/null and b/images/5.Logistic/LR_13.png differ