第5章 Logistic回归
假设现在有一些数据点,我们用一条直线对这些点进行拟合(该线称为最佳拟合直线),这个拟合过程就称作回归。利用Logistic回归进行分类的主要思想是:根据现有数据对分类边界线建立回归公式,以此进行分类。
Logistic回归的一般流程
- 收集数据:可以使用任何方法。
- 准备数据:由于需要进行距离计算,因此要求数据类型为数值型。另外,结构化数据格式则最佳。
- 分析数据:采用任意方法对数据进行分析。
- 训练算法:大部分时间将用于训练,训练的目的是为了找到最佳的分类回归系数。
- 测试算法:一旦训练步骤完成,分类将会很快。
- 使用算法:首先,我们需要输入一些数据,并将其转换成对应的结构化数值;接着,基于训练好的回归系数就可以对这些数值进行简单的回归计算,判定它们属于哪个类别;在这之后,我们就可以在输出的类别上做一些其他分析工作。
5.1 基于 Logistic 回归和 Sigmoid 函数的分类
Logistic回归:
- 优点:计算代价不高,易于理解和实现。
- 缺点:容易欠拟合,分类精度可能不高。
- 适用数据类型:数值型和标称型数据。
我们想要的函数应该是,能接受所有的输入然后预测出类别。例如,在两个类的情况下,上述函数输出0或1。Sigmoid函数具体的计算公式如下:
σ
(
z
)
=
1
1
+
e
−
z
\sigma(z)=\frac{1}{1+e^{-z}}
σ(z)=1+e−z1
当x为0时,Sigmoid函数值为0.5。随着x的增大,对应的Sigmoid值将逼近于1;而随着x的减小,Sigmoid值将逼近于0。
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
def sigmoid(x):
return 1 / (1 + np.exp(-x))
x = np.linspace(-10, 10, 100)
plt.title("Logical Function")
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(-12, 12)
plt.ylim(0.0, 1.1)
plt.plot(x, sigmoid(x))
plt.plot([-12, 12], [0.5, 0.5])
plt.plot([0, 0], [0, 1])
plt.show()
为了实现Logistic回归分类器,我们可以在每个特征上都乘以一个回归系数,然后把所有的结果值相加,将这个总和代入Sigmoid函数中,进而得到一个范围在0~1之间的数值。任何大于0.5的数据被分入1类,小于0.5即被归入0类。所以,Logistic回归也可以被看成是一种概率估计。
5.2 基于最优化方法的最佳回归系数确定
Sigmoid函数的输入记为z,由下面公式得出:
z
=
w
0
x
0
+
w
1
x
1
+
w
2
x
2
+
⋯
+
w
n
x
n
z=w_{0}x_{0}+w_{1}x_{1}+w_{2}x_{2}+\cdots+w_{n}x_{n}
z=w0x0+w1x1+w2x2+⋯+wnxn
如果采用向量的写法,上述公式可以写成 z = w T x z=w^{T}x z=wTx,它表示将这两个数值向量对应元素相乘然后全部加起来即得到z值。其中的向量x是分类器的输入数据,向量w也就是我们要找到的最佳参数(系数),从而使得分类器尽可能地精确。为了寻找该最佳参数,需要用到最优化理论的一些知识。
5.2.1 梯度上升法
梯度上升法基于的思想是:要找到某函数的最大值,最好的方法是沿着该函数的梯度方向探寻。如果梯度记为∇,则函数f(x,y)的梯度由下式表示:
∇
f
(
x
,
y
)
=
(
d
f
(
x
,
y
)
d
x
d
f
(
x
,
y
)
d
y
)
\nabla f(x,y)= \left(\begin{matrix} \frac{df(x,y)}{dx} \\ \frac{df(x,y)}{dy} \end{matrix} \right)
∇f(x,y)=(dxdf(x,y)dydf(x,y))
这是机器学习中最易造成混淆的一个地方,但在数学上并不难,需要做的只是牢记这些符号的意义。这个梯度意味着要沿x的方向移动
d
f
(
x
,
y
)
d
x
\frac{df(x,y)}{dx}
dxdf(x,y),沿y的方向移动
d
f
(
x
,
y
)
d
y
\frac{df(x,y)}{dy}
dydf(x,y)。其中,函数f(x,y)必须要在待计算的点上有定义并且可微。
梯度上升算法到达每个点后都会重新估计移动的方向。从P0开始,计算完该点的梯度,函数就根据梯度移动到下一点P1。在P1点,梯度再次被重新计算,并沿新的梯度方向移动到P2。如此循环迭代,直到满足停止条件。迭代的过程中,梯度算子总是保证我们能选取到最佳的移动方向。
梯度上升算法沿梯度方向移动了一步。可以看到,梯度算子总是指向函数值增长最快的方向。这里所说的是移动方向,而未提到移动量的大小。该量值称为步长,记做
α
\alpha
α。用向量来表示的话,梯度算法的迭代公式如下:
w
:
=
w
+
α
∇
w
f
(
w
)
w:=w+\alpha\nabla_{w} f(w)
w:=w+α∇wf(w)
该公式将一直被迭代执行,直至达到某个停止条件为止,比如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。
梯度下降算法:
你最经常听到的应该是梯度下降算法,它与这里的梯度上升算法是一样的,只是公式中的加法需要变成减法。因此,对应的公式可以写成:
w : = w − α ∇ w f ( w ) w:=w-\alpha\nabla_{w} f(w) w:=w−α∇wf(w)
梯度上升算法用来求函数的最大值,而梯度下降算法用来求函数的最小值。
5.2.2 训练算法:使用梯度上升找到最佳参数
def load_data():
data = pd.read_csv('./data/05testSet.txt', header=None, sep='\t', names=['x1', 'x2', 'class'])
return data
def showData(myData):
fig = plt.figure()
ax = fig.add_subplot(111)
plt.title('Scatter plot of training data')
plt.xlabel('x1')
plt.ylabel('x2')
positive = myData[myData['class'].isin([1])]
negative = myData[myData['class'].isin([0])]
ax.scatter(positive['x1'], positive['x2'], s=30, c='purple', marker='+', label='Class 0')
ax.scatter(negative['x1'], negative['x2'], s=30, c='b', marker='o', label='Class 1')
plt.legend(loc='best')
plt.show()
mydata = load_data()
showData(mydata)
图中有100个样本点,每个点包含两个数值型特征:X1和X2。在此数据集上,我们将通过使用梯度上升法找到最佳回归系数,也就是拟合出Logistic回归模型的最佳参数。
梯度上升法的伪代码如下:
每个回归系数初始化为1
重复R次:
计算整个数据集的梯度
使用alpha × gradient更新回归系数的向量
返回回归系数
def loadDataSet():
dataMat = []
labelMat = []
fr = open('./data/05testSet.txt')
for line in fr.readlines():
lineArr = line.strip().split()
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
labelMat.append(int(lineArr[2]))
return dataMat, labelMat
def sigmoid(z):
return 1 / (1 + np.exp(-z))
def gradAscent(dataMat, classLabels, alpha, iteration):
"""
梯度上升方法找最佳参数
:param dataMat: 训练集
:param classLabels: 标签
:param alpha: 学习率
:param iteration: 迭代次数
:return:
"""
dataMatrix = np.matrix(dataMat)
labelMat = np.matrix(classLabels).transpose()
m, n = dataMatrix.shape
weights = np.ones((n, 1))
for k in range(iteration):
h = sigmoid(dataMatrix*weights)
error = (labelMat - h)
weights = weights + alpha * dataMatrix.T * error
return weights
data, lables = loadDataSet()
alpha = 0.001
iteration = 500
my_weights = gradAscent(data, lables, alpha, iteration)
print(my_weights)
[[ 4.12414349]
[ 0.48007329]
[-0.6168482 ]]
5.2.3 分析数据:画出决策边界
上面已经解出了一组回归系数,它确定了不同类别数据之间的分隔线。那么怎样画出该分隔线,从而使得优化的过程便于理解呢?
def show_boundary(data, weights):
w = weights.getA()
fig = plt.figure()
ax = fig.add_subplot(111)
plt.title('Decision Boundary')
plt.xlabel('x1')
plt.ylabel('x2')
positive = data[data['class'].isin([1])]
negative = data[data['class'].isin([0])]
ax.scatter(positive['x1'], positive['x2'], s=30, c='purple', marker='+', label='Class 0')
ax.scatter(negative['x1'], negative['x2'], s=30, c='b', marker='o', label='Class 1')
x = np.arange(-3.0, 3.0, 0.1)
y = (-w[0] - w[1]*x) / w[2]
ax.plot(x, y, color='g')
plt.legend(loc='best')
plt.show()
show_boundary(mydata, my_weights)
5.2.4 训练算法:随机梯度上升
梯度上升算法在每次更新回归系数时都需要遍历整个数据集,该方法在处理100个左右的数据集时尚可,但如果有数十亿样本和成千上万的特征,那么该方法的计算复杂度就太高了。一种改进方法是一次仅用一个样本点来更新回归系数,该方法称为随机梯度上升算法。由于可以在新样本到来时对分类器进行增量式更新,因而随机梯度上升算法是一个在线学习算法。与“在线学习”相对应,一次处理所有数据被称作是“批处理”。
随机梯度上升算法可以写成如下的伪代码:
所有回归系数初始化为1
对数据集中每个样本
计算该样本的梯度
使用alpha × gradient更新回归系数值
返回回归系数值
def stocgradAscent(dataMat, classLabels, alpha):
"""
随机梯度上升
:param dataMat: 训练集
:param classLabels: 标签
:param alpha: 学习率
:return:
"""
dataMat = np.array(dataMat)
m, n = dataMat.shape
weights = np.ones(n)
for i in range(m):
h = sigmoid(sum(dataMat[i] * weights))
error = classLabels[i] - h
weights = weights + alpha * error * dataMat[i]
return weights
alpha = 0.01
my_weights = stocgradAscent(data, lables, alpha)
my_weights = np.matrix(my_weights.reshape([3, 1]))
print(my_weights)
show_boundary(mydata, my_weights)
[[ 1.01702007]
[ 0.85914348]
[-0.36579921]]
可以看到,随机梯度上升算法与梯度上升算法在代码上很相似,但也有一些区别:第一,后者的变量h和误差error都是向量,而前者则全是数值;第二,前者没有矩阵的转换过程,所有变量的数据类型都是NumPy数组。
图5-6展示了随机梯度上升算法在200次迭代过程中回归系数的变化情况。其中的系数2,也就是图5-5中的X2只经过了50次迭代就达到了稳定值,但系数1和0则需要更多次的迭代。另外值得注意的是,在大的波动停止后,还有一些小的周期性波动。不难理解,产生这种现象的原因是存在一些不能正确分类的样本点(数据集并非线性可分),在每次迭代时会引发系数的剧烈改变。我们期望算法能避免来回波动,从而收敛到某个值。另外,收敛速度也需要加快。对于图5-6存在的问题,可以通过修改程序的随机梯度上升算法来解决。
def stocgradAscent02(dataMat, classLabels, alpha, iteration):
"""
改进后的随机梯度上升
:param dataMat: 训练集
:param classLabels: 标签
:param alpha: 学习率
:param iteration: 迭代次数
:return: 回归系数值
"""
dataMat = np.array(dataMat)
m, n = dataMat.shape
weights = np.ones(n)
for i in range(iteration):
dataIndex = list(range(m))
for j in range(m):
# alpha每次迭代时需要调整
alpha = 4 / (1.0 + i + j) + 0.01
# 随机选取样本来更新回归系数
randIndex = int(np.random.uniform(0, len(dataIndex)))
h = sigmoid(sum(dataMat[randIndex] * weights))
error = classLabels[randIndex] - h
weights = weights + alpha * error * dataMat[randIndex]
del(dataIndex[randIndex])
return weights
alpha = 0.01
iteration = 150
my_weights = stocgradAscent02(data, lables, alpha, iteration)
my_weights = np.matrix(my_weights.reshape([3, 1]))
print(my_weights)
show_boundary(mydata, my_weights)
[[13.77472812]
[ 0.89055145]
[-2.17190635]]
比较图5-7和图5-6可以看到两点不同。第一点是,图5-7中的系数没有像图5-6里那样出现周期性的波动,这归功于stocGradAscent02()里的样本随机选择机制;第二点是,图5-7的曲线比图5-6收敛的快很多,这是由于stocGradAscent02()可以收敛得更快。这次我们仅仅对数据集做了20次遍历,而之前的方法是500次。
5.3 示例:从疝气病症预测病马的死亡率
本节将使用Logistic回归来预测患有疝病的马的存活问题。这里的数据包含368个样本和28个特征。我并非育马专家,从一些文献中了解到,疝病是描述马胃肠痛的术语。然而,这种病不一定源自马的胃肠问题,其他问题也可能引发马疝病。该数据集中包含了医院检测马疝病的一些指标,有的指标比较主观,有的指标难以测量,例如马的疼痛级别。
另外需要说明的是,除了部分指标主观和难以测量外,该数据还存在一个问题,数据集中有30%的值是缺失的。下面将首先介绍如何处理数据集中的数据缺失问题,然后再利用Logistic回归和随机梯度上升算法来预测病马的生死。
5.3.1 准备数据:处理数据中的缺失值
数据中的缺失值是个非常棘手的问题,有很多文献都致力于解决这个问题。那么,数据缺失究竟带来了什么问题?假设有100个样本和20个特征,这些数据都是机器收集回来的。若机器上的某个传感器损坏导致一个特征无效时该怎么办?此时是否要扔掉整个数据?这种情况下,另外19个特征怎么办?它们是否还可用?答案是肯定的。因为有时候数据相当昂贵,扔掉和重新获取都是不可取的,所以必须采用一些方法来解决这个问题。
下面给出了一些可选的做法:
- 使用可用特征的均值来填补缺失值;
- 使用特殊值来填补缺失值,如-1;
- 忽略有缺失值的样本;
- 使用相似样本的均值添补缺失值;
- 使用另外的机器学习算法预测缺失值。
现在,我们对下一节要用的数据集进行预处理,使其可以顺利地使用分类算法。在预处理阶段需要做两件事:第一,所有的缺失值必须用一个实数值来替换,因为我们使用的NumPy数据类型不允许包含缺失值。这里选择实数0来替换所有缺失值,恰好能适用于Logistic回归。这样做的直觉在于,我们需要的是一个在更新时不会影响系数的值。回归系数的更新公式如下:
w
e
i
g
h
t
s
=
w
e
i
g
h
t
s
+
a
l
p
h
a
∗
e
r
r
o
r
∗
d
a
t
a
M
a
t
r
i
x
[
r
a
n
d
I
n
d
e
x
]
weights=weights+alpha*error*dataMatrix[randIndex]
weights=weights+alpha∗error∗dataMatrix[randIndex]
如果dataMatrix的某特征对应值为0,那么该特征的系数将不做更新,即:
w
e
i
g
h
t
s
=
w
e
i
g
h
t
s
weights=weights
weights=weights
另外,由于sigmoid(0)=0.5,即它对结果的预测不具有任何倾向性,因此上述做法也不会对误差项造成任何影响。基于上述原因,将缺失值用0代替既可以保留现有数据,也不需要对优化算法进行修改。此外,该数据集中的特征取值一般不为0,因此在某种意义上说它也满足“特殊值”这个要求。
预处理中做的第二件事是,如果在测试数据集中发现了一条数据的类别标签已经缺失,那么我们的简单做法是将该条数据丢弃。这是因为类别标签与特征不同,很难确定采用某个合适的值来替换。采用Logistic回归进行分类时这种做法是合理的,而如果采用类似kNN的方法就可能不太可行。原始的数据集经过预处理之后保存成两个文件:horseColicTest.txt和horseColicTraining.txt。
5.3.2 测试算法:用 Logistic 回归进行分类
使用Logistic回归方法进行分类并不需要做很多工作,所需做的只是把测试集上每个特征向量乘以最优化方法得来的回归系数,再将该乘积结果求和,最后输入到Sigmoid函数中即可。如果对应的Sigmoid值大于0.5就预测类别标签为1,否则为0。
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
# 使用pandas加载数据
def load_data(file_name):
data = pd.read_csv(file_name, header=None, sep='\t')
return data
def sigmoid(z):
return 1 / (1 + np.exp(-z))
def stocgradAscent02(dataMat, classLabels, alpha, iteration):
"""
改进后的随机梯度上升
:param dataMat: 训练集
:param classLabels: 标签
:param alpha: 学习率
:param iteration: 迭代次数
:return: 回归系数值
"""
dataMat = np.array(dataMat)
m, n = dataMat.shape
weights = np.ones(n)
for i in range(iteration):
dataIndex = list(range(m))
for j in range(m):
# alpha每次迭代时需要调整
alpha = 4 / (1.0 + i + j) + 0.01
# 随机选取样本来更新回归系数
randIndex = int(np.random.uniform(0, len(dataIndex)))
h = sigmoid(sum(dataMat[randIndex] * weights))
error = classLabels[randIndex] - h
weights = weights + alpha * error * dataMat[randIndex]
del(dataIndex[randIndex])
return weights
def classifyVector(inX, weights):
prob = sigmoid(sum(inX * weights))
if prob >= 0.5:
return 1.0
else:
return 0.0
def colicTest(weights):
errorCount = 0
numTestVect = 0.0
frTest = open('./data/05horseColicTest.txt')
for line in frTest.readlines():
numTestVect += 1.0
currLine = line.strip().split('\t')
lineArr = []
for i in range(21):
lineArr.append(float(currLine[i]))
if int(classifyVector(np.array(lineArr), weights)) != int(currLine[21]):
errorCount += 1
errorRate = (float(errorCount) / numTestVect)
print('the error rate of this test is: %f' % errorRate)
return errorRate
def multiTest(weights):
numTest = 10
errorSum = 0.0
for k in range(numTest):
errorSum += colicTest(weights)
print('after %d iterations the average error rate is: %f' % (numTest, errorSum/float(numTest)))
if __name__ == '__main__':
file_name = './data/05horseColicTraining.txt'
trainData = load_data(file_name)
m, n = trainData.shape
X = trainData.iloc[:, : n - 1]
labels = trainData.iloc[:, -1:]
alpha = 0.01
iters = 500
X = np.array(X)
labels = np.array(labels)
weights = stocgradAscent02(X, labels, alpha, iters)
print(weights)
multiTest(weights)
[ 23.92422463 -3.06121102 4.34721833 -2.32422927 1.60858913
-6.25905843 6.09438072 -12.47080956 -17.9416587 -16.17191663
37.75196802 -44.55228352 50.35555788 21.60621382 -25.27671323
5.56857495 -2.9613017 -0.7922193 0.11365851 -12.42177038
-5.30166081]
the error rate of this test is: 0.328358
the error rate of this test is: 0.328358
the error rate of this test is: 0.328358
the error rate of this test is: 0.328358
the error rate of this test is: 0.328358
the error rate of this test is: 0.328358
the error rate of this test is: 0.328358
the error rate of this test is: 0.328358
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.328358
从上面的结果可以看到,10次迭代之后的平均错误率为35%。事实上,这个结果并不差,因为有30%的数据缺失。当然,如果调整colicTest()中的迭代次数和stochGradAscent02()中的步长,平均错误率可以降到20%左右。
5.4 本章小结
Logistic回归的目的是寻找一个非线性函数Sigmoid的最佳拟合参数,求解过程可以由最优化算法来完成。在最优化算法中,最常用的就是梯度上升算法,而梯度上升算法又可以简化为随机梯度上升算法。
随机梯度上升算法与梯度上升算法的效果相当,但占用更少的计算资源。此外,随机梯度上升是一个在线算法,它可以在新数据到来时就完成参数更新,而不需要重新读取整个数据集来进行批处理运算。