引言

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

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

"回归"一词源于最佳拟合，表示要找到最佳拟合参数集。训练分类器时的做法就是寻找最佳拟合参数，使用的是最优化算法。

Logistic回归的一般过程

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

一、基于Logistic回归和Sigmoid函数的分类

Logistic回归

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

海维塞德阶跃函数（Heaviside step function）或单位阶跃函数：能接受所有的输入然后预测出类别。问题在于：该函数在跳跃点上从0瞬间跳跃到1，这个瞬间跳跃过程有时很难处理。

另一个函数也有类似的性质，且数学上更易处理，这就是Sigmoid函数。

为了实现Logistic回归分类器，我们可以在每个特征上都乘以一个回归系数，然后把所有的结果值相加，将这个总和带入Sigmoid函数中，进而得到一个范围在0～1之间的数值。任何大于0.5的数据被分入1类，小于0.5即被归于0类。所以，Logistic回归也可以被看成是一种概率估计。

阶跃函数：在数学中，如果实数域上的某个函数可以用半开区间上的指示函数的有限次线性组合来表示。

指示函数：定义在某集合X上的函数，表示其中有哪些元素属于某一子集A。

二、基于最优化方法的最佳回归系数确定

Sigmoid函数的输入记作z，由下面公式得出：

z = w0x0+w1x1+w2x2+...+wnxn 或 z = wTx (向量写法)

其中，向量x是分类器的输入数据，向量w就是要找的最佳参数（系数），从而使得分类器尽可能地精确。

2.1 梯度上升法

梯度上升法基于的思想是：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。如果梯度记为▽，则函数f(x,y)的梯度由下式表示：

这个梯度意味着要沿x的方向移动

沿y的方向移动

其中，函数f(x,y)必须要在待计算的点上有定义并且可微。

梯度算子总是指向函数值增长最快的方向。这里所说的是移动方向，而未提到移动量的大小。该量值称为步长，记作α。用向量来表示的话，梯度算法的迭代公式如下：

该公式将一直被迭代执行，直至达到某个停止条件为止，比如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。

梯度下降算法，与这里的梯度上升算法是一样的，只是公式中的加法需要变成减法。用来求函数的最小值。

基于上面的内容，看一个Logistic回归分类器的应用例子。

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

上图简单数据集中有100个样本点，每个点包含两个数值型特征：X1和X2.在此数据集上，我们将通过使用梯度上升法找到最佳回归系数，也就是拟合出Logistic回归模型的最佳参数。

梯度上升法的伪代码如下：

每个回归系数初始化为1
重复R次：计算整个数据集的梯度使用alpha✖️gradient更新回归系数的向量返回回归系数

Logistic回归梯度上升优化算法

from numpy import *def loadDataSet():dataMat = []; labelMat = []fr = open('./datas/testSet.txt')# 读取文件，文件包含100行，3列数据：X1，X2，labelfor line in fr.readlines():lineArr = line.strip().split()dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])# 设置X0的值为1.0labelMat.append(int(lineArr[2]))return dataMat, labelMatdef sigmoid(inX):return 1.0/(1+exp(-inX))def gradAscent(dataMatIn, classLabels):# 转换为Numpy矩阵数据类型# dataMatIn为100✖️3矩阵（X0,X1,X2）dataMatrix = mat(dataMatIn)  # 2维Numpy数组转为矩阵labelMat = mat(classLabels).transpose() # 转置，行向量转为列向量m, n = shape(dataMatrix)  # m=100,n=3# 梯度上升算法所需参数alpha = 0.001  # 向目标移动的步长maxCycles = 500  # 迭代次数weights = ones((n, 1))  # 3✖️1for k in range(maxCycles):# 矩阵相乘h = sigmoid(dataMatrix*weights)  # 100✖️1 列向量error = (labelMat - h)  # 计算真实类别与预测类别的差值 weights = weights + alpha * dataMatrix.transpose() * error  # 3X100✖️100X1=3✖️1return weights

dataArr, labelMat = loadDataSet()
weights = gradAscent(dataArr, labelMat)
print(weights)

[[ 4.12414349][ 0.48007329][-0.6168482 ]]

2.3 分析数据：画出决策边界

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

def plotBestFit(weights):import matplotlib.pyplot as plt
#     weights = wei.getA()  # 转换为numpy矩阵dataMat, labelMat = loadDataSet()dataArr = array(dataMat)n = shape(dataArr)[0]  # 行数 100xcord1 = []; ycord1 = []  # 属于分类1xcord2 = []; ycord2 = []  # 属于分类0for i in range(n):if int(labelMat[i]) == 1:xcord1.append(dataArr[i,1])ycord1.append(dataArr[i,2])else:xcord2.append(dataArr[i,1])ycord2.append(dataArr[i,2])fig = plt.figure()ax = fig.add_subplot(111)# 绘制散点图# s：点大小，默认20； marker：形状，默认为oax.scatter(xcord1, ycord1, s=30, c='red', marker='s')ax.scatter(xcord2, ycord2, s=30, c='green')# 最佳拟合直线# 设置sigmoid函数为0。0是两个分类的分界处。# 设定 0=w0x0+w1x1+w2x2，然后解出X2和X1的关系式，即分隔线的方程（X0=1）# 解得 X2=（-w0-w1x1）/ w2x = arange(-3.0, 3.0, 0.1)  # start, end step 返回array对象，且可使用float型数据y = (-weights[0]-weights[1]*x)/weights[2]ax.plot(x, y)plt.xlabel('X1'); plt.ylabel('X2')
#     plt.title("The best Logistic regression fitting line obtained by the gradient rising algorithm after 500 iterations")plt.show()

plotBestFit(weights.getA())

2.4 训练算法：随机梯度上升

梯度上升算法在每次更新回归系数时都需要遍历整个数据集，该方法在处理100个左右的数据集时尚可，但如果有数十亿样本和成千上万的特征，那么该方法的计算复杂度就太高了。

一种改进方法是一次仅用一个样本点来更新回归系数，该方法称为随机梯度上升算法。由于可以在新样本到来时对分类器进行增量式更新，因而随机梯度上升算法是一个在线学习算法。与“在线学习”相对应，一次处理所有数据被称作是“批处理”。

随机梯度上升算法伪代码

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

def stocGradAscent0(dataMatrix, classLabels):m, n = shape(dataMatrix)  # 100,3alpha = 0.01weights = ones(n)  # [1. 1. 1.]for i in range(m):h = sigmoid(sum(dataMatrix[i]*weights))  # 每个样本对应相乘再求和error = classLabels[i] - hweights = weights + alpha * error * dataMatrix[i]return weights

dataArr, labelMat = loadDataSet()
weights = stocGradAscent0(array(dataArr), labelMat)
plotBestFit(weights)

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

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

改进的随机梯度上升算法

def stocGradAscent1(dataMatrix, classLabels, numIter=150):m, n = shape(dataMatrix)weights = ones(n)for j in range(numIter):dataIndex = list(range(m))# alpha每次迭代时需要调整，会缓解数据波动或者高频波动# alpha会随着迭代次数不断减小，但永远不会减小到0.# 是为了保证在多次迭代之后新数据仍然具有一定的影响。# 如果要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获得更大的回归系数。for i in range(m):# 在降低alpha的函数中，alpha每次减少1/(j+i)。这样当j<<max(i)时，alpha就不是严格下降的。# 避免参数的严格下降也常见于模拟退火等其他优化算法中。alpha = 4 / (1.0+j+i) + 0.01# 随机选取更新，然后删掉该值再进行下一次迭代randIndex = int(random.uniform(0, len(dataIndex)))  # uniform() 方法将随机生成下一个实数，它在 [x, y] 范围内h = sigmoid(sum(dataMatrix[randIndex]*weights))error = classLabels[randIndex] - hweights = weights + alpha * error * dataMatrix[randIndex]del(dataIndex[randIndex])return weights

dataArr, labelMat = loadDataSet()
weights = stocGradAscent1(array(dataArr), labelMat)
plotBestFit(weights)

dataArr, labelMat = loadDataSet()
# 默认迭代次数为150次，这里改为500次
weights = stocGradAscent1(array(dataArr), labelMat, 500)
plotBestFit(weights)

三、示例：从疝气病症预测病马的死亡率

示例：使用Logistic回归估计马疝病的死亡率

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

该数据集中包含368个样本和28个特征。来自2010年1月11日的UCI机器学习数据库（http://archive.ics.uci.edu/ml/datasets/Horse+Colic）

该数据集中包含了医院检测马疝病的一些指标，有的指标比较主观，有的指标难以测量，例如马的疼痛级别。

另外，数据集中有30%的值是缺失的。

3.1 准备数据：处理数据中的缺失值

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

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

299个训练样本，22个特征

所作处理有：

1. 选择实数0替换所有缺失值
2. 如果一条数据的类别标签已经缺失，丢弃该条数据

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

把测试集上的每个特征向量乘以最优化方法得来的回归系数，再将该乘积结果求和，最后输入到Sigmoid函数中即可。大于0.5预测类别标签为1，否则为0

Logistic回归分类函数

def classifyVector(inX, weights):prob = sigmoid(sum(inX*weights))if prob > 0.5:return 1.0else:return 0.0def colicTest():frTrain = open('./datas/horseColicTraining.txt')frTest = open('./datas/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 = stocGradAscent1(array(trainingSet), trainingLabels, 400)errorCount = 0; numTestVec = 0.0for line in frTest.readlines():numTestVec += 1.0currLine = 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 += 1errorRate = (float(errorCount) / numTestVec)print("the error rate of this test is: %f" % errorRate)return errorRatedef multiTest():numTests = 10; errorSum = 0.0for k in range(numTests):errorSum += colicTest()print("after %d iterations the average error rate is: %f" % (numTests, errorSum/float(numTests)))

multiTest()

/Users/mac/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:16: RuntimeWarning: overflow encountered in expapp.launch_new_instance()the error rate of this test is: 0.373134
the error rate of this test is: 0.373134
the error rate of this test is: 0.358209
the error rate of this test is: 0.328358
the error rate of this test is: 0.358209
the error rate of this test is: 0.358209
the error rate of this test is: 0.208955
the error rate of this test is: 0.402985
the error rate of this test is: 0.328358
the error rate of this test is: 0.328358
after 10 iterations the average error rate is: 0.341791