简介
逻辑回归应该算得上是机器学习领域必须掌握的经典算法之一,并且由于其简单有效、可并行化、可解释性强等优点,至今仍然也是分类问题中最基础和最受欢迎的算法之一。尽管它的名字里面有“回归”两字,但是实际上它是用来做分类的,“逻辑”两字是其英文名字“Logistic”的音译,实际上是来自于该方法中使用的Logit函数。逻辑回归可以用于一些常见的分类问题,比如垃圾邮件过滤、网上虚假交易识别、肿瘤恶性或者良性的判断等。逻辑回归处理的问题可以简单地看成判断是“0”还是“1”,其输出是一个介于0~1之间的值,表明预测结果属于0或者1的概率。
LR的类型
LR主要可以被分成以下3种:
- 二值逻辑回归
预测的目标只有2种可能的输出,比如肿瘤的恶性和良性、是否是垃圾邮件等。 - 多项式逻辑回归
预测的目标类别有三种及其以上,并且不考虑顺序。比如,预测那种食物更受欢迎(苹果,香蕉,橘子等)。 - 有序逻辑回归
预测的目标类别有三种及其以上,但是顺序很重要。比如,预测用户对电影的评分,从1到5,分值越大说明越受用户喜欢。
决策边界
LR主要是用来处理分类问题,其输出是一个介于0~1之间的概率值,为了预测一个数据究竟属于哪个类别,我们可以设置一个阈值。基于这个阈值,我们就可以来对数据进行分类。下面是逻辑回归过程中决策平面的变化过程:
线性回归和逻辑回归的区别
线性回归的输出是输入的加权和。逻辑回归并没有直接输出输入的加权和,而是通过了一个激活函数将加权和映射到0~1之间。通过下面的图例来展示他们之间的区别:
逻辑回归的数学定义
根据上节的内容,我们知道逻辑回归就是在线性回归模型的基础上增加了激活函数,因此我们先来看下线性回归模型的定义:
上节已经介绍了Sigmoid激活函数,LR其实就是将线性回归的表达式再放进Sigmoid激活函数中,仅此而已。于是我们就可以得到LR的完整数学表达式:
逻辑回归的损失函数
损失函数是用来衡量模型输出与真实输出之间的差异,它可以指导我们朝着将差异最小化的方向去优化我们的模型。逻辑回归处理的问题的标签就只有"0"和"1"两种,故我们可以令预测结果为正样本(类别为1)的概率为,预测结果为负样本(类别为0)的概率为-,可有如下表达式:
其中是通过上一节介绍的LR的表达式计算出来的:
我们将预测为“1”和”0“的概率综合起来,可以写成如下形式:
论数学的简洁性。😜
上述表达式可以表示单个样本的标签为的概率。假设我们的训练数据包含个已经标记好的样本,设为,这些都是已经发生过的事实,我们需要根据这些样本数据来估计出逻辑回归模型的参数。根据极大似然估计的原理,我们可以写出似然函数:
似然函数连乘的形式不便于求导,我们对上式两侧取对数将其变成累加的形式,可以得到:
注意,上图中的对应我们公式中的,即模型的预测值。
逻辑回归的训练
得到了逻辑回归的目标函数之后,接下来我们需要求出的梯度,以便后续使用SGD等算法进行优化。在对求导之前,我们先做一些准备工作,即先对求导,回想一下概率的计算公式如下:
下面正式对求导,求导过程如下:
关于梯度下降算法可以参考 深入浅出--梯度下降法及其实现- 简书。
两个问题
1. 为什么不用线性回归来做分类?
假设我们有一个关于肿瘤的数据集,它包含一个特征,即肿瘤大小,标签是肿瘤是否是恶性,我们可以将其画出来,如下:
我们在的位置画一条线,与线性回归拟合的直线相交,再做与轴的垂线,交点为,如上图中黄色所示。令所有所有处于点左侧的点为负样本,令右侧的为正样本,看起来好像也可以对样本点进行很好的分类。但是考虑训练数据中可能会有一些异常值,这些值可能会影响到最终的预测结果。还是设置阈值为0.5,假如现在多加了一个异常样本点进去训练,那么情况可能会变成下图这样:
2. 为什么逻辑回归不用MSE作为损失函数?
我们知道在线性回归中,是使用的MSE(最小均方误差)来作为损失函数的。而在逻辑回归中却变成了对数损失函数,这是为什么呢?我们先写出它们各自的表达式:
- 对数损失函数
- MSE损失函数
在上述式子中,代表真实的标签,代表模型预测的标签,代表标签的数量,我们假设,即标签要么是“0”,要么是“1”。
我们来计算一下,当真实的标签与模型预测的标签不一致时,这两种损失函数的损失值分别是多少,以及对数损失究竟比MSE好在哪里。
例子
假如我们有一个样本数据,其真实标签为“1”,模型预测的标签为“0”。
使用MSE损失函数计算出来的损失值为:
使用对数损失函数计算出来的损失值为:
因为对数函数的曲线如下:
综上所述,可以看到,MSE损失函数的值与对数损失函数值相比,不值一提。因此,当真实值与模型输出值不一致的时候,对数损失函数对逻辑回归模型预测错误的惩罚力度是非常大的。
当然,在模型预测值与真实值一致的情况下,这两个损失函数计算出来的损失值都是一样的,都是0。我们可以看出,MSE对于二分类问题并不是一个很好的选择,是因为在模型分错的情况下,损失很小,惩罚力度不够。同样,在多分类情况下,即标签是通过one-hot进行编码时,MSE仍然不是一个好的选择。 在分类场景下,我们经常使用基于梯度的方法(比如拟牛顿法,梯度下降等)来最小化损失函数,从而找到参数的最优解。然后,如果损失函数是非凸的,这类方法不能保证我们找到全局最优解,相反很可能陷入局部最小值。凸函数和非凸函数如下:
上图中蓝色的点便是函数的极小值点,对于右边的非凸函数,可以找到多个极小值点,这并不是我们希望的。
注意,这里的凸函数和我们直观理解上的意义相反。
我们先来了解一下什么是凸函数。先给出维基百科上关于凸函数的定义:
凸函数是具有如下特性的一个定义在某个向量空间的凸子集(区间)上的实值函数:对其定义域上的任意两点,总有
如下图所示:
如果对于所有的,均有,那么就是凸函数。因此,如果我们能够证明我们的损失函数的二阶导数始终大于等于0,那么我们就可以证明它是一个凸函数,那就意味着一定有一个全局最小值。接下来,我们从数学的角度来证明一下,MSE损失函数是非凸函数的,而对数损失函数是凸函数。为了简化计算,我们假设样本数据只有一个特征和一个二值标签。
-
MSE损失函数
由上述前提可得,MSE损失函数可写成:
令为的一阶导数,计算可得:
当时,
当时,
综上所述,如果逻辑回归使用MSE损失函数的话,损失函数是非凸的,模型训练的时候可能陷入局部最优解,因此是不推荐的。 -
对数损失函数
根据上述前提,先写成其损失函数:
代码实践
以下代码是《机器学习实战》中的一个例子,使用逻辑回归来预测患有疝病的马的存活率问题。代码如下:
import numpy as np
import matplotlib.pyplot as plt
class LogisticRegression():
def __init__(self, filename, alpha=0.001, MaxCycle=500):
self.filename = filename
self.maxCycles = MaxCycle
self.alpha = alpha #learning rate
def load_data_set(self):
dataMat = []; labelMat = []
fr = open(self.filename, 'r')
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(self, inx):
return 1.0/(1+np.exp(-inx))
def gradientAscent(self, dataMatIn, classLabels):
dataMatrix = np.mat(dataMatIn)
labelMat = np.mat(classLabels).transpose()
m,n = np.shape(dataMatrix)
weights = np.ones((n, 1))
for k in range(self.maxCycles):
h = self.sigmoid(dataMatrix*weights)
error = labelMat - h #equals to Yi-w^TXi
#updata weights
weights = weights + self.alpha*dataMatrix.transpose()*error
return weights
def stoc_grad_ascent0(self, data_mat, class_labels):
"""
随机梯度上升,只使用一个样本点来更新回归系数
:param data_mat: 输入数据的数据特征(除去最后一列),ndarray
:param class_labels: 输入数据的类别标签(最后一列数据)
:return: 得到的最佳回归系数
:随机梯度上升的效果不怎么好,大概错分了三分之一的样本,但是迭代的次数与m一样,即100次,而梯度下降发迭代了500次,并且没有矩阵
:转置的运算,减少了计算量
"""
m, n = np.shape(data_mat)
alpha = 0.01
weights = np.ones(n)
for i in range(m):
# sum(data_mat[i]*weights)为了求 f(x)的值, f(x)=a1*x1+b2*x2+..+nn*xn,
# 此处求出的 h 是一个具体的数值,而不是一个矩阵
h = self.sigmoid(sum(data_mat[i] * weights))
error = class_labels[i] - h
# 还是和上面一样,这个先去看推导,再写程序
weights = weights + alpha * error * data_mat[i]
return weights
def stoc_grad_ascent1(self, data_mat, class_labels, num_iter=150):
"""
改进版的随机梯度上升,使用随机的一个样本来更新回归系数, 为了解决随机梯度上升算法中,由于数据集并非线性可分的问题,
在每次迭代的时候都会导致系数的剧烈变化的问题,我们期望算法能够避免来回波动,从而收敛到某个值,通过改变alpha学习率因子
:param data_mat: 输入数据的数据特征(除去最后一列),ndarray
:param class_labels: 输入数据的类别标签(最后一列数据
:param num_iter: 迭代次数
:return: 得到的最佳回归系数
"""
m, n = np.shape(data_mat)
weights = np.ones(n)
for j in range(num_iter):#默认迭代次数为150
# 这里必须要用list,不然后面的del没法使用
data_index = list(range(m))
for i in range(m):
# i和j的不断增大,导致alpha的值不断减少,但是不为0
alpha = 4 / (1.0 + j + i) + 0.01
# 随机产生一个 0~len()之间的一个值
# random.uniform(x, y) 方法将随机生成下一个实数,它在[x,y]范围内,x是这个范围内的最小值,y是这个范围内的最大值。
rand_index = int(np.random.uniform(0, len(data_index)))
h = self.sigmoid(np.sum(data_mat[data_index[rand_index]] * weights))
error = class_labels[data_index[rand_index]] - h
weights = weights + alpha * error * data_mat[data_index[rand_index]]
del (data_index[rand_index]) #每次随机选一个值,然后使用完之后删掉
return weights
def plot_best_fit(self, weights):
dataMat, labelMat = self.load_data_set()
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()
def colic_test(self):
"""
打开测试集和训练集,并对数据进行格式化处理,其实最主要的的部分,比如缺失值的补充(真的需要学会的),人家已经做了
:return:
"""
f_train = open('./5.Logistic/HorseColicTraining.txt', 'r')
f_test = open('./5.Logistic/HorseColicTest.txt', 'r')
training_set = []
training_labels = []
# 解析训练数据集中的数据特征和Labels
# trainingSet 中存储训练数据集的特征,trainingLabels 存储训练数据集的样本对应的分类标签
for line in f_train.readlines():
curr_line = line.strip().split('\t')
if len(curr_line) == 1:
continue # 这里如果就一个空的元素,则跳过本次循环
line_arr = [float(curr_line[i]) for i in range(21)] #每一行21个数据
training_set.append(line_arr)
training_labels.append(float(curr_line[21]))
# 使用 改进后的 随机梯度下降算法 求得在此数据集上的最佳回归系数 trainWeights
train_weights = self.stoc_grad_ascent1(np.array(training_set), training_labels, 500)
error_count = 0
num_test_vec = 0.0
# 读取 测试数据集 进行测试,计算分类错误的样本条数和最终的错误率
for line in f_test.readlines():
num_test_vec += 1
curr_line = line.strip().split('\t')
if len(curr_line) == 1:
continue # 这里如果就一个空的元素,则跳过本次循环
line_arr = [float(curr_line[i]) for i in range(21)]
if int(self.classify_vector(np.array(line_arr), train_weights)) != int(curr_line[21]):
error_count += 1
error_rate = error_count / num_test_vec
print('the error rate is {}'.format(error_rate))
return error_rate
def multi_test(self):
"""
调用 colicTest() 10次并求结果的平均值
:return: nothing
"""
num_tests = 10
error_sum = 0
for k in range(num_tests):
error_sum += self.colic_test()
print('after {} iteration the average error rate is {}'.format(num_tests, error_sum / num_tests))
# do the classification work
def classify_vector(self, inx, weights):
prob = self.sigmoid(np.sum(inx * weights))
if prob > 0.5:
return 1.0
return 0.0
def test(self, method=1):
data_arr, class_labels = self.load_data_set()
if method == 1:
# 注意,这里的grad_ascent返回的是一个 matrix, 所以要使用getA方法变成ndarray类型
weights = self.gradientAscent(data_arr, class_labels).getA()
elif method == 2:
weights = self.stoc_grad_ascent0(np.array(data_arr), class_labels)
else:
weights = self.stoc_grad_ascent1(np.array(data_arr), class_labels)
self.plot_best_fit(weights)
if __name__ == '__main__':
lr = LogisticRegression('./5.Logistic/TestSet.txt');
lr.test(3)
lr.multi_test()
以下是程序拟合出来的决策平面:关于代码使用的疝病数据集,可以从这里下载。https://github.com/apachecn/AiLearning/blob/master/docs/ml/5.md
上述代码是手动实现了包括梯度下降等训练方法,实际上Sklearn库为我们高度封装了这些基础算法。下面以莺尾花分类例子来展示如何在Sklearn中使用逻辑回归模型。
首先导入必要的包
# 导入必要的几个包
import matplotlib.pyplot as plt
import numpy as np
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
sklearn.datasets中的鸢尾花数据集一共包含4个特征变量,1个类别变量。共有150个样本,这里存储了其萼片和花瓣的长宽,共4个属性,鸢尾植物分三类,种类分别为山鸢尾、杂色鸢尾、维吉尼亚鸢尾。加载数据集并且进行划分:
# 载入数据集,Y的值有0,1,2三种情况,每种特征50个样本
iris = load_iris()
X = iris.data[:, :2] #获取花卉两列数据集
Y = iris.target
# 划分训练集和测试集
x_train, x_test, y_train, y_test = train_test_split(X,Y, test_size = 0.3, random_state = 0)
构造逻辑回归模型实例,并且进行训练:
#逻辑回归模型,C=1e5表示目标函数。
lr = LogisticRegression(C=1e5)
lr = lr.fit(X,Y)
对模型进行评估:
print("Logistic Regression模型训练集的准确率:%.3f" %lr.score(x_train, y_train))
print("Logistic Regression模型测试集的准确率:%.3f" %lr.score(x_test, y_test))
输出:在测试集上进行测试:
from sklearn import metrics
y_hat = lr.predict(x_test)
accuracy = metrics.accuracy_score(y_test, y_hat) #错误率,也就是np.average(y_test==y_pred)
print("Logistic Regression模型正确率:%.3f" %accuracy)
结果:对数据进行可视化:
# Plot the decision boundary. For that, we will assign a color to each
# point in the mesh [x_min, x_max]x[y_min, y_max].
x1_min, x1_max = X[:, 0].min() - .5, X[:, 0].max() + .5 # 第0列的范围
x2_min, x2_max = X[:, 1].min() - .5, X[:, 1].max() + .5 # 第1列的范围
h = .02
x1, x2 = np.meshgrid(np.arange(x1_min, x1_max, h), np.arange(x2_min, x2_max, h)) # 生成网格采样点
grid_test = np.stack((x1.flat, x2.flat), axis=1) # 测试点
grid_hat = lr.predict(grid_test) # 预测分类值
# grid_hat = lr.predict(np.c_[x1.ravel(), x2.ravel()])
grid_hat = grid_hat.reshape(x1.shape)
plt.figure(1, figsize=(6, 5))
# 预测值的显示, 输出为三个颜色区块,分布表示分类的三类区域
plt.pcolormesh(x1, x2, grid_hat,cmap=plt.cm.Paired)
# plt.scatter(X[:, 0], X[:, 1], c=Y,edgecolors='k', cmap=plt.cm.Paired)
plt.scatter(X[:50, 0], X[:50, 1], marker = '*', edgecolors='red', label='setosa')
plt.scatter(X[50:100, 0], X[50:100, 1], marker = '+', edgecolors='k', label='versicolor')
plt.scatter(X[100:150, 0], X[100:150, 1], marker = 'o', edgecolors='k', label='virginica')
plt.xlabel('Sepal length')
plt.ylabel('Sepal width')
plt.legend(loc = 2)
plt.xlim(x1.min(), x1.max())
plt.ylim(x2.min(), x2.max())
plt.title("Logistic Regression classification result", fontsize = 15)
plt.xticks(())
plt.yticks(())
plt.grid()
plt.show()
结果:参考
- 《机器学习实战》-- Peter Harrington
- 《深度学习推荐系统》-- 王喆
- https://zhuanlan.zhihu.com/p/44591359
- http://www.siyuanblog.com/?p=784
- http://www.whjtop.cn/iris/iris.html
- https://towardsdatascience.com/logistic-regression-detailed-overview-46c4da4303bc
- https://towardsdatascience.com/building-a-logistic-regression-in-python-301d27367c24
- https://towardsdatascience.com/introduction-to-logistic-regression-66248243c148
- https://towardsdatascience.com/understanding-logistic-regression-9b02c2aec102
- https://towardsdatascience.com/introduction-to-linear-regression-and-polynomial-regression-f8adc96f31cb
- https://towardsdatascience.com/why-not-mse-as-a-loss-function-for-logistic-regression-589816b5e03c