第 5 章 Logistic 回归
[TOC]
本章内容
- Sigmoid 函数和 Logistic 回归分类器
- 最优化理论初步
- 梯度下降最优化算法
- 数据中的缺失项处理
假设现在有一些数据点，我们用一条直线对这些点进行拟合（该线称为最佳拟合直线），这个拟合过程就称为回归。利用 Logistic 回归进行分类的主要思想是：根据现有数据对分类边界建立回归公式，以此进行分类。
Logistic 回归的一般过程：
- 收集数据：采用任意方法收集数据。
- 准备数据：由于需要进行距离计算，因此要求数据类型为数值型。另外，结构化数据格式则最佳。
- 分析数据：采用任意方法对数据进行分析。
- 训练算法：大部分时间将用于训练，训练的目的是为了找到最佳的分类回归系数。
- 测试算法：一旦训练步骤完成，分类将会很快。
- 使用算法：
- 首先，我们需要输入一些数据，并将其转换成对应的结构化数值；
- 接着，基于训练好的回归系数就可以对这些数值进行简单的额回归计算，判定它们属于哪个类别；
- 在这之后，我们就可以在输出的类别上做一些其他分析工作。
1. 基于 Logistic 回归和 SIgmoid 函数的分类
Logistic 回归：
- 优点：计算代价不高，易于理解和实现
- 缺点：容易欠拟合，分类精度可能不高
- 适用数据类型：数值型和标称型数据。
我们想要的函数应该是，能接受所有的输入然后预测出类别。例如，在两个雷的情况下，上述函数输出 0 或 1，该函数称为 海维塞德阶跃函数（Heaviside step function），或者直接称为 ** 单位阶跃函数**。
海维塞德阶跃函数的问题在于：该函数在跳跃点上从 0 瞬间跳跃到 1，这个瞬间跳跃过程有时很难处理。考虑另一个函数： Sigmoid 函数：
图5-1 给出 Sigmoid 函数在不同坐标尺度下的两条曲线图。当 x 为 0 时， Sigmoid 函数值为 0.5.随着 x 的增大，对应的 Sigmoid 值将逼近于 1；而随着 x 的减小，Sigmoid 值将逼近于 0。
问题：最佳回归系数是多少？如何确定它们的大小？
2. 基于最优化方法的最佳回归系数确定
Sigmoid 函数的输入记为 z：
写成向量形式：
表示将这两数值向量对应元素相乘然后全部加起来即得到 z 值。其中的向量 X 是分类器的输入数据，向量 W 也就是我们要找到的最佳参数（系数）
2.1 梯度上升法
梯度上升法基于的思想是：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。如果 梯度记为 $$，则函数 f(x,y) 的梯度：
这个梯度意味着要沿 x 的方向移动
，沿 y 方向移动
。其中，函数 f(x,y) 必须在待计算的点互换有定义并且可微。
一个具体的函数例子见图5-2
图5-2 中的梯度上升算法沿梯度方向移动了一步。可以看到，梯度算法总是指向函数值增长最快的方向这里所说的是移动方向，而未提到移动量的大小。该量值称为步长，记作 $\alpha$。用向量来表示的话，梯度算法的迭代公式：
该公式将一直被迭代执行，直至达到某个停止条件为止，比如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。
梯度下降算法：与这里的梯度上升算法是一样的，只是公式中的加号变成减号：
梯度上升算法是用来求函数的最大值，梯度下降算法用来求函数的最小值
基于上面的内容，我们来看一个 Logistic 回归分类器的应用例子，从图5-3 可以看出我们采用的数据集：
2.2 训练算法：使用梯度上升找到最佳参数
图5-3 中有100个样本点，每个点包含两个数值型特征：X1 和 X2。在此数据集上，我们将通过使用梯度上升法找到最佳回归系数，也就是拟合出 Logistic 回归模型的最佳参数。
梯度上升法的伪代码：
每个回归系数初始化为1
重复 R 次：
计算整个数据集的梯度
使用 alpha * gradient 更新回归系数的向量
返回回归系数
添加 loadDataSet() 函数：
def loadDataSet():
"""
读取数据
:return:
"""
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() 函数：
def sigmoid(inX):
"""
sigmoid 函数
:param inX:
:return:
"""
return 1.0/(1+np.exp(-inX))
添加 gradAscent() 函数：
def gradAscent(dataMatIn, classLabels):
"""
梯度上升算法
:param dataMatIn: 2维数组，每列分别代表每个不同的特征，每行代表每个训练样本
:param classLabels: 类别标签
:return:
"""
dataMatrix = np.mat(dataMatIn)
labelMat = np.mat(classLabels).transpose() # 行向量转置成列向量
m,n = np.shape(dataMatrix)
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
2.3 分析数据：画出决策边界
上面已经解出一组回归系数，它确定了不同类别数据之间的分隔线。那么怎样画出该分隔线，从而使得优化的过程便于理解呢？
添加 plotBestFit() 函数：
def plotBestFit(wei):
"""
画出数据集和 Logistic 回归最佳拟合直线的函数
:param wei:
:return:
"""
import matplotlib.pyplot as plt
weights = wei.getA()
# weights = wei
dataMat, labelMat = loadDataSet()
dataArr = np.array(dataMat)
n = np.shape(dataArr)[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])
else:
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')
ax.scatter(xcord2, ycord2, s=30, c='green')
x = np.arange(-3.0, 3.0, 0.1)
y = (-weights[0] - weights[1] * x)/weights[2]
ax.plot(x,y)
plt.xlabel('X1')
plt.ylabel('X2')
plt.show()
测试：
dataArr,labelMat = loadDataSet()
weights=gradAscent(dataArr,labelMat)
print(weights) #[[ 4.12414349] [ 0.48007329] [-0.6168482 ]]
plotBestFit(weights)
这个分类结果相当不错，从图上看只错分了两到四个点。但是尽管例子简单且数据量很小，却需要大量的计算（300次乘法）。下一节对算法进行改进。
2.4 训练算法：随机梯度上升
随机梯度上升：一次仅用一个样本点来更新回归系数。由于可以在新样本到来时对分类器进行增量式更新，因此也是一个在线学习算法。
随机梯度上升的伪代码：
所有回归系数初始化为1
对数据集中每个样本
计算该样本的梯度
使用 alpha*gradient 更新回归系数值
返回回归系数值
以为是随机梯度上升算法的实现代码：
def stocGradAscent0(dataMatrix, classLabels):
"""
随机梯度上升算法
:param dataMatrix:
:param classLabels:
:return:
"""
m,n = np.shape(dataMatrix)
alpha = 0.01
weights = np.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
可以看到，随机梯度上升算法与梯度上升算法在代码上很相似，有的区别是：第一，后者的变量 h 和 误差 error 都是向量，而前者则是数值；第二，前者没有矩阵的转换过程，所有变量的数据类型都是 numpy 数组。
测试代码：
dataArr,labelMat = loadDataSet()
weights2 = stocGradAscent0(np.array(dataArr),labelMat)
print('weights2:',weights2) # [ 1.01702007 0.85914348 -0.36579921]
plotBestFit(weights2)
一个判断优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到了稳定值，是否还会不断地变化？。对此，我们在 stocGradAscent0() 函数上做了些修改，使其在整个数据集上运行 200 次。最终绘制的三个回归系数的变化情况如图5-6.
图5-6 展示了随机梯度上升算法在 200 次迭代过程中回归系数的变化情况。其中的系数2，也就是图5-5 中的 X2 只经过了 50 次迭代就达到稳定值，但系数 1 和 0 则需要更多次的迭代。值得注意的是，在大的波动停止后，还以一些小的周期性波动。产生这种现象的原因是存在一些不能正确分类的样本点（数据集并非线性可分），在每次迭代时会引发系数的剧烈改变
改进的随机算法 stocGradAscent1：
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
"""
改进的随机梯度上升算法
:param dataMatrix:
:param classLabels:
:return:
"""
m,n = np.shape(dataMatrix)
weights = np.ones(n)
for j in range(numIter):
dataIndex = range(m)
for i in range(m):
alpha = 4/(1.0+j+i)+0.01 # alpha 每次迭代时需要调整
randIndex = int(np.random.uniform(0,len(dataIndex))) # 随机选取更新
h = sigmoid(sum(dataMatrix[i]*weights))
error = classLabels[i] - h
weights = weights + alpha * error * dataMatrix[randIndex]
del(dataIndex[randIndex])
return weights
alpha 在每次迭代的时候都会调整，这回缓解图5-6 上的数据波动或者高频波动。另外虽然 alpha 会随着迭代次数不断减小，但永远不会减小到 0 ，这是因为调整式子中还有一个常数项。必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。如果要处理的问题是动态变化的，那么可以适当加大常数项，来确保新的值获得更大的回归系数。还有，在降低 alpha 的函数中，alpha 每次减少 1/(j+1)，其中 j 是迭代次数，i 是样本点的下标。这样当 j<<max(i) 时，alpha 就不是严格下降的，避免参数的严格下降也常见于模拟退火算法等其他优化算法中。
比较图5-7 和图5-6 可以看出两点不同： 1.图5-7中的系数没有像图5-6那样出现周期性的波动，这归功于 stocGradAscent1 里的样本随机选择机制；2.图5-7的水平轴比图5-6 短了很多，这是由于 stocGradAscent1 可以收敛得更快。这次我们仅仅对数据集做了20次遍历，而之前是500次。
程序运行后得到图5-8的结果图，该分隔线达到了与 gradAscent() 差不多的效果，但是所使用的计算量更少。
3. 示例：从疝气病症预测病马的死亡率
使用 Logistic 回归来预测患有疝气的马的存活问题。
示例：使用 Logistic 回归估计马疝病的死亡率
- 收集数据：给定数据文件。
- 准备数据：用 python 解析文本文件并填充缺失值。
- 分析数据：可视化并观察数据。
- 训练算法：使用优化算法，找到最佳的系数。
- 测试算法：为了量化回归的效果，需要观察错误率。根据错误率决定是否回退到训练阶段，通过改变迭代的次数和步长等参数得到更好的回归系数。
- 使用算法：实现一个简单的命令行程序来收集马的症状并输出预测结果。
注意：除了部分指标主观和难以测量外，该数据还存在一个问题，数据集中有 30% 的值是缺失的。
3.1 准备数据：处理数据中的缺失值
数据缺失究竟带来了什么问题？假设有100个样本和20 个特征，这些数据都是机器收集回来的。若机器的某个传感器损坏导致一个特征无效时该怎么办？此时是否要扔掉整个数据？这种情况下，另外19个特征怎么办？它们是否还柯永松？答案是肯定的、因为有时候数据相当昂贵，扔掉和重新获取都是不可取的，所以必须采取一些方法来解决这个问题：
- 使用可用特征的均值来填补缺失值。
- 使用特殊值来填补缺失值，如 -1。
- 忽略缺失值的样本。
- 使用相似样本的均值填补缺失值。
- 使用另外的机器学习算法预测缺失值。
在预处理阶段需要做两件事：
- 所有的缺失值必须用一个实数值来替换，因为我们使用的 Numpy 数据类型不允许包含缺失值。这里选择实数 0 来替换所有缺失值，恰好能适用于 Logistic 回归。
- 如果在测试数据集中发现了一条数据的类别标签已经缺失，那么我们的简单做法是将该条数据丢弃。这是因为类别标签与特征不同，很难确定采用某个合适的值来替换。
3.2 测试算法：用 Logistic 回归进行分类
使用 Logistic 回归方法进行分类并不需要做很多工作，所需做的只是把测试集上每个特征向量乘以最优化方法得来的回归系数，再将该乘积结果求和，最后输入到 Sigmoid 函数中即可。如果对应的 Sigmoid 值大于 0.5 就预测类别标签为 1，否则为 0.
添加 classifyVector 函数：
def classifyVector(inX, weights):
"""
根据 Sigmoid 值判断分类
:param inX:
:param weights:
:return:
"""
prob = sigmoid(sum(inX*weights))
if prob > 0.5:return 1.0
else: return 0.0
添加 colicTest() 函数：
def colicTest():
"""
打开训练集和测试集，并对数据进行格式化处理，以及训练、测试
:return:
"""
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 = stocGradAscent1(np.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(np.array(lineArr), trainWeights)) != int(currLine[21]):
errorCount += 1
errorRate = (float(errorCount)/numTestVec)
print('the error rate of this test is : %f' % errorRate)
return errorRate
添加 multiTest() 函数：
def multiTest():
"""
重复测试，并求取平均值
:return:
"""
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)))
测试：
multiTest()
如果调整 colicTest() 中的迭代次数和 stocGradAscent1() 中的步长，平均错误率可以降到 20% 左右。
4. 本章小结
Logistic 回归的目的是寻找一个非线性函数 Sigmoid 的最佳拟合参数，求解过程可以由最优化算法来完成。在最优化算法中，最常用的就是梯度上升算法，而梯度上升算法又可以简化为随机梯度上升算法。
