【作者主页】Francek Chen
【专栏介绍】$\lceil$Python机器学习$\rfloor$ 机器学习是一门人工智能的分支学科，通过算法和模型让计算机从数据中学习，进行模型训练和优化，做出预测、分类和决策支持。Python成为机器学习的首选语言，依赖于强大的开源库如Scikit-learn、TensorFlow和PyTorch。本专栏介绍机器学习的相关算法以及基于Python的算法实现。
【GitCode】专栏资源保存在我的GitCode仓库：https://gitcode.com/Morse_Chen/Python_machine_learning。

---

  在介绍了机器学习中相关的基本概念和技巧后，本章我们继续讲解参数化模型中的线性模型。有了前文的基础，我们可以先来对KNN算法和线性回归进行比较，进一步回答“什么是参数化模型”这一问题。对于机器学习算法来说，其目标通常可以抽象为得到某个从输入空间到输出空间的映射 $f:\mathcal{X}\to \mathcal{Y}$，对每个输入数据 $x\in \mathcal{X}$，该映射给出预测的标签 $y=f(x)$。而对于映射的形式，不同算法有不同的假设。像KNN这样，不对$f$的形式做先验假设、在学习中可以得到其任意形式的模型，称为非参数化模型（nonparametric model）。而与之相对，有一类算法会先假设$f$具有某种特定的形式。例如，我们可以假设输入与输出一定满足某个二次函数关系，即 $y=a{x}^2+bx+c$。这时，我们就只需要对映射的参数 $a$、$b$、$c$ 进行学习。像这样对$f$进行先验假设的模型，称为参数化模型（parametric model）。

  相比于非参数化模型，参数化模型由于限制了$f$可能的集合，学习难度相对较低。以上面的二次函数为例，非参数化模型需要在所有可能映射的集合 F = { f : X → Y } \mathcal F=\{f:\mathcal X\to\mathcal Y\} F={f:X→Y} 中寻找合适的$f$，而参数化模型的搜索范围只有 F = { f : X → Y ∣ f ( x ) = a x 2 + b x + c } \mathcal F=\{f:\mathcal X\to\mathcal Y | f(x)=ax^2+bx+c\} F={f:X→Y∣f(x)=ax2+bx+c}，对应于一个3维的实数参数空间 $\{(a,b,c)\in {\mathbb{R}}^3\}$，降低了搜索难度和时间开销。然而，对搜索空间的限制也成为了参数化模型的缺点。如果真实的输入输出关系是 $y=e^x$，那么以二次函数为基础的参数化模型，显然无法在大范围的输入上得到都误差较小的结果。因此，参数化模型通常需要对已知数据和问题特性进行分析，确定合适的参数化假设，才能得到理想的学习结果。除此之外，参数化模型由于假设了数据的分布，其参数的数量通常和数据集的大小无关。因为无论数据集中有多少数据，我们都认为它们是由同一个分布生成的。而非参数化模型也并非“没有参数的模型”，但其参数的数量通常会随数据集的大小而变化，例如KNN算法的参数可以看成是每一个数据的特征和标签，它与数据集始终是同等规模的。

  在线性回归中，我们利用参数化的线性假设解决了回归问题，而分类问题作为机器学习任务中的另一大类别，其与回归问题既有相似也有不同。通常来说，回归问题的输出是连续的，而分类问题的输出是离散的。设输入数据 $x\in {\mathbb{R}}^d$，输出标签 $y\in \mathcal{C}$。其中 $\mathcal{C}\subseteq \mathbb{N}$ 是一有限的离散集合，其每一个元素表示一个不同的类别。事实上，当$\mathcal{C}$是一有限的离散集合时，多分类问题与二分类问题等价。因为我们总可以先判断“$x$是否属于第一类”，再判断“$x$是否属于第二类”，以此类推，从而可以用至多 $|\mathcal{C}|-1$ 次二分类来完成$|\mathcal{C}|$分类。因此，大多数时候我们可以只考虑最简单的二分类问题，其中样本标签 y ∈ { 0 , 1 } y\in\{0,1\} y∈{0,1}。对于分类问题来说，同样可以利用线性假设对问题建模，这一模型就是逻辑斯谛回归（logistic regression），又称为对数几率回归。下面，我们将详细介绍如何利用参数化模型逻辑斯谛回归来处理分类问题。

一、逻辑斯谛函数下的线性模型

  与线性回归类似，对二分类问题，我们同样可以作线性假设。设学习到的映射为 f  ⁣ : R d → { 0 , 1 } f \colon \mathbb{R}^d \to \{0, 1\} f:Rd→{0,1}，参数为$\mathbf{\theta }$。在线性回归中，我们直接计算输入样本与$\mathbf{\theta }$的乘积${\mathbf{\theta }}^{\mathrm{T}}\mathbf{x}$。然而，直接使用该乘积在二分类问题中有两个问题。第一，该乘积是连续的，并不能拟合离散变量；第二，该乘积的取值范围是$\mathbb{R}$，与我们期望的 { 0 , 1 } \{0,1\} {0,1}相距很远。为了解决这两个问题，最简单的方法是再引入阈值$z$，定义$f$如下：
$$f(\mathbf{x})=\{\begin{array}{rlr}& 0& ({\mathbf{\theta }}^{\mathrm{T}}\mathbf{x}\le z)\\ & 1& ({\mathbf{\theta }}^{\mathrm{T}}\mathbf{x}>z)\end{array}$$

  似乎刚刚的两个问题都得到了解决，但是，这一方法又为模型训练带来了新的困难。在线性回归中我们已经介绍过，无论是解析方法还是梯度下降，都需要以函数的梯度为基础。然而，这样的分类方法太“硬”了，使得$f$在阈值处出现了跳跃，从而不再可导；而在阈值之外可导的地方，其导数又始终为$0$。因此，硬分类得到的$f$难以直接训练。

  我们不妨换一个角度来考虑二分类问题。如果把样本$\mathbf{x}$的类别$y$看作是有0和1两种取值的随机变量，我们只需要判断 $P(y=0|\mathbf{x})$ 和 $P(y=1|\mathbf{x})$ 之间大小关系，再将$\mathbf{x}$归为概率较大的一类即可。事实上，我们并不一定需要$f$的输出必须是0或1之中的一个。如果$f$能够给出样本$\mathbf{x}$的类别$y$的概率分布，即 $f(\mathbf{x})=P(y=1|\mathbf{x})$，同样可以达到分类的效果。相比于硬分类，概率分布可以用连续函数建模，从而可以对$f$求梯度。并且，在很多决策问题中，给出每个分类的概率信息比直接给出最后的分类结果要更有用。

  至此我们已经解决了用连续函数给出离散分类和硬分类函数不可导的两个问题。但我们最开始已经提到，${\mathbf{\theta }}^{\mathrm{T}}\mathbf{x}$的取值范围是$\mathbb{R}$，而概率分布的取值范围应当是$[0,1]$。因此，我们需要某种从$\mathbb{R}$到$[0,1]$的映射来确保其取值范围相同。在实践中，我们通常采用逻辑斯谛函数（logistic function）$\sigma$，其定义为$$\sigma (x)=\frac{1}{1+e^{-x}}=\frac{e^x}{e^x+1}$$ 该函数的图像如图1所示。

图1 逻辑斯谛函数的图像

  逻辑斯谛函数是一种常见的sigmoid函数。sigmoid函数是一类形状类似于‘S’型的函数，但在实践中如无特殊说明，一般就指逻辑斯谛函数。逻辑斯谛函数有许多优秀的性质。首先，它关于$(0,0.5)$点对称，从而有 $\sigma (x)+\sigma (-x)=1$，这意味着 $P(y=0|x)=P(y=1|-x)$，即当$x$相反时，其概率分布也正好相反。其次，从图像中可以看出，$\sigma (x)$在$x$偏离0时会迅速收敛到0或1。例如在$x=6$时，$\sigma (x)\approx 0.9975$，与1已经非常接近。这一性质使得其对$x$的变化较为敏感，适合作为分类函数。最后，逻辑斯谛函数在$\mathbb{R}$上连续、单调递增且可导，具有良好的分析性质。其导数为
$$\begin{array}{rl}\frac{\mathrm{d}\sigma (x)}{\mathrm{d}x}& =-\frac{1}{(1+e^{-x}{)}^2}\cdot e^{-x}\cdot (-1)\\ & =\frac{1}{1+e^{-x}}\cdot \frac{e^{-x}}{1+e^{-x}}\\ & =\sigma (x)(1-\sigma (x))\end{array}$$

  综上所述，我们只需要用逻辑斯谛函数对${\mathbf{\theta }}^{\mathrm{T}}\mathbf{x}$进行变换，就可以得到符合要求的映射 $f_{\mathbf{\theta }}(\mathbf{x})=\sigma ({\mathbf{\theta }}^{\mathrm{T}}\mathbf{x})$。虽然这一映射不再是线性的，但它依然是以线性函数${\mathbf{\theta }}^{\mathrm{T}}\mathbf{x}$为基础，再经过某种变换得到的。这样的模型属于广义线性模型，有关广义线性模型可以参考本文的拓展阅读。

二、最大似然估计

  确定了逻辑斯谛回归的数学模型之后，我们接下来还需要确定优化目标。对于有关概率分布的问题，我们常常使用最大似然估计（maximum likelihood estimation，MLE）的思想来优化模型，也即是寻找逻辑斯谛回归的参数$\mathbf{\theta }$使得模型在训练数据上预测出正确标签的概率最大。

  设共有$N$个样本${\mathbf{x}}_1,\dots ,{\mathbf{x}}_N$，类别分别是$y_1,\dots ,y_N\in [0,1]$。对于样本${\mathbf{x}}_i$，如果 $y_i=0$，那么模型预测正确的概率为 $1-f_{\mathbf{\theta }}({\mathbf{x}}_i)$；如果 $y_i=1$，那么概率为$f_{\mathbf{\theta }}({\mathbf{x}}_i)$。将两者综合起来，可以得到模型正确的概率为 $f_{\mathbf{\theta }}({\mathbf{x}}_i{)}^{y_i}(1-f_{\mathbf{\theta }}({\mathbf{x}}_i){)}^{1-y_i}$。假设样本之间是两两独立的，那么模型将所有样本的分类都预测正确的概率就等于单个样本概率的乘积：$$L(\mathbf{\theta })=\prod _{i=1}^N{f}_{\mathbf{\theta }}({\mathbf{x}}_i{)}^{y_i}(1-f_{\mathbf{\theta }}({\mathbf{x}}_i){)}^{1-y_i}$$ 该函数也称为似然函数（likelihood function）。为了使模型的预测尽可能准确，我们需要寻找使似然函数最大的参数$\mathbf{\theta }$。但是，该函数的连乘形式使得求导和优化都很困难，在计算机上直接计算甚至很容易造成浮点数越界。因此，我们一般将似然值取对数，也即是对数似然（log-likelihood），将连乘转化为求和：$$l(\mathbf{\theta })=\log L(\mathbf{\theta })=\sum _{i=1}^N\left[y_i\log f_{\mathbf{\theta }}({\mathbf{x}}_i)+(1-y_i)\log (1-f_{\mathbf{\theta }}({\mathbf{x}}_i))\right]$$ 由于对数函数是单调递增的，优化$l(\mathbf{\theta })$和$L(\mathbf{\theta })$可以得到相同的结果。于是，我们的优化目标为$$\underset{\mathbf{\theta }}{\max }l(\mathbf{\theta })$$ 对$l(\mathbf{\theta })$求梯度，得到
$$\begin{array}{rl}\nabla l(\mathbf{\theta })& =\sum _{i=1}^N\left[\frac{\nabla f_{\mathbf{\theta }}}{f_{\mathbf{\theta }}}{y}_i-\frac{\nabla f_{\mathbf{\theta }}}{1-f_{\mathbf{\theta }}}(1-y_i)\right]\end{array}$$ 再把 $f_{\mathbf{\theta }}(\mathbf{x})=\sigma ({\mathbf{\theta }}^{\mathrm{T}}\mathbf{x})$ 代入，并利用 ${\nabla }_{\mathbf{\theta }}\sigma ({\mathbf{\theta }}^{\mathrm{T}}\mathbf{x})=\sigma ({\mathbf{\theta }}^{\mathrm{T}}\mathbf{x})(1-\sigma ({\mathbf{\theta }}^{\mathrm{T}}\mathbf{x}))\mathbf{x}$ 就得到
$$\begin{array}{rl}\nabla l(\mathbf{\theta })& =\sum _{i=1}^N\left[(1-\sigma ({\mathbf{\theta }}^{\mathrm{T}}{\mathbf{x}}_i))y_i{\mathbf{x}}_i-\sigma ({\mathbf{\theta }}^{\mathrm{T}}{\mathbf{x}}_i)(1-y_i){\mathbf{x}}_i\right]\\ & =\sum _{i=1}^N(y_i-\sigma ({\mathbf{\theta }}^{\mathrm{T}}{\mathbf{x}}_i)){\mathbf{x}}_i\end{array}$$

  注意到我们的优化目标是最大化$l$，因此定义损失函数 $J(\mathbf{\theta })=-l(\mathbf{\theta })$，将优化目标转化为最小化$J(\mathbf{\theta })$。根据梯度下降算法，设学习率为$\eta$，则$\mathbf{\theta }$的更新公式为
$$\begin{array}{rl}\mathbf{\theta }& \leftarrow \mathbf{\theta }+\eta \nabla l(\mathbf{\theta })\\ & =\mathbf{\theta }+\eta \sum _{i=1}^N(y_i-\sigma ({\mathbf{\theta }}^{\mathrm{T}}{\mathbf{x}}_i)){\mathbf{x}}_i\end{array}$$

  需要注意，在实际计算时，我们通常会对样本取平均，来消除样本规模对学习率$\eta$的影响。为了得到更简洁和易于计算的形式，定义样本矩阵 $\mathbf{X}=({\mathbf{x}}_1,\dots ,{\mathbf{x}}_N{)}^{\mathrm{T}}$，标签向量 $\mathbf{y}=(y_1,\dots ,y_N{)}^{\mathrm{T}}$，以及向量形式的逻辑斯谛函数 $\sigma (\mathbf{u})=(\sigma (u_1),\cdots ,\sigma (u_N){)}^{\mathrm{T}}$，就可以将梯度写成矩阵形式：$$\nabla J(\mathbf{\theta })=-\nabla l(\mathbf{\theta })={\mathbf{X}}^{\mathrm{T}}(\mathbf{y}-\sigma (\mathbf{X}\mathbf{\theta }))$$ 从而梯度下降算法的矩阵形式为$$\mathbf{\theta }\leftarrow \mathbf{\theta }+\eta {\mathbf{X}}^{\mathrm{T}}(\mathbf{y}-\sigma (\mathbf{X}\mathbf{\theta }))$$

  除此之外，我们再为模型加入$L_2$正则化约束，设正则系数为$\lambda$，则完整的优化目标与迭代公式为$$\underset{\mathbf{\theta }}{\min }J(\mathbf{\theta })=\underset{\mathbf{\theta }}{\min }-l(\mathbf{\theta })+\frac{\lambda }{2}\Vert \mathbf{\theta }{\Vert }_2^2$$$$\nabla J(\mathbf{\theta })=-{\mathbf{X}}^{\mathrm{T}}(\mathbf{y}-\sigma (\mathbf{X}\mathbf{\theta }))+\lambda \mathbf{\theta }$$$$\mathbf{\theta }\leftarrow (1-\lambda \eta )\mathbf{\theta }+\eta {\mathbf{X}}^{\mathrm{T}}(\mathbf{y}-\sigma (\mathbf{X}\mathbf{\theta }))$$

  逻辑斯谛回归的损失函数也是凸函数，图2展示了在逻辑斯谛回归模型 $f_{\mathbf{\theta }}(\mathbf{x})=\sigma ({\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2)$ 中，带有正则化约束的损失函数等值线。为了能在二维平面中展示参数的变化，${\theta }_0$已经固定，仅有${\theta }_1$和${\theta }_2$变化。图中绿色的曲线表示从不同初始参数开始进行梯度下降的轨迹。这些等值线已经不是椭圆形，但进行梯度下降时，凸函数的性质仍然保证$\theta$可以收敛到最优解。

图2 逻辑斯谛回归损失函数的梯度下降

  由于逻辑斯谛函数只适用于二分类的情况，对于多分类任务，我们需要寻找新的映射函数。这时，我们通常采用柔性最大值（softmax）函数来将线性部分的预测${\mathbf{\theta }}^{\mathrm{T}}\mathbf{x}$映射到多分类概率上。该函数的定义为$$\sigma (\mathbf{z}{)}_i=\frac{e^{z_i}}{\sum _{j=1}^K{e}^{z_j}}$$ 其中，$\mathbf{z}$是一个$K$维的得分向量，每一维取值 $z_j\in (-\infty ,+\infty )$ 代表第$j$类的得分。softmax函数将$K$维的向量映射到新的$K$维向量$\sigma (\mathbf{z})$，容易验证，映射后的向量各个元素取值在$(0,1)$，并且它们取值之和为1。因此，$\sigma (\mathbf{z})$可以看作一个有$K$种取值的离散随机变量的概率分布。这一特性使得softmax函数可以在多分类问题上作为映射函数。设参数${\mathbf{\theta }}_1,\dots ,{\mathbf{\theta }}_K$分别用来预测样本$\mathbf{x}$属于第1类到第$K$类的未归一化概率，记 $\mathbf{z}=({\mathbf{\theta }}_1^{\mathrm{T}}\mathbf{x},\dots ,{\mathbf{\theta }}_K^{\mathrm{T}}\mathbf{x})$。那么，$\sigma (\mathbf{z}{)}_i$就给出了$\mathbf{x}$属于第$i$类的概率。

  下面我们来考察softmax函数的一些性质。直观上看，在原始的向量$\mathbf{z}$中越大的元素，映射后仍然越大。因此，softmax函数不改变向量元素的大小顺序。此外，与逻辑斯谛函数类似，softmax函数避免了直接取最大元素导致的不可导问题。事实上，如果在上式中令 $K=2$，softmax函数就退化为逻辑斯谛函数。或者说，逻辑斯谛函数只是softmax函数在二分类前提下的特殊情况。所以，我们同样用$\sigma$来表示softmax函数。

  简单起见，在本文的后续部分我们仍然考虑二分类任务。

三、分类问题的评价指标

  在动手实现逻辑斯谛回归之前，我们还需要解决一个问题：如何判断分类模型的好坏呢？对于一个二分类问题，最简单的评估模型好坏的方式就是测试模型的正确率。这时，无论是用直接给出类别的硬分类，还是给出不同分类概率的软分类，我们都需要考虑分类的阈值。例如，如果模型判断某样本的类别是正类 $y=1$ 的概率是0.6，而预设的阈值是0.5，就应当将该样本归为正类；如果阈值是0.7，虽然该样本是正类的概率比是负类的概率大，我们仍然应当将它归为负类。在我们的默认思维中，一般会通过比较正负类的概率大小给出最后的判断，这实际上隐含了“阈值为0.5”的假设。

  大家或许有疑问，为什么不能遵循最自然的大小关系，将0.5作为阈值呢？这与实际应用中分类任务的特点有关。例如，当前的任务是进行汽车的出厂检测，根据检测的数据判断汽车是否有质量问题，没有问题是负类，有问题是正类。在这样的任务中，如果模型把有问题的汽车错判为没有问题，就可能让质量不合格的汽车出场，其后果比把没有问题的汽车错判为有问题要严重得多。因此，我们需要把正类的阈值设置的很高，比如0.99。可以发现，虽然抽象的分类问题中，正类负类之间的地位是对等的，因此将正类错判为负类和将负类错判为正类没有什么区别；但是实际场景中，这两者往往并不对称，我们必须根据任务的特点设置合适的阈值。并且，只有正确率这一标准是不够的，还需要分别考虑正负类分别被错判的比例。

  在机器学习中，我们通常用图3的混淆矩阵（confusion matrix）来统计不同的分类结果。其中，真阳性（true positive，TP）表示将真实的正类判别为正类，真阴性（true negative，TN）表示将真实的负类判别为负类，这两者属于判断正确的结果。假阴性（false negative，FN）表示将真实的正类判别为负类，假阳性（false positive，FP）表示将真实的负类判别为正类，属于判断错误的结果。我们常说的精度（accuracy），也称准确率，就是判断正确的样本数量除以样本总数的比例：$$\mathrm{Acc}=\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{FP}+\mathrm{FN}+\mathrm{TN}}$$

  在上面汽车质检的例子中，我们更关心找出有问题的汽车占真正有问题汽车的比例，即真阳性率（true positive rate），也称查全率或召回率（recall）：$$\mathrm{Rec}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}$$

  此外，我们常常还希望假阳性尽可能少，即真阳性占模型判为阳性的样本比例尽可能高。这一指标称为查准率（precision），也称精确率：$$\mathrm{Prec}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}$$

  由于样本总数固定，混淆矩阵的4个指标并非相互独立。因此，我们常用查全率（召回率）和查准率（精确率）这两个标准作为代表，其他指标的变化也可以通过它们来间接反映。

图3 混淆矩阵的组成

  或许大家已经发现，混淆矩阵导出的数个统计指标是不可兼得的，它们都与预设的阈值有关。例如，如果我们将阈值设得很低，就会将大量样本判为正类，提高真阳性的数量，从而使召回率上升。然而，这种做法同时会使许多实际上为负类的样本被归为正类，造成大量假阳性，反而可能让精确率下降。因此，我们应当根据任务具体的特点，先确定不同指标的重要程度，再选取合适的阈值。为了同时反映各个指标的情况，达到平衡点，我们引入新的指标F1分数（F1 score），定义为精确率和召回率的调和平均值：$$F_1=\frac{2}{1/\mathrm{Prec}+1/\mathrm{Rec}}=\frac{2\times \mathrm{Prec}\times \mathrm{Rec}}{\mathrm{Prec}+\mathrm{Rec}}$$

  可以看出，F1分数最小为0，在精确率或召回率有一方为0时取到；最大为1，在精确率和召回率相等时取到。因此，单方面增大或减小精确率和召回率中的一个无法使F1分数增高，必须要让两者得到平衡。如果在实际问题中，精确率和召回率的重要程度不同，我们还可以把F1分数扩展为F-beta分数：$$F_{\beta }=(1+{\beta }^2)\times \frac{\mathrm{Prec}\times \mathrm{Rec}}{{\beta }^2\times \mathrm{Prec}+\mathrm{Rec}}$$

  上式相当于给精确率添加了权重 $1/(1+{\beta }^2)$，给召回率添加了权重 ${\beta }^2/(1+{\beta }^2)$。$\beta$越大，召回率就越重要；反之，$\beta$越小，精确率就越重要。从数学角度来看，当$F_{\beta }$的定义式中 $\beta =0$ 时，它就退化为精确率；当 $\beta \to \infty$ 时，它就退化为召回率。在实践中，$F_2$和$F_{0.5}$都是较为常用的指标。

  对于更普遍的情况来说，我们希望衡量模型在不同阈值下的整体表现，根据不同指标随阈值的变化来得出与阈值无关的模型本身的特性。因此，我们选择假阳性率（false positive rate，FPR）与真阳性率（召回率）这两个随阈值变化趋势相同的指标，把它们的值绘制成曲线。其中，FPR定义为$$\mathrm{FPR}=\frac{\mathrm{FP}}{\mathrm{FP}+\mathrm{TN}}$$

  我们来具体分析一下FPR与TPR随阈值的变化趋势。当阈值为0时，模型会把所有样本归为正类。这时不存在假阴性和真阴性，因此FPR和TPR都为1。当阈值逐渐增大，模型将把越来越多的样本归为负类，所以真阳性和假阳性减少，假阴性和真阴性增多，从而FPR和TPR都减小。而阈值增大到1时，所有样本都被归为负类，情况与阈值为0时正好相反，不存在真阳性和假阳性，FPR和TPR都减小到0。因此，这两个统计指标都随阈值的增大而减小。

  那么，它们的变化趋势为何能反映模型的好坏呢？为了方便理解，我们在图4中绘制了TPR随FPR的变化曲线，该曲线称为受试者操作特征（receiver operating characteristic，ROC）曲线。从右上角到左下角，阈值由0增大到1。需要注意，这两个指标并非同步变化的。例如，某样本真实类别为负类，模型判断它是正类的概率为0.305。并且阈值从0.3增大到0.31的过程中，其他样本的类别预测没有变化，仅有这一个样本由正类变为负类。那么在这个过程中，真阳性和假阴性没有变化，假阳性减少，真阴性增加，从而TPR不变、FPR减小。该过程反映在ROC上是一条水平线。相反，如果上例的样本真实类别是正类，那么真阳性减少、假阴性增加、假阳性和真阴性不变，从而TPR减小，FPR不变，在ROC线上是一条竖直线。由于数据集是离散的，而模型预测的概率值连续，绝大多数情况下不会有两个样本的预测值相同。因此，ROC曲线一般是交替由水平线和竖直线组成，呈阶梯状，仅当不少于两个且类别不同的样本被模型赋予了完全相同预测值，ROC曲线才会出现斜线的部分。

图4 ROC与AUC示意

  对于模型来说，我们期望它能将尽可能多的样本分类正确、同时还要少出错，也即TPR尽可能高、FPR尽可能低。最理想的情况下，当阈值合适的时候，所有的正样本都归为正类，所有的负样本都归为负类，这时 $\mathrm{FPR}=0$，$\mathrm{TPR}=1$，在图中是左上角的顶点。由于ROC曲线的形状是固定的，该曲线必然只有 $(0,0)-(0,1)-(1,1)$ 的两段折线。而一般来说，由于数据集中噪声和模型自身能力的限制，虽然不一定存在完美分割的阈值，但我们仍然希望ROC曲线能尽可能偏向左上方。图4中还有一条从左下到右上的黑色虚线，它表示模型完全随机预测时的ROC。平均意义上，由于每个样本在任何阈值下被归为正类和负类的概率都是0.5，FPR和TPR将同步变化，形成了这条直线。

  “曲线偏向左上方”是一个定性的描述，为了定量衡量ROC曲线表示的模型好坏，我们通常计算ROC曲线与$x$轴和直线 $x=1$ 围成的面积，称为曲线下面积（area under the curve，AUC）。最好情况下，ROC经过左上角，AUC为1；随机情况下AUC是黑色虚线下的三角形面积，为0.5；如果模型比随机预测还要差，AUC会小于0.5，但这种情况事实上不会发生。

  AUC除了ROC下的面积，还有另一层含义。我们以一个包含5个样本的数据集来说明，这些样本依次记为 $n_1,n_2,p_1,p_2,p_3$，其中 $n_1,n_2$ 是负类，$p_1,p_2,p_3$ 是正类。我们将样本按模型预测的正类概率从小到大排序。理想情况下，存在将正负类完美分隔的阈值，即所有正类的概率大于所有负类的概率，如表1所示。

表 1 可以完美分隔的情况

$n_1$	$n_2$	$p_1$	$p_2$	$p_3$
$0.1$	$0.2$	$0.6$	$0.7$	$0.8$

  显然，如果我们任取一个负类、任取一个正类，负类排在正类前面的概率是$1$，此时AUC也是$1$。假如模型并不理想，给出的概率如表2所示。

表2 无法完美分隔的情况

$n_1$	$p_1$	$n_2$	$p_2$	$p_3$
$0.1$	$0.5$	$0.6$	$0.7$	$0.8$

同样任取一个负类$n$、任取一个正类$p$，这时负类排在正类前面的概率为$$P=P(n=n_1)P(p=p_1,p_2,p_3)+P(n=n_2)P(p=p_2,p_3)=\frac{1}{2}\times 1+\frac{1}{2}\times \frac{2}{3}=\frac{5}{6}$$ 再来计算此时的AUC。绘制出形如图5的ROC曲线，容易算出AUC为$$\mathrm{AUC}=\frac{2}{3}\times \frac{1}{2}+1\times \frac{1}{2}=\frac{5}{6}$$

图5 表2对应的ROC曲线

  与上面得到的负类排在正类前面的概率相等！事实上，这一结论对于更多样本的场合也是成立的。从这里可以看出，模型将真实的正类和负类分得越清楚，就有越多的负类排在正类前面，AUC就越大。如果模型分类完全随机，那么任取一个负类和一个正类，它们的位置也是随机的，负类排在正类前面的概率为0.5，所以得到的AUC也是0.5，同样与上面的结论相符。这一结论告诉我们，虽然我们是从阈值的变化中推导出ROC与AUC，但是AUC的值事实上与阈值是的选取无关的，只与模型本身对正负类的预测结果有关。在下一节的动手应用中，我们就用本节讲述的各种标准来作为模型的评价指标。

  此外需要提醒的是，本节介绍的分类指标大都是针对二分类问题的。对于多分类问题，准确率（$\mathrm{Acc}$）仍然是一个最直接的分类指标，也即是样本预测类别和真实类别相同的概率。而召回率、精确率和F1分数皆为特定分类而计算的，例如在二分类中它们是针对正例而计算的。因此，在多分类任务中，我们可以针对关注的类别来计算召回率、精确率和F1分数。更进一步，我们可以对所有类别的F1分数做平均，从而得到该多分类任务的一个总体的F1分数。具体来说，如果是对所有类别的F1分数做直接平均，那么得到的最终指标称为宏观F1分数（macro-F1 score）；如果是先计算所有类别总的TP、FP和FN，再由此计算精确率、召回率和F1分数，那么得到的最终指标称为微观F1分数（micro-F1 score）。

四、动手实现逻辑斯谛回归

  下面，我们动手实现梯度下降算法求解逻辑斯谛回归。本文所用的数据集lr_dataset.csv包含了二维平面上的一些点，这些点按位置的不同分为两类。表3展示了数据集中的几条样本，每条样本依次包含横坐标、纵坐标和类别标签。我们的任务是训练逻辑斯谛回归模型，使模型对点的类别预测尽可能准确。

表3 逻辑斯谛回归数据集中的样本示例

横坐标	纵坐标	标签
0.4304	0.2055	1
0.0898	-0.1527	1
-0.8257	-0.9596	0

  首先，我们读入并处理数据集，并将其在平面上的分布展示出来。图中的每个红色圆圈代表一个正类，蓝色叉号代表一个负类。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator

# 从源文件中读入数据并处理
lines = np.loadtxt('lr_dataset.csv', delimiter=',', dtype=float)
x_total = lines[:, 0:2]
y_total = lines[:, 2]
print('数据集大小:', len(x_total))

# 将得到的数据在二维平面上制图，不同的类别染上不同的颜色以便于观察样本点的分布
pos_index = np.where(y_total == 1)
neg_index = np.where(y_total == 0)
plt.scatter(x_total[pos_index, 0], x_total[pos_index, 1], marker='o', color='coral', s=10)
plt.scatter(x_total[neg_index, 0], x_total[neg_index, 1], marker='x', color='blue', s=10)
plt.xlabel('X1 axis')
plt.ylabel('X2 axis')
plt.show()

# 划分训练集与测试集
np.random.seed(0)
ratio = 0.7
split = int(len(x_total) * ratio)
idx = np.random.permutation(len(x_total))
x_total = x_total[idx]
y_total = y_total[idx]
x_train, y_train = x_total[:split], y_total[:split]
x_test, y_test = x_total[split:], y_total[split:]

  接下来，我们实现上节中所描述的几个评价指标。这里选用准确率和AUC两种。对于AUC以外的指标，我们需要样本的真实标签以及某个阈值下模型的预测标签。而AUC中已经包含了阈值变化，所以需要的是模型原始的预测概率。

def acc(y_true, y_pred):
    return np.mean(y_true == y_pred)

def auc(y_true, y_pred):
    # 按预测值从大到小排序，越靠前的样本预测正类概率越大
    idx = np.argsort(y_pred)[::-1]
    y_true = y_true[idx]
    y_pred = y_pred[idx]
    # 把y_pred中不重复的值当作阈值，依次计算FP样本和TP样本数量
    # 由于两个数组已经排序且位置对应，直接从前向后累加即可
    tp = np.cumsum(y_true) 
    fp = np.cumsum(1 - y_true)
    tpr = tp / tp[-1]
    fpr = fp / fp[-1]
    # 依次枚举FPR，计算曲线下的面积
    # 方便起见，给FPR和TPR最开始添加(0,0)
    s = 0.0
    tpr = np.concatenate([[0], tpr])
    fpr = np.concatenate([[0], fpr])
    for i in range(1, len(fpr)):
        s += (fpr[i] - fpr[i - 1]) * tpr[i]
    return s

  由于本节所用的数据集大小较小，我们不再每次取小批量进行迭代，而是直接用完整的训练集计算梯度。接下来，我们按照之前推导的梯度下降公式定义训练函数，并设置学习率和迭代次数，查看训练结果。除了训练曲线之外，我们还将模型参数$\mathbf{\theta }$表示的直线 ${\mathbf{\theta }}^{\mathrm{T}}\mathbf{x}=0$ 在平面上展示出来。在直线同一侧的点被模型判断成同一类。

# 逻辑斯谛函数
def logistic(z):
    return 1 / (1 + np.exp(-z))

def GD(num_steps, learning_rate, l2_coef):
    # 初始化模型参数
    theta = np.random.normal(size=(X.shape[1],))
    train_losses = []
    test_losses = []
    train_acc = []
    test_acc = []
    train_auc = []
    test_auc = []
    for i in range(num_steps):
        pred = logistic(X @ theta)
        grad = -X.T @ (y_train - pred) + l2_coef * theta
        theta -= learning_rate * grad
        # 记录损失函数
        train_loss = - y_train.T @ np.log(pred) \
                     - (1 - y_train).T @ np.log(1 - pred) \
                     + l2_coef * np.linalg.norm(theta) ** 2 / 2
        train_losses.append(train_loss / len(X))
        test_pred = logistic(X_test @ theta)
        test_loss = - y_test.T @ np.log(test_pred) \
                    - (1 - y_test).T @ np.log(1 - test_pred)
        test_losses.append(test_loss / len(X_test))
        # 记录各个评价指标，阈值采用0.5
        train_acc.append(acc(y_train, pred >= 0.5))
        test_acc.append(acc(y_test, test_pred >= 0.5))
        train_auc.append(auc(y_train, pred))
        test_auc.append(auc(y_test, test_pred))
        
    return theta, train_losses, test_losses, train_acc, test_acc, train_auc, test_auc

  最后，我们定义各个超参数，进行训练并绘制出训练损失、准确率和AUC随训练轮数的变化曲线，以及模型学到的分割直线。

# 定义梯度下降迭代的次数，学习率，以及L2正则系数
num_steps = 250
learning_rate = 0.002
l2_coef = 1.0
np.random.seed(0)

# 在x矩阵上拼接1
X = np.concatenate([x_train, np.ones((x_train.shape[0], 1))], axis=1)
X_test = np.concatenate([x_test, np.ones((x_test.shape[0], 1))], axis=1) 

theta, train_losses, test_losses, train_acc, test_acc, train_auc, test_auc = GD(num_steps, learning_rate, l2_coef)

# 计算测试集上的预测准确率
y_pred = np.where(logistic(X_test @ theta) >= 0.5, 1, 0)
final_acc = acc(y_test, y_pred)
print('预测准确率：', final_acc)
print('回归系数：', theta)

plt.figure(figsize=(13, 9))
xticks = np.arange(num_steps) + 1
# 绘制训练曲线
plt.subplot(221)
plt.plot(xticks, train_losses, color='blue', label='train loss')
plt.plot(xticks, test_losses, color='red', ls='--', label='test loss')
plt.gca().xaxis.set_major_locator(MaxNLocator(integer=True))
plt.xlabel('Epochs')
plt.ylabel('Loss')
plt.legend()

# 绘制准确率
plt.subplot(222)
plt.plot(xticks, train_acc, color='blue', label='train accuracy')
plt.plot(xticks, test_acc, color='red', ls='--', label='test accuracy')
plt.gca().xaxis.set_major_locator(MaxNLocator(integer=True))
plt.xlabel('Epochs')
plt.ylabel('Accuracy')
plt.legend()

# 绘制AUC
plt.subplot(223)
plt.plot(xticks, train_auc, color='blue', label='train AUC')
plt.plot(xticks, test_auc, color='red', ls='--', label='test AUC')
plt.gca().xaxis.set_major_locator(MaxNLocator(integer=True))
plt.xlabel('Epochs')
plt.ylabel('AUC')
plt.legend()

# 绘制模型学到的分隔直线
plt.subplot(224)
plot_x = np.linspace(-1.1, 1.1, 100)
# 直线方程：theta_0 * x_1 + theta_1 * x_2 + theta_2 = 0
plot_y = -(theta[0] * plot_x + theta[2]) / theta[1]
pos_index = np.where(y_total == 1)
neg_index = np.where(y_total == 0)
plt.scatter(x_total[pos_index, 0], x_total[pos_index, 1], marker='o', color='coral', s=10)
plt.scatter(x_total[neg_index, 0], x_total[neg_index, 1], marker='x', color='blue', s=10)
plt.plot(plot_x, plot_y, ls='-.', color='green')
plt.xlim(-1.1, 1.1)
plt.ylim(-1.1, 1.1)
plt.xlabel('X1 axis')
plt.ylabel('X2 axis')
plt.savefig('output_16_1.png')
plt.savefig('output_16_1.pdf')
plt.show()

  为了进一步展示模型学到的概率信息，我们把模型预测的每个数据点的概率展示在图6(a)中，其中的数据点的坐标是$(x_1,x_2)$，纵轴是经过逻辑斯谛函数后模型预测的概率，绿色的平面是把分隔直线从二维伸展到三维的结果。样本点越靠上，说明模型预测该样本点属于正类的概率越大。可以看出，这些点分布出的曲面很像空间中逻辑斯谛函数的图像。事实上，如果沿着分隔平面的视角看过去，这些点恰好会组成逻辑斯谛函数的曲线，如图6(b)所示。这是因为逻辑斯谛回归模型给出的分类结果本质上还是概率结果，在分隔平面附近的点模型很难判断它们是正类还是负类，给出的概率接近0.5；远离分隔平面的点模型比较有信心，给出的概率接近0或者1。由于逻辑斯谛回归中线性预测部分${\mathbf{\theta }}^{\mathrm{T}}\mathbf{x}$是由逻辑斯谛函数$\sigma$映射为概率的，这些点的概率在空间中自然也就符合逻辑斯谛函数的曲线。

图6 逻辑斯谛回归预测值的三维示意

五、使用sklearn中的逻辑斯谛回归模型

  使用scikit-learn库中linear_model模块的LogisticRegression类可以建立逻辑回归模型，其语法格式和参数说明如下。

sklearn.linear_model.LogisticRegression(penalty='l2', dual=False, tol=0.0001, C=1.0, fit_intercept=True, intercept_scaling=1, class_weight=None, random_state=None, solver='liblinear', max_iter=100, multi_class='ovr', verbose=0, warm_start=False, n_jobs=1)

参数名称	说明
penalty	接收str。表示正则化选择参数，可选l1或l2。默认为l2。‘l1’：L1 正则化(Lasso)，可以用于特征选择；‘l2’：L2 正则化(Ridge)，在处理多重共线性时通常更为稳定；‘none’：不使用正则化。
dual	接收bool，默认值=False。是否使用对偶形式的优化。在penalty='l2’和solver='liblinear’时可以使用。对于n_samples > n_features的数据集，通常dual=False是更有效的。
tol	接收float，默认值=1e-4。停止准则的容忍度。优化算法在收敛时的容忍度，较低的值意味着更严格的停止条件。
C	接收float，默认值=1.0。正则化的倒数。较小的值表示较强的正则化，较大的值表示较弱的正则化。C的值越大，正则化的效果越弱。
fit_intercept	接收bool，默认值=True。是否计算截距（偏置）。如果设置为False，则模型将不会包括截距，假定数据已中心化。
intercept_scaling	接收float，默认值=1。当solver='liblinear’和fit_intercept=True时，拯救截距的缩放系数。只有当fit_intercept=True时有效。此参数在solver='liblinear’时特别有用。
class_weight	接收dict, list, str，默认值=None。类别权重。可以是dict类型，指定每个类别的权重，或是‘balanced’，使类别权重与样本频率成反比，或者None，表示不使用权重。
solver	接收str，默认值=‘lbfgs’。求解器类型，用于优化。可选值包括：‘newton-cg’：牛顿共轭梯度法；‘lbfgs’：拟牛顿法(常用且稳定)；‘liblinear’：线性分类器；‘saga’：随机平均梯度下降法(适合大规模数据集)。
max_iter	接收int，默认值=100。最大迭代次数，用于求解优化算法的迭代过程。
multi_class	接收str，默认值=‘auto’。多分类策略。可选值包括：‘auto’：自动选择，通常会选择‘ovr’；‘ovr’(one-vs-rest)：将多分类问题拆分为多个二分类问题；‘multinomial’：多项式分类，适用于整个类别的情况。
verbose	接收int，默认值=0。决定是否输出详细信息，调试时可能有用。
warm_start	接收bool，默认值=False。是否使用上次调用fit方法的结果来初始化当前模型，从而加快计算速度。
n_jobs	接收int，默认值=1。用于设置计算时的并行作业数。1表示单线程，-1表示使用所有可用的处理器核心。

  与线性回归相似，sklearn同样提供了封装好的逻辑斯谛回归模型LogisticRegression。我们直接使用该工具求解逻辑斯谛回归问题，并与自己实现的梯度下降方法进行比较。sklearn中的方法默认添加了正则系数 $\lambda =1.0$ 的$L_2$约束，与我们的参数设置相同。最终得到的回归系数与我们动手实现得到的结果基本是一致的。

from sklearn.linear_model import LogisticRegression

# 使用线性模型中的逻辑回归模型在数据集上训练
# 其提供的liblinear优化算法适合在较小数据集上使用
# 默认使用系数为1.0的L2正则化约束
# 其他可选参数请参考官方文档
lr_clf = LogisticRegression(solver='liblinear')
lr_clf.fit(x_train, y_train)
print('回归系数：', lr_clf.coef_[0], lr_clf.intercept_)

# 在数据集上用计算得到的逻辑回归模型进行预测并计算准确度
y_pred = lr_clf.predict(x_test)
print('准确率为：',np.mean(y_pred == y_test))

六、交叉熵与最大似然估计

  在训练逻辑斯谛回归模型时，我们从概率分布出发，采用最大似然估计的思想得到了损失函数。而从信息论的角度，我们也能得到相同的结果。在信息论中，当一个随机事件发生时，它就会提供一定的信息。而事件发生概率越小，其发生时所提供的信息量也就越大。例如，连续抛一枚硬币$100$次，如果出现正面和反面的次数接近，那么这很正常；如果连续出现$100$次正面，那么我们就不免怀疑硬币上是否有什么机关。用数学语言描述，设事件$X_i$发生的概率为$P(X_i)$，那么$X_i$发生所能提供的信息是$$I(X_i)=-\log P(X_i)$$ 从上式中可以看出，确定事件发生不会提供任何信息，这也符合我们的直观感受。而当许多事件互相影响时，我们还需要对这些事件整体的性质进行研究。

  下面，我们只考虑事件离散且有限的情况。设有$n$个事件$X_1,\dots ,X_n$，其发生的概率分别为 $P(X_1),\dots ,P(X_n)$，满足 $\sum _{i=1}^{n}P(X_i)=1$，且任意两个事件都互斥，即 $\forall i\ne j,P(X_i\cap X_j)=0$。我们可以用一个随机变量$X$来表示这些事件，$X=i$ 表示事件$X_i$发生，并用 $p(x)=P(X=x)$ 来表示这些事件的概率分布。在每一时刻，这些事件中有且仅有一个会发生。可以发现，预测每一时刻发生事件的难度取决于分布的整体性质。如果该分布中有某个事件发生的概率很大，预测的难度就较低；反过来，如果各个事件发生的概率都很接近，那么就很难预测到底哪一个事件会发生，也就是分布的不确定性更大。例如，我们可以几乎确定太阳明天会从东边升起，却无法预测抛一枚均质硬币会得到正面还是反面，因为两者出现的概率几乎都是$1/2$。模仿物理学中衡量系统无序程度的熵（entropy）的概念，在信息论中，我们也用熵来衡量分布的不确定程度。上述分布的熵$H(p)$定义为$$H(p)={\mathbb{E}}_{X\sim p(x)}[I(X)]=\sum _{i=1}^{n}P(X_i)I(X_i)=-\sum _{i=1}^{n}P(X_i)\log P(X_i)$$ 经过一些数学推导可以得到，当某个事件发生的概率为$1$、其他事件发生概率为$0$时，分布的熵最小，为 $H=0$；当所有事件发生的概率都相等，即 $P(X_i)=1/n$ 时，分布的熵最大，为 $H=\log n$。

  更进一步，如果关于随机变量$X$存在两个概率分布$p(x)$和$q(x)$，我们可以用相对熵（relative entropy）来衡量这两个分布的距离。相对熵又称为库尔贝克-莱布勒散度（Kullback-Leibler divergence），简称KL散度，其定义为$$D_{\mathrm{KL}}(p\Vert q)={\mathbb{E}}_{X\sim p(x)}\left[\log \frac{p(x)}{q(x)}\right]$$

  KL散度是一种比较特殊的距离度量。我们知道，判断两个数$a$和$b$的关系，既可以计算它们的差 $a-b$，也可以计算它们的比值$a/b$。$a/b>1$ 说明$a$较大，反之同理，而比值越接近1，则说明$a$与$b$越接近。我们观察它的定义式，其期望的内部是$\begin{array}{r}\log \frac{p(x)}{q(x)}\end{array}$，是在衡量$x$处两个分布之间的距离。同时，期望是以分布$p(x)$为基准计算的。因此，KL散度可以理解为以$p$为权重的、分布$p$与$q$的加权平均距离。从定义就可以看出，它不满足我们一般需要距离满足的对称性，即 $D_{\mathrm{KL}}(p\Vert q)\ne D_{\mathrm{KL}}(q\Vert p)$。但是，对任意两个分布，KL散度始终是非负的。

  我们以图7为例定性地说明这一点，图的左半部分展示了两个连续变量的概率密度函数$p$和$q$，右半部分的绿色曲线是 $p\log (p/q)$，其下方的面积就是KL散度的值。观察右边绿色曲线大于0和小于0的部分与左边$p$、$q$之间大小的对应关系可以发现，在 $p>q$ 的地方，$\log (p/q)$总是正数，其曲线下的面积也是正数；在 $p<q$ 的地方，虽然$\log (p/q)$是负数，但是其权重$p$较小，加权后的总面积总是小于$p>q$的部分。由于概率密度需要满足归一化性质，即其曲线下方的面积必须为1，不可能出现$q$始终在$p$上方的情况。因此，KL散度计算时，曲线下的正负面积抵消后总是非负，并且在 $p=q$ 时取到最小值0。

图7 KL散度示意图

  对于离散随机变量，KL散度可以进一步拆分为
$$\begin{array}{rl}{D}_{\mathrm{KL}}(p\Vert q)& ={\mathbb{E}}_{X\sim p(X)}\left[\log \frac{p(X)}{q(X)}\right]\\ & =\sum _{i=1}^n(p(X_i)\log p(X_i)-p(X_i)\log q(X_i))\\ & =-H(p)-\sum _{i=1}^{n}p(X_i)\log q(X_i)\\ & =-H(p)+H(p,q)\end{array}$$ 其中，$H(p,q)=-\sum _{i=1}^{n}p(X_i)\log q(X_i)$ 就称为分布$p$与$q$的交叉熵（cross entropy）。在二分类问题中，随机变量$X$对应样本$\mathbf{x}$的类别$y$，只有0和1两种取值。令$p(X)$等于样本$\mathbf{x}$的类别是$X$的概率；$q(X)$等于模型预测的样本类别为$X$的概率，即 $q(X=1)=f_{\mathbf{\theta }}(\mathbf{x}),q(X=0)=1-f_{\mathbf{\theta }}(\mathbf{x})$。我们期望模型预测的概率尽可能接近真实类别，因此要最小化$p$与$q$之间的距离$D_{\mathrm{KL}}(p\Vert q)$。而在离散化KL散度的定义中，$-H(p)$只与样本真实类别有关，无法通过模型优化。因此，我们只需要最小化交叉熵$H(p,q)$。将$p$与$q$用上述定义代入，可得
$$\begin{array}{rl}H(p,q)& =-\sum _{i=1}^{n}p(X_i)\log q(X_i)\\ & =-p(X=1)\log q(X=1)-p(X=0)\log q(X=0)\\ & =-y\log f_{\mathbf{\theta }}(\mathbf{x})-(1-y)\log (1-f_{\mathbf{\theta }}(\mathbf{x}))\end{array}$$ 如果再对所有样本的交叉熵求和，就得到总的交叉熵为$$H(p,q)=-\sum _{i=1}^N{y}_i\log f_{\mathbf{\theta }}({\mathbf{x}}_i)+(1-y_i)\log (1-f_{\mathbf{\theta }}({\mathbf{x}}_i))$$

  可以发现，总交叉熵恰好等于负的对数似然函数。因此，逻辑回归问题中，最大化对数似然函数与最小化交叉熵是等价的。事实上，无论是离散情况还是连续情况这一结论都成立，所以我们也经常称交叉熵为逻辑回归的损失函数。交叉熵在涉及到概率分布的模型中十分重要，我们在后面的树模型中还会再次用到这一概念。

七、拓展：广义线性模型

  本文介绍的逻辑斯谛回归模型的参数化假设为 $f_{\mathbf{\theta }}(\mathbf{x})=\sigma ({\mathbf{\theta }}^{\mathrm{T}}\mathbf{x})$，并不是传统意义上的线性模型，而是对线性映射${\mathbf{\theta }}^{\mathrm{T}}\mathbf{x}$又做了某种变换。该变换的目的是使线性的${\mathbf{\theta }}^{\mathrm{T}}\mathbf{x}$符合问题要求的性质。事实上，除了分类问题之外，还有许多线性模型无法直接应用的场景。例如在研究某一货物的销量和价格的关系时，如果用线性模型去拟合，可能得到“货物价格每上升1元，销量减少10个”这样的结论。然而，如果该货物原本的售价是10元、销量是100个，我们用该模型就会得到“货物售价30元时，销量是-100个”这样不现实的结论。这是由于货物销量$y$关于售价$x$并不完全是线性关系。于是，为了拓展线性模型的适用范围，我们在其基础上引入了广义线性模型（generalized linear model，GLM）的概念。前面已经讲过的线性回归模型和逻辑斯谛回归模型，都属于GLM的范畴。

  在逻辑斯谛回归中，我们假设$y$服从伯努利分布，且 $P(y=1)$ 关于$x$线性。事实上，普通的线性回归中也隐含了$y$服从正态分布的假设，并且其均值关于$x$线性，从而我们可以用最小二乘法去拟合回归直线。如果我们把$y$服从的分布拓展到指数分布族，就得到了GLM。指数分布族的定义如下：$$P(\mathbf{y}|\mathbf{\eta })=b(\mathbf{y})e^{{\mathbf{\eta }}^{\mathrm{T}}\mathbf{T}(\mathbf{y})-A(\mathbf{\eta })}$$ 其中，$\mathbf{\eta }=\mathbf{X}\mathbf{\theta }$ 称为线性预测子，标量函数$b(\mathbf{y})$、向量函数$\mathbf{T}(\mathbf{y})$和标量函数$A(\mathbf{\eta })$都是已知的。简单起见，下面我们采用其标量形式，并记 $\eta ={\mathbf{\theta }}^{\mathrm{T}}\mathbf{x}$：$$P(y|\eta )=b(y)e^{\eta T(y)-A(\eta )}$$ 通过选取不同的 $b$，$T$，$A$ 函数就可以得到不同的概率分布。例如，令$$b(y)=\frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{-\frac{y^2}{2{\sigma }^2}},T(y)=y,A(\eta )=\frac{{\eta }^2}{2}$$ 就得到高斯分布；令$$b(y)=1,T(y)=y,A(\eta )=\ln (1+e^{\eta })$$ 就得到伯努利分布。

在假设了$\mathbf{y}$服从的分布后，我们还希望计算$\mathbf{y}$的期望和方差，从而可以建立预测模型。指数分布族有非常好的性质，我们省略推导，直接给出其期望和方差分别为：$$\mathbb{E}[\mathbf{y}|\mathbf{\eta }]={\nabla }_{\mathbf{\theta }}A(\mathbf{\eta }),\mathrm{Var}[\mathbf{y}|\mathbf{\eta }]={\nabla }_{\mathbf{\theta }}^{2}A(\mathbf{\eta })$$ 接下来，我们建立参数化模型预测$\mathbf{y}$的期望，即 $f_{\mathbf{\theta }}(\mathbf{X})={\nabla }_{\mathbf{\theta }}A(\mathbf{\eta })$，再通过真实的训练数据来优化参数$\mathbf{\theta }$。这里，我们将高斯分布和伯努利分布代入，来还原线性回归和逻辑斯谛回归的模型。对于高斯分布，有$$f_{\mathbf{\theta }}(\mathbf{X})={\nabla }_{\mathbf{\theta }}A(\mathbf{\eta })={\nabla }_{\mathbf{\theta }}\left(\frac{{\mathbf{\eta }}^{\mathrm{T}}\mathbf{\eta }}{2}\right)=\mathbf{\eta }=\mathbf{X\theta }$$ 得到的结果与线性回归模型相同。对于伯努利分布，有$$f_{\mathbf{\theta }}(\mathbf{X})={\nabla }_{\mathbf{\theta }}A(\mathbf{\eta })={\nabla }_{\mathbf{\theta }}\ln (1+e^{\mathbf{\eta }})=\frac{e^{\mathbf{\eta }}}{1+e^{\mathbf{\eta }}}=\frac{1}{1+e^{-\mathbf{X\theta }}}=\sigma (\mathbf{X\theta })$$ 得到的结果与逻辑斯谛回归模型相同。可以发现，这里逻辑斯谛函数是由伯努利分布和GLM假设自然导出的。

  由线性模型到广义线性模型是一个由具象到抽象的过程。普通的线性模型的应用场景非常有限，但将其中的内核抽象出来，得到广义线性模型后，就可以解决所有适用于指数分布族的问题。除了我们详细介绍的线性回归模型和逻辑斯谛回归模型之外，泊松分布、多项分布、分类分布等都可以导出不同的模型，可以在不同的任务中发挥作用。

附：以上文中的数据集及相关资源下载地址：
链接：https://pan.quark.cn/s/ad0527701cee
提取码：VHL2