逻辑回归算法的原理 逻辑回归算法原理实例

逻辑回归算法的名字里虽然带有“回归”二字，但实际上逻辑回归算法是用来解决分类问题的算法。线性回归和逻辑回归相当于一对“孪生兄弟”，本文将从二分类入手，介绍逻辑回归算法的预测函数、损失函数（成本函数）和梯度下降算法公式，然后由二分类延伸到多分类的问题，接下来介绍正则化，即通过数学的手段来解决模型过拟合问题，最后用一个乳腺癌检测的实例及其模型性能优化来结束全文。各位朋友在看这篇博客的时候，根据提示不断联想线性回归与逻辑回归的区别与联系。

1 逻辑回归算法的原理

假设有一场球赛，我们有两支球队的所有出场球员信息、历史交锋成绩、比赛时间、主客场、裁判和天气等信息，根据这些信息预测球队的输赢。假设比赛结果记为${y}$，赢球标记为1，输球标记为0，这就是一个典型的二元分类问题，可以用逻辑回归算法来解决。
从这个例子里可以看出，逻辑回归算法的输出${y \in \{0,1}\}$是个离散值，这是与线性回归算法的最大区别。

1.1 预测函数

需要找出一个预测函数模型，使其值输出在${[0,1]}$之间。然后选择一个基准值，如${0.5}$，如果算出来的预测值大于${0.5}$，就认为其预测值为1，反之则其预测值为0。我们选择${g(z)= \frac{1}{1+e^{-z}}}$来作为预测函数。函数${g(z)}$称为${Sigmoid}$函数，也称为${Logistic}$函数。图像如下：

当${z=0}$时，${g(z)=0.5}$。

当${z>0}$时，${g(z)>0.5}$，当${z}$越来越大时，${g(z)}$无限接近于${1}$。

当${z<0}$时，${g(z)<0.5}$，当${z}$越来越小时，${g(z)}$无限接近于${0}$。

这正是我们想要的针对二元分类算法的预测函数。

问题来了，怎样把输入特征和预测函数结合起来呢？

结合线性回归函数的预测函数${h_{\theta}(x)={\theta}^Tx}$，假设令${z(x)={\theta}^Tx}$，则逻辑回归算法的预测函数如下：

${h_{\theta}(x)=g(z)=g({\theta}^Tx)=\frac{1}{1+e^{-{\theta}^Tx}}}$。下面解读预测函数。

${h_{\theta}(x)}$表示在输入值为${x}$，参数为${\theta}$的前提下${y=1}$的概率。用概率论的公式可以写成：${h_{\theta}(x)=P(y=1|x,\theta)}$，即在输入${x}$及参数${\theta}$条件下${y=1}$的概率。由条件概率公式可以推导出

${P(y=1|x,\theta)+P(y=0|x,\theta)=1}$

对二分类来说，这是一个非黑即白的世界。

1.2 判定边界

逻辑回归算法的预测函数由以下两个公式给出：
${h_{\theta}(x)=g({\theta}^Tx)}$ 和 ${g(z)= \frac{1}{1+e^{-z}}}$
假定${y=1}$的判定条件是${h_{\theta}(x)\geq0.5}$，${y=0}$的判定条件是${h_{\theta}(x)\leq 0.5}$，所以${\theta^Tx=0}$就是我们的判定边界。
假定有两个变量${x_1，x_2}$，其逻辑回归预测函数是${h_{\theta}(x)=g({\theta}_0+{\theta}_1{x_1}+{\theta}_2{x_2})}$。假设给定参数${\theta=[-3\ 1 \ 1]^T}$，那么得到判定边界为${-3+x_1+x_2=0}$。
如果预测函数是多项式${h_{\theta}(x)=g({\theta}_0+{\theta}_1{x_1}+{\theta}_2{x_2}+{\theta}_3{x_1}^2+{\theta}_4{x_2}^2)}$，给定${\theta=[-1\ 0 \ 0 \ 1 \ 1]^T}$，则得到判定边界函数为${{x_1}^2+{x_2}^2=1}$。这是二阶多项式的情况，更一般的多阶多项式可以表达出更复杂的判定边界。

1.3 损失函数（成本函数）

我们不能使用线性回归模型的损失函数来推导逻辑回归的损失函数，因为那样的损失函数太复杂，最终很可能会导致无法通过迭代找到损失函数值最小的点。为了容易地求出损失函数的最小值，我们分成${y=1}$和${y=0}$两种情况分别考虑其预测函数值与真实值的误差。我们先考虑最简单的情况，即计算某一个样本${(x,y)}$的预测值与真实值的误差，损失函数如下：

${Cost(h_{\theta}(x), y)={-log(h_{\theta}(x)), \ if \ y=1 \brace -log(1-h_{\theta}(x)), \ if \ y=0}}$，其中${h_{\theta}(x)}$表示预测为${1}$的概率。

回顾成本/损失的定义：预测值与真实值的差异。差异越大，损失越大，模型受到的“惩罚”也越严重。当${y=1}$时，随着${h_{\theta}(x)}$的值（预测为1的概率）越来越大，预测值越来越接近真实值，其成本/损失越来越小。当${y=0}$时，随着${h_{\theta}(x)}$的值（预测为1的概率）越来越大，预测值越来越偏离真实值，其成本/损失越来越大。

符合上述规律的函数模型有很多，为什么我们要选择自然对数函数作为损失函数呢？
逻辑回归模型的预测函数是${sigmoid}$函数，而${sigmoid}$函数里有${e}$的${n}$次方运算，自然对数刚好是其逆运算，比如${log(e^n)=n}$，最终会推导出形式优美的逻辑回归模型参数的迭代函数，而不需要去涉及对数运算和指数函数运算。这就是我们选择自然对数函数作为损失函数的原因。

下面将损失函数统一写法。分开表述的成本计算公式始终不方便，合并成一个公式岂不完美。然后就有了下面的公式：
${Cost(h_{\theta}(x), y)=-y*{log(h_{\theta}(x))-(1-y)*log(1-h_{\theta}(x))}}$
由于${y \in \{0,1}\}$是离散值，当${y=1}$时，${1-y=0}$，上式的后半部分为${0}$；当${y=0}$时，上式的前半部分为${0}$。因此上式与分开表达的成本计算公式是等价的。

介绍到这里，损失函数就要隆重登场了。根据一个样本的成本计算公式，很容易写出所有样本的成本/损失平均值，即损失函数：
${J(\theta)=- \frac{1}{m} \lbrack \sum_{i=1}^{m} \ y^{(i)} * log(h_\theta(x^{(i)}))+(1-y^{(i)}) * log(1-h_\theta(x^{(i)}))\rbrack}$

1.4 梯度下降算法

和线性回归类似，我们使用梯度下降算法来求解逻辑回归模型参数。根据梯度下降算法的定义，可以得出：
${\theta_j = \theta_j - \alpha \frac{\partial J(\theta)}{\partial \theta_j}}$，这里的关键同样是求解损失函数的偏导数。最终推导出来的梯度下降算法公式为：
${\theta_j = \theta_j - \alpha \frac{1}{m} \sum_{i=1}^{m} (h_\theta(x^{(i)})-y^{(i)})x_j^{(i)}}$
这个公式的形式和线性回归算法的参数迭代公式是一样。当然，由于这里${h_{\theta}(x)=\frac{1}{1+e^{-{\theta}^Tx}}}$，而线性回归算法里${h_{\theta}(x)={\theta}^Tx}$。所以，两者的形式是一样，但是数值计算方法则完全不同。
至此，我们把逻辑回归算法的相关原理解释清楚了。

2 多元分类

逻辑回归模型可以解决二分类问题，即${y=\{0,1\}}$，能否用来解决多元分类问题呢？答案是肯定的。针对多元分类问题，${y=\{0,1,2,...,n \}}$，总共有${n+1}$个类别。其解决思路是，首先将问题转化为二分类问题，即${y=0}$是一个类别，${y=\{ 1,2,...,n\}}$作为另一个类别，然后计算这两个类别的概率；接着，把${y=1}$作为一个类别，把${y=\{ 0,2,...,n\}}$作为另一个类别，再计算这两个类别的概率。由此推广开来，总共需要${n+1}$个预测函数（分类器）。预测出来的概率最高的那个类别，就是样本所属的类别。

3 正则化

回顾线性回归算法梳理，过拟合是指模型很好地拟合了训练样本，但对新样本预测地准确性很差，因为模型太复杂了。解决办法是减少输入特征的个数，或者获取更多的训练样本。正则化也是用来解决模型过拟合问题的一个方法。

• 保留所有的特征，减小特征的权重${\theta_j}$的值。确保所有的特征对预测值都要少量的贡献。
• 当每个特征${x_i}$对预测值${y}$都有少量贡献时，这样的模型可以良好地工作，这就是正则化的目的，可以用它解决特征过多时地过拟合问题。

3.1 线性回归模型正则化

我们先来看线性回归模型的损失函数是如何正则化的：
${J(\theta)=\frac{1}{2m} \lbrack \sum_{i=1}^{m}(h_{\theta}(x^{(i)})-y^{(i)})^2 \rbrack+\lambda \sum_{j=1}^{n} \theta_j^2}$
公式前半部分是线性回归模型的损失函数，后半部分为加入的正则化项。其中${\lambda}$有两个目的，既要维持对训练样本的拟合，又要避免对训练样本的过拟合。如果${\lambda}$值太大，则能确保不出现过拟合，但可能会导致欠拟合。

从数学角度来看，损失函数增加了一个正则项后，损失函数不再唯一地由预测值与真实值的误差所决定了，还和参数${\theta}$的大小有关。有了这个限制之后，要实现损失函数最小的目的，${\theta}$就不能随便取值了。比如某个比较大的${\theta}$值可能会让预测值与真实值的误差${(h_{\theta}(x^{(i)})-y^{(i)})^2}$值很小，但会导致${\theta_j^2}$很大，最终结果就是损失函数太大。这样，通过调节参数${\lambda}$，就可以控制正则项的权重，从而避免线性回归算法过拟合。

梯度下降的时候为什么要对${\theta_j^2}$进行收缩呢？因为加入正则项的损失函数和${\theta_j^2}$成正比，所以迭代的时候需要不断地努力缩小${\theta_j^2}$的值。

3.2 逻辑回归模型正则化

使用相同的思路，我们对逻辑回归模型的损失函数进行正则化，其方法也是在原来的损失函数的基础上加上正则项：
${J(\theta)=- \frac{1}{m} \lbrack \sum_{i=1}^{m} \ y^{(i)} * log(h_\theta(x^{(i)}))+(1-y^{(i)}) * log(1-h_\theta(x^{(i)}))\rbrack + \frac{\lambda}{2m} \sum_{j=1}^{n}\theta_j^2}$
相应地，正则化后的参数迭代公式为：
${\theta_j=\theta_j-\alpha \frac{\partial J(\theta)}{\theta_j}}$
${\theta_j=\theta_j-\alpha \lbrack \frac{1}{m} \sum_{i=1}^{m}(h_\theta(x^{(i)})-x_j^{(i)}) + \frac{\lambda}{m} \theta_j \rbrack}$
需要注意的是，上式中${j \geq 1}$，因为${\theta_0}$没有参与正则化。另外要留意，逻辑回归和线性回归的参数迭代算法看起来形式是一样的，但其实它们的算法不一样，因为两个式子的预测函数${h_\theta(x)}$不一样。针对线性回归${h_\theta(x) = \theta^Tx}$，而针对逻辑回归${h_\theta(x) = \frac{1}{1+e^{-\theta^Tx}}}$。

4 实例：乳腺癌检测

使用逻辑回归算法解决乳腺癌检测问题，我们需要先采集肿瘤病灶造影图片，然后对图片进行分析，从图片中提取特征，再根据特征来训练模型。最终使用模型来检测新采集到的肿瘤病灶造影，以便判断肿瘤是良性还是恶性的。这是一个典型的二分类问题。

4.1 数据采集和特征提取

在工程应用中，数据采集和特征提取工作往往决定着项目的成败。为了简单起见，作为示例，我们直接加载scikit-learn自带的一个乳腺癌数据集。这个数据集是已经采集后的数据：

# 载入数据
from sklearn.datasets import load_breast_cancer

cancer = load_breast_cancer()
X = cancer.data
y = cancer.target
print('data type: {0}; no.positive: {1}; no. negative: {2}'.format(X.shape, y[y==1].shape, y[y==0].shape))
print(cancer.data[0])

"""输出：
data type: (569, 30); no.positive: (357,); no. negative: (212,)
[1.799e+01 1.038e+01 1.228e+02 1.001e+03 1.184e-01 2.776e-01 3.001e-01
 1.471e-01 2.419e-01 7.871e-02 1.095e+00 9.053e-01 8.589e+00 1.534e+02
 6.399e-03 4.904e-02 5.373e-02 1.587e-02 3.003e-02 6.193e-03 2.538e+01
 1.733e+01 1.846e+02 2.019e+03 1.622e-01 6.656e-01 7.119e-01 2.654e-01
 4.601e-01 1.189e-01]
 
 我们可以看到，数据集中共有569个样本，每个样本30个特征，其中357个阳性（y=1）样本，212个阴性样本。

"""

print(cancer.feature_names) # 查看特征名称

"""输出：
['mean radius' 'mean texture' 'mean perimeter' 'mean area'
 'mean smoothness' 'mean compactness' 'mean concavity'
 'mean concave points' 'mean symmetry' 'mean fractal dimension'
 'radius error' 'texture error' 'perimeter error' 'area error'
 'smoothness error' 'compactness error' 'concavity error'
 'concave points error' 'symmetry error' 'fractal dimension error'
 'worst radius' 'worst texture' 'worst perimeter' 'worst area'
 'worst smoothness' 'worst compactness' 'worst concavity'
 'worst concave points' 'worst symmetry' 'worst fractal dimension']
 
"""

为了强调特征提取工作的重要性，这里介绍一下这些特征值的物理含义，我们也可以思考一下，如果让我们来提取特征，会怎么做？
这个数据集总共从病症造影图片中提取了以下10个关键属性。

• radius：半径，即病症中心点离边界的平均距离。
• texture：纹理，灰度值的标准偏差。
• perimeter：周长，即病灶的大小。
• area：面积，也是反映病症大小的一个指标。
• smoothness：平滑度，即半径的变化幅度。
• compactness：密实度，周长的平方除以面积的商，再减1。
• concavity：凹度，凹陷部分轮廓的严重程度。
• concave points：凹点，凹陷轮廓的数量。
• symmetry：对称性。
• fractal dimension：分形维度。

从这些指标里可以看出，有些指标属于“复合”指标，即由其他指标经过运算得到的。比如密实度。
不要小看这种运算构建出来的新特征，这是事物内在逻辑关系的体现。
划重点哈：提取特征时，不妨从事物的内在逻辑关系入手，分析已有特征之间的关系，从而构造出新的特征。

4.2 模型训练

from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression

# 把数据集分为训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
print(X_train.shape, y_train.shape)
print(X_test.shape, y_test.shape)
"""
(455, 30) (455,)
(114, 30) (114,)
"""

# 使用LogisticRegression模型来训练，并计算训练集的评分数据和测试集的评分数据。
# 模型训练
model = LogisticRegression()
model.fit(X_train, y_train)

train_score = model.score(X_train, y_train)
test_score = model.score(X_test, y_test)
print('train score: {train_score: .6f}; test score: {test_score: .6f}'.format(train_score=train_score, test_score=test_score))
"""
train score:  0.958242; test score:  0.947368
"""

# 样本预测
import numpy as np #导入numpy，因为下面要用到它（np。equal()函数）
y_pred = model.predict(X_test)
print('matches: {0}/{1}'.format(np.equal(y_pred, y_test).shape[0], y_test.shape[0]))

"""
matches: 114/114
总共114个测试样本，全部预测正确。这里有个疑问，为什么全部预测正确，而test_score却只有0.973684，而不是1呢？
答案是：scikit-learn不是使用这个测试集数据来计算分数，因为这个数据不能完全反映误差情况，而是使用预测概率数据来计算模型评分。
"""

针对二分类问题，${LogisticRegression}$模型会针对每个样本输出两个概率，即为0的概率和为1的概率，哪个概率高就预测为哪个类别。
我们可以找出针对测试数据集，模型预测的“自信度”低于90%的样本。怎么找出这些样本呢？
我们先计算出测试集里的每个样本的预测概率数据，针对每个样本，它会有两个数据，一是预测其为阳性的概率，另外一个是预测其为阴性的概率；
接着找出预测为阴性的概率大于0.1的样本，然后在结果集里找出预测为阳性的概率也大于0.1的样本，这样就找出了模型预测“自信度”低于90%的样本。
这是因为所有类别的预测概率加起来一定是100%。
我们看一下概率数据：

# 预测概率：找出预测概率低于90%的样本
y_pred_proba = model.predict_proba(X_test) # 计算每个测试样本的预测概率
# 打印第一个样本的数据，以便我们了解数据形式
print('sample of predict probability: {0}'.format(y_pred_proba[0]))

"""输出：
sample of predict probability: [0.00588337 0.99411663]
"""

# 找出第一列，即预测为阴性的概率大于0.1的样本，保存在result里
y_pred_proba_0 = y_pred_proba[:, 0] > 0.1
result = y_pred_proba[y_pred_proba_0]

# 在result结果集里，找到第二列，即预测为阳性的概率大于0.1的样本
y_pred_proba_1 = result[:, 1] > 0.1
print(result[y_pred_proba_1])

"""
[[0.22451742 0.77548258]
 [0.10006953 0.89993047]
 [0.10750749 0.89249251]
 [0.7861653  0.2138347 ]
 [0.16501894 0.83498106]
 [0.16245788 0.83754212]
 [0.40120013 0.59879987]
 [0.54515479 0.45484521]
 [0.80401218 0.19598782]
 [0.34730441 0.65269559]
 [0.27727652 0.72272348]
 [0.27517093 0.72482907]
 [0.65010429 0.34989571]
 [0.23931504 0.76068496]
 [0.31196411 0.68803589]
 [0.75666569 0.24333431]
 [0.65202448 0.34797552]
 [0.43434309 0.56565691]
 [0.40553876 0.59446124]
 [0.40382138 0.59617862]]
"""

我们使用${model.predict\_proba()}$来计算概率，同时找出那些预测“自信度”低于90%的样本。可以看出，最没有把握的样本是[0.54515479 0.45484521]，即只有45.48%的概率是阳性。大家在自己实现这个例子的时候，输出结果可能会有差别，因为训练集和测试集是随机分配的。

模型优化

我们使用${LogisticRegression}$模型的默认参数训练出来的模型，准确性看起来还是挺高的。问题是，有没有优化空间？如果有，往哪个方向优化呢？
我们先尝试增加多项式特征。首先我们使用${Pipeline}$来增加多项式特征：

from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline

# 增加多项式预处理
def polynomial_model(degree=1, **kwarg):
    polynomial_features = PolynomialFeatures(degree=degree, include_bias=False)
    logistic_regression = LogisticRegression(**kwarg)
    pipeline = Pipeline([("polynomial_features", polynomial_features), ("logistic_regression", logistic_regression)])
    return pipeline

# 接着，增加二项多项式特征，创建并训练模型
import time
model = polynomial_model(degree=2, penalty='l1') #我们使用L1范数作为正则项（参数penalty='l1'）

start = time.clock()
model.fit(X_train, y_train)

train_score = model.score(X_train, y_train)
cv_score = model.score(X_test, y_test)
print('elaspe: {0:.6f}; train_score: {1: 0.6f}; cv_score: {2: .6f}'.format(time.clock()-start, train_score, cv_score))

"""输出：
elaspe: 0.369242; train_score:  0.993407; cv_score:  0.964912
"""

可以看到训练集评分和测试集评分都增加了。
为什么使用${L1}$范数作为正则项呢？因为它可以实现参数的稀疏化，即自动帮助我们选择出那些对模型有关联的特征。
我们可以观察一下有多少个特征没有被丢弃，即其对应的模型参数${\theta_j}$非0

logistic_regression = model.named_steps['logistic_regression']
print('model parameters shape: {0}; count of non-zero element: {1}'.format(
    logistic_regression.coef_.shape, np.count_nonzero(logistic_regression.coef_)))

"""输出：
model parameters shape: (1, 495); count of non-zero element: 113
"""

逻辑回归模型的${coef\_}$属性里保存的就是模型参数。从输出结果可以看到，
增加二阶多项式特征后，输入特征由原来的30个增加到了495个，最终大多数特征都被丢弃，只保留了94个有效特征。

到此，总算把逻辑回归算法的思路整理出来了，当初线性回归和逻辑回归这两个“孪生兄弟”把我折磨得欲生欲死，飘飘欲仙。

一点小扩展

我们的预测函数是写成向量形式的：
${h_\theta=g(z)=g(\theta^Tx)=\frac{1}{1+e^{-\theta^Tx}}}$
这个预测函数一次只计算一个样本的预测值，我们要怎样一次性计算出所有样本的预测值呢？
${h=g(X\theta)}$
上述公式即可达到目的。其中${g(x)}$为${sigmoid}$函数。${X}$为${m * n}$的矩阵，即数据集的矩阵表达。
损失函数也有对应的矩阵形式：
${J(\theta)=\frac{1}{m}(-y^T*log(h)-(1-y)^T*log(1-h))}$，其中，${y}$为目标值向量，${h}$为一次性计算出来的所有样本的预测值。