数据科学学习指南（六）

原文：zh.annas-archive.org/md5/9950f82165bee426119f99afc3ab612d

译者：飞龙

协议：CC BY-NC-SA 4.0

第十八章：案例研究：如何称量一只驴

驴子在肯尼亚农村扮演着重要角色。人们需要它们来运输作物、水和人，以及耕种田地。当一只驴子生病时，兽医需要知道它的体重以便开对量的药物。但是在肯尼亚农村，许多兽医没有称重器，因此他们需要猜测驴子的体重。药物用量过少可能会导致感染再次发作；药物过多则可能导致有害的过量。肯尼亚有超过 180 万头驴子，因此估算驴子体重的方法简单而准确非常重要。

在这个案例研究中，我们追随凯特·米尔纳和乔纳森·鲁吉尔 的工作，创建一个模型，供肯尼亚农村的兽医使用，以准确估计驴子的体重。像往常一样，我们遵循数据科学生命周期的步骤，但这次我们的工作偏离了这本书迄今为止涵盖的基础知识。您可以把这个案例研究看作是反思许多数据处理核心原则的机会，并理解如何扩展这些原则以应对具体情况。我们直接评估测量误差的来源，设计一个反映对过量用药关注的特殊损失函数，构建一个以适用性为重点的模型，并使用相对于驴子大小的特殊标准评估模型预测。

我们从数据的范围开始。

驴子研究问题与范围

我们的动力问题是：在没有称重器的情况下，兽医如何准确估算驴子的体重？让我们考虑一下他们更容易得到的信息。他们可以携带卷尺，并找出驴子在其他尺寸上的大小，比如它的高度。他们可以观察动物的性别，评估其一般状态，并询问驴子的年龄。因此，我们可以将我们的问题进一步细化为：兽医如何从易于获取的测量数据中准确预测驴子的体重？

为了解决这个更精确的问题，驴避难所 在肯尼亚农村的 17 个移动驱虫点进行了研究。

在范围方面（第二章），目标人群是肯尼亚农村的驴群。访问框架是所有被带到驱虫点的驴子。样本包括 2010 年 7 月 23 日至 8 月 11 日被带到这些点的所有驴子，但有一些注意事项：如果在一个地点有太多的驴子需要测量，科学家们会选择一部分来测量，并且任何怀孕或明显生病的驴子都被排除在研究之外。

为了避免不小心对同一只驴进行两次称重，每只驴在称重后都会标记。为了量化测量误差并评估称重过程的重复性，我们对 31 只驴进行了两次测量，而工作人员并不知道这是一只已经被重新称重的驴。

考虑到这一抽样过程，此数据的潜在偏差源包括：

覆盖偏差

这 17 个地点位于肯尼亚东部雅塔区和大裂谷内纳瓦沙区周围的地区。

选择偏差

只有被送往庇护所的驴才能参加这项研究，当一个地点有太多的驴时，会选择一个非随机样本。

测量偏差

除了测量误差外，称可能存在偏差。理想情况下，在使用场地之前和之后，天平应该校准（第十二章）。

尽管存在这些潜在的偏差源，但从肯尼亚农村地区拥有关心动物健康的主人那里获取驴的访问框架似乎是合理的。

我们的下一步是清理数据。

数据整理和转换

我们首先看一下数据文件的内容。为此，我们打开文件并检查前几行（第八章）：

`from` `pathlib` `import` `Path`

`# Create a Path pointing to our datafile`
`insp_path` `=` `Path``(``'``data/donkeys.csv``'``)`

`with` `insp_path``.``open``(``)` `as` `f``:`
    `# Display first five lines of file`
    `for` `_` `in` `range``(``5``)``:`
        `print``(``f``.``readline``(``)``,` `end``=``'``'``)`

BCS,Age,Sex,Length,Girth,Height,Weight,WeightAlt
3,<2,stallion,78,90,90,77,NA
2.5,<2,stallion,91,97,94,100,NA
1.5,<2,stallion,74,93,95,74,NA
3,<2,female,87,109,96,116,NA

由于文件是 CSV 格式的，我们可以轻松地将其读入数据框架：

`donkeys` `=` `pd``.``read_csv``(``"``data/donkeys.csv``"``)`
`donkeys`

	BCS	年龄	性别	长度	胸围	身高	体重	WeightAlt
0	3.0	<2	种马	78	90	90	77	NaN
1	2.5	<2	种马	91	97	94	100	NaN
2	1.5	<2	种马	74	93	95	74	NaN
...	...	...	...	...	...	...	...	...
541	2.5	10-15	种马	103	118	103	174	NaN
542	3.0	2-5	种马	91	112	100	139	NaN
543	3.0	5-10	种马	104	124	110	189	NaN

544 rows × 8 columns

超过五百只驴参与了调查，每只驴进行了八次测量。根据文档，粒度是单个驴（第九章）。表 18-1 提供了这八个特征的描述。

表 18-1. 驴研究代码手册

特征	数据类型	特征类型	描述
BCS	float64	有序	身体条件评分：从 1（消瘦）到 3（健康）到 5（肥胖），每 0.5 单位增加。
年龄	string	有序	年龄（年），小于 2 岁、2-5 岁、5-10 岁、10-15 岁、15-20 岁和 20 岁以上
性别	string	名义	性别类别：种马、阉马、母驴
长度	int64	数值	体长（厘米），从前腿肘到骨盆后部
胸围	int64	数值	身体围长（厘米），在前腿后面测量
身高	int64	数值	身体高度（厘米），到颈部连接背部的点
体重	int64	数值	体重（千克）
WeightAlt	float64	数值	在一部分驴身上进行的第二次称重测量

图 18-1 是将驴子理想化为带有颈部和附加的腿的圆柱体的图示。身高是从地面到肩膀上方的颈部底部的测量；胸围是围绕身体，就在后腿后面；长度是从前肘到骨盆后部。

图 18-1。驴的胸围、长度和高度的图示，被描述为对圆柱体上的测量

我们的下一步是对数据进行一些质量检查。在上一节中，我们基于范围列出了一些潜在的质量问题。接下来，我们检查测量和它们的分布的质量。

让我们首先比较对一小部分驴子进行的两次体重测量，以检查秤的一致性。我们为这 31 只被称重两次的驴子制作了这两个测量值之间的差异的直方图：

`donkeys` `=` `donkeys``.``assign``(``difference``=``donkeys``[``"``WeightAlt``"``]` `-` `donkeys``[``"``Weight``"``]``)`

`px``.``histogram``(``donkeys``,` `x``=``"``difference``"``,` `nbins``=``15``,`
    `labels``=``dict``(`
        `difference``=``"``Differences of two weighings (kg)<br>on the same donkey``"`
    `)``,`
    `width``=``350``,` `height``=``250``,`
`)`

这些测量值彼此之间都在 1 公斤以内，大部分都是完全相同的（四舍五入到最接近的公斤）。这让我们对测量的准确性有信心。

接下来，我们寻找体况评分中的异常值：

`donkeys``[``'``BCS``'``]``.``value_counts``(``)`

BCS
3.0    307
2.5    135
3.5     55
      ... 
1.5      5
4.5      1
1.0      1
Name: count, Length: 8, dtype: int64

从这个输出中，我们看到只有一个消瘦的（BCS = 1）和一个肥胖的（BCS = 4.5）驴子。让我们看看这两只驴子的完整记录：

`donkeys``[``(``donkeys``[``'``BCS``'``]` `==` `1.0``)` `|` `(``donkeys``[``'``BCS``'``]` `==` `4.5``)``]`

	BCS	年龄	性别	长度	胸围	身高	体重	体重备用
291	4.5	10-15	女性	107	130	106	227	NaN
445	1.0	>20	女性	97	109	102	115	NaN

由于这些 BCS 值非常极端，我们要谨慎考虑将这两只驴子纳入我们的分析范围。在这两个极端类别中，我们每个类别只有一只驴子，因此我们的模型可能无法延伸到 BCS 为 1 或 4.5 的驴子。我们将这两条记录从数据框中删除，并注意到我们的分析可能不适用于消瘦或肥胖的驴子。总的来说，我们在删除数据框中的记录时要小心。稍后，如果这五只 BCS 为 1.5 的驴子在我们的分析中看起来不正常，我们也可能决定将它们删除，但目前，我们将它们保留在我们的数据框中。一般来说，我们需要一个充分的理由来排除数据，并且我们应该记录这些操作，因为它们可能会影响我们的发现。如果我们删除与模型不符的任何记录，可能会导致过度拟合。

我们接下来删除这两个异常值：

`def` `remove_bcs_outliers``(``donkeys``)``:`
    `return` `donkeys``[``(``donkeys``[``'``BCS``'``]` `>``=` `1.5``)` `&` `(``donkeys``[``'``BCS``'``]` `<``=` `4``)``]` 

`donkeys` `=` `(``pd``.``read_csv``(``'``data/donkeys.csv``'``)`
           `.``pipe``(``remove_bcs_outliers``)``)`

现在，我们检查体重值的分布，以查看是否存在任何质量问题：

`px``.``histogram``(``donkeys``,` `x``=``'``Weight``'``,` `nbins``=``40``,` `width``=``350``,` `height``=``250``,`
            `labels``=``{``'``Weight``'``:``'``Weight (kg)``'``}``)`

看起来有一只非常轻的驴子，体重不到 30 公斤。接下来，我们检查体重和身高之间的关系，以评估用于分析的数据质量。

`px``.``scatter``(``donkeys``,` `x``=``'``Height``'``,` `y``=``'``Weight``'``,` `width``=``350``,` `height``=``250``,`
          `labels``=``{``'``Weight``'``:``'``Weight (kg)``'``,` `'``Height``'``:``'``Height (cm)``'``}``)`

小驴离主要驴群较远，会对我们的模型产生过大影响，因此我们将其排除。同样，我们要注意，如果有一两匹重驴对我们未来的模型拟合产生过大影响，我们也可能要将它们排除：

`def` `remove_weight_outliers``(``donkeys``)``:`
    `return` `donkeys``[``(``donkeys``[``'``Weight``'``]` `>``=` `40``)``]`

`donkeys` `=` `(``pd``.``read_csv``(``'``data/donkeys.csv``'``)`
           `.``pipe``(``remove_bcs_outliers``)`
           `.``pipe``(``remove_weight_outliers``)``)`

`donkeys``.``shape`

(541, 8)

总之，根据我们的清理和质量检查，我们从数据框中删除了三个异常观测。现在我们几乎可以开始我们的探索性分析了。在继续之前，我们将一些数据设置为测试集。

我们讨论了在第十六章中将测试集与训练集分开的重要性。最佳实践是在分析早期就将测试集分开，这样我们在详细探索数据之前就开始做出关于适合哪种模型以及在模型中使用哪些变量的决定。很重要的一点是，我们的测试集不参与这些决策，以便模拟我们的模型在全新数据上的表现。

我们将数据分成 80/20 的比例，其中 80%用于探索和建模。然后，我们用设置的 20%来评估模型。我们使用简单随机抽样将数据框分为测试集和训练集。首先，我们随机打乱数据框的索引：

`np``.``random``.``seed``(``42``)`
`n` `=` `len``(``donkeys``)`
`indices` `=` `np``.``arange``(``n``)`
`np``.``random``.``shuffle``(``indices``)`
`n_train` `=` `int``(``np``.``round``(``(``0.8` `*` `n``)``)``)`

接下来，我们将数据框的前 80%分配给训练集，剩余的 20%分配给测试集：

`train_set` `=` `donkeys``.``iloc``[``indices``[``:``n_train``]``]`
`test_set` `=` `donkeys``.``iloc``[``indices``[``n_train``:``]``]`

现在我们准备探索训练数据，寻找有用的关系和分布，为我们的建模提供信息。

探索中

我们查看数据框中的形状和关系特征，以便进行转换和模型制作（第十章）。我们首先看看年龄、性别和体况这些分类特征如何与体重相关联：

`f1` `=` `px``.``box``(``train_set``,` `x``=``"``Age``"``,` `y``=``"``Weight``"``,` 
            `category_orders` `=` `{``"``Age``"``:``[``'``<2``'``,` `'``2-5``'``,` `'``5-10``'``,` 
                                      `'``10-15``'``,` `'``15-20``'``,` `'``>20``'``]``}``)`
`f2` `=` `px``.``box``(``train_set``,` `x``=``"``Sex``"``,` `y``=``"``Weight``"``)`

`# We wrote the left_right function as a shorthand for plotly's make_subplots`
`fig` `=` `left_right``(``f1``,` `f2``,` `column_widths``=``[``0.7``,` `0.3``]``)`

`fig``.``update_xaxes``(``title``=``'``Age (yr)``'``,` `row``=``1``,` `col``=``1``)`
`fig``.``update_xaxes``(``title``=``'``Sex``'``,` `row``=``1``,` `col``=``2``)`
`fig``.``update_yaxes``(``title``=``'``Weight (kg)``'``,` `row``=``1``,` `col``=``1``)`

`fig` `=` `px``.``box``(``train_set``,` `x``=``"``BCS``"``,` `y``=``"``Weight``"``,` `points``=``"``all``"``,`
             `labels``=``{``'``Weight``'``:``'``Weight (kg)``'``,` `'``BCS``'``:``'``Body condition score``'``}``,`
             `width``=``550``,` `height``=``250``)`
`fig`

注意，我们绘制了身体状况评分的点和箱线图，因为我们之前看到评分为 1.5 的观测值只有少数几个，所以我们不希望过多解读仅有少量数据点的箱线图（第十一章）。看起来，体重中位数随着体况评分增加而增加，但增长并非简单线性。另一方面，三种性别类别的体重分布看起来大致相同。至于年龄，一旦一匹驴达到五岁，其体重分布似乎不会有太大变化。但两岁以下的驴和两至五岁的驴的体重普遍较低。

接下来，让我们检查定量变量。我们在散点图矩阵中绘制所有定量变量的成对关系：

骡子的身高、长度和腰围都与体重以及彼此之间线性相关。这并不太令人惊讶；只要知道骡子的一个维度，我们就可以大致猜测其他维度。腰围与体重的相关性最高，这在相关系数矩阵中得到了验证：

`train_numeric``.``corr``(``)`

	权重	长度	腰围	身高
体重	1.00	0.78	0.90	0.71
长度	0.78	1.00	0.66	0.58
腰围	0.90	0.66	1.00	0.70
身高	0.71	0.58	0.70	1.00

我们的探索揭示了数据的几个可能与建模相关的方面。我们发现骡子的腰围、长度和身高都与体重以及彼此之间线性相关，腰围与体重的线性关系最强。我们还观察到，体况评分与体重呈正相关；骡子的性别似乎与体重无关；对于 5 岁以上的骡子，年龄也与体重无关。在下一节中，我们将利用这些发现来构建我们的模型。

模拟骡子的体重

我们想建立一个简单的模型来预测骡子的体重。这个模型应该易于兽医在现场仅使用手算器时实现。模型也应该易于解释。

我们还希望模型能够依赖于兽医的情况——例如，他们是否正在开具抗生素或麻醉药。为了简洁起见，我们只考虑开具麻醉药的情况。我们第一步是选择一个反映这种情况的损失函数。

麻醉药开具的损失函数

麻醉药的过量可能比不足更糟。兽医很容易看出骡子麻醉药不足（它会抱怨），并且兽医可以给骡子再多点。但另一方面，麻醉药过多可能会有严重后果，甚至可能致命。因此，我们需要一个非对称的损失函数：对于体重的过高估计，它的损失应该大于对低估的损失。这与我们到目前为止在本书中使用的其他所有损失函数不同，它们都是对称的。

为了这个目的，我们创建了一个损失函数 anes_loss(x)：

`def` `anes_loss``(``x``)``:`
    `w` `=` `(``x` `>``=` `0``)` `+` `3` `*` `(``x` `<` `0``)`
    `return` `np``.``square``(``x``)` `*` `w`

相对误差为 100 ( y − y ^ ) / y ^ ，其中 y 为真实值，y ^ 为预测值。我们可以通过一个图示来展示损失函数的非对称性：

请注意，x 轴上值为-10 反映了 10%的过高估计。

接下来，让我们使用这个损失函数拟合一个简单 接下来，让我们使用这个损失函数拟合一个简单的线性模型。

拟合一个简单的线性模型

我们发现，腰围在我们的训练集中与体重的相关性最高。所以我们拟合了一个形式为的模型：

θ 0 + θ 1 Girth

要找到数据的最佳拟合θ 0和θ 1，我们首先创建一个具有腰围和截距项的设计矩阵。我们还创建了观察到的驴子体重的y向量：

`X` `=` `train_set``.``assign``(``intr``=``1``)``[``[``'``intr``'``,` `'``Girth``'``]``]`
`y` `=` `train_set``[``'``Weight``'``]`
`X`

	intr	Girth
230	1	116
74	1	117
354	1	123
...	...	...
157	1	123
41	1	103
381	1	106

433 rows × 2 columns

现在我们想要找到最小化数据上平均麻醉损失的θ 0和θ 1。为此，我们可以像在第十五章中一样使用微积分，但在这里我们将使用scipy包中的minimize方法，该方法执行数值优化（见第二十章）：

`from` `scipy``.``optimize` `import` `minimize`

`def` `training_loss``(``X``,` `y``)``:`
    `def` `loss``(``theta``)``:`
        `predicted` `=` `X` `@` `theta`
        `return` `np``.``mean``(``anes_loss``(``100` `*` `(``y` `-` `predicted``)` `/` `predicted``)``)`
    `return` `loss`

`results` `=` `minimize``(``training_loss``(``X``,` `y``)``,` `np``.``ones``(``2``)``)`
`theta_hat` `=` `results``[``'``x``'``]`

After fitting:
θ₀ = -218.51
θ₁ =    3.16

让我们看看这个简单模型的效果如何。我们可以使用模型来预测训练集上的驴子体重，然后找到预测中的误差。接下来的残差图显示了模型误差占预测值的百分比。相对于驴子的大小来说，预测误差较小更为重要，因为对于 100 公斤的驴子来说，10 公斤的误差比对于 200 公斤的驴子来说要糟糕得多。因此，我们计算每个预测的相对误差：

predicted = X @ theta_hat
resids = 100 * (y - predicted) / predicted

让我们来检查一下相对误差的散点图：

resid = pd.DataFrame({
    'Predicted weight (kg)': predicted, 'Percent error': resids})
px.scatter(resid, x='Predicted weight (kg)', y='Percent error',
           width=350, height=250)

使用最简单的模型，一些预测偏差达到了 20%至 30%。让我们看看稍微复杂一点的模型是否改善了预测。

拟合多元线性模型

让我们考虑进一步的模型，将其他数字变量纳入考虑。我们有三个数字变量来衡量驴子的腰围、长度和高度，有七种组合这些变量的模型：

[['Girth'],
 ['Length'],
 ['Height'],
 ['Girth', 'Length'],
 ['Girth', 'Height'],
 ['Length', 'Height'],
 ['Girth', 'Length', 'Height']]

对于这些变量组合中的每一个，我们可以使用我们的特殊损失函数来拟合一个模型。然后我们可以查看每个模型在训练集上的表现：

`def` `training_error``(``model``)``:`
    `X` `=` `train_set``.``assign``(``intr``=``1``)``[``[``'``intr``'``,` `*``model``]``]`
    `theta_hat` `=` `minimize``(``training_loss``(``X``,` `y``)``,` `np``.``ones``(``X``.``shape``[``1``]``)``)``[``'``x``'``]`
    `predicted` `=` `X` `@` `theta_hat`
    `return` `np``.``mean``(``anes_loss``(``100` `*` `(``y` `-` `predicted``)``/` `predicted``)``)`

`model_risks` `=` `[`
    `training_error``(``model``)`
    `for` `model` `in` `models`
`]`

	model	mean_training_error
0	[Girth]	94.36
1	[Length]	200.55
2	[Height]	268.88
3	[Girth, Length]	65.65
4	[Girth, Height]	86.18
5	[Length, Height]	151.15
6	[Girth, Length, Height]	63.44

正如我们之前所述，驴子的腰围是体重的最佳单一预测因子。然而，腰围和长度的组合的平均损失要比仅有腰围的损失要小得多，而且这个特定的双变量模型几乎和包含所有三个变量的模型一样好。由于我们想要一个简单的模型，我们选择了双变量模型而不是三变量模型。

接下来，我们使用特征工程将分类变量纳入模型中，这将改善我们的模型。

将定性特征纳入模型

在我们的探索性分析中，我们发现驴体条件和年龄的箱线图可能包含有助于预测体重的信息。由于这些是分类特征，我们可以将它们转换为 0-1 变量，采用独热编码，如第十五章所述。

独热编码允许我们调整模型中每个类别组合的截距项。我们当前的模型包括数值变量周长和长度：

θ 0 + θ 1 Girth + θ 2 Length

如果我们将年龄特征清理为三个类别——年龄<2、年龄 2-5和年龄>5——年龄的独热编码将创建三个 0-1 特征，每个类别一个。将独热编码特征包括在模型中得到：

θ 0 + θ 1 Girth  + θ 2 Length  + θ 3 Age<2  + θ 4 Age2-5 

在这个模型中，年龄<2表示小于 2 岁的驴为 1，否则为 0。类似地，年龄 2-5表示 2 至 5 岁的驴为 1，否则为 0。

我们可以将这个模型看作适合三个相同的线性模型，唯一不同的是常数的大小，因为该模型等同于：

( θ 0 + θ 3 ) + θ 1 Girth + θ 2 Length for a donkey under 2 ( θ 0 + θ 4 ) + θ 1 Girth + θ 2 Length for a donkey between 2 and 4 θ 0 + θ 1 Girth + θ 2 Length for a donkey over 5

现在让我们对我们的三个分类变量（体质、年龄和性别）都应用独热编码：

`X_one_hot` `=` `(`
    `train_set``.``assign``(``intr``=``1``)`
    `[``[``'``intr``'``,` `'``Length``'``,` `'``Girth``'``,` `'``BCS``'``,` `'``Age``'``,` `'``Sex``'``]``]`
    `.``pipe``(``pd``.``get_dummies``,` `columns``=``[``'``BCS``'``,` `'``Age``'``,` `'``Sex``'``]``)`
    `.``drop``(``columns``=``[``'``BCS_3.0``'``,` `'``Age_5-10``'``,` `'``Sex_female``'``]``)`
`)`
`X_one_hot`

	intr	长度	周长	BCS_1.5	...	年龄 _<2	年龄 _>20	性别 _ 阉割	性别 _ 种马
230	1	101	116	0	...	0	0	0	1
74	1	92	117	0	...	0	0	0	1
354	1	103	123	0	...	0	1	0	0
...	...	...	...	...	...	...	...	...	...
157	1	93	123	0	...	0	0	0	1
41	1	89	103	0	...	1	0	0	0
381	1	86	106	0	...	0	0	0	0

433 rows × 15 columns

对于每个分类特征，我们删除了一个虚拟变量。由于BCS、Age和Sex分别有六、六和三个类别，因此我们为设计矩阵添加了 12 个虚拟变量，总计 15 列，包括截距项、周长和长度。

让我们看看哪些分类变量，如果有的话，能够改进我们的双变量模型。为此，我们可以拟合包括所有三个分类特征的模型，以及周长和长度的模型：

`results` `=` `minimize``(``training_loss``(``X_one_hot``,` `y``)``,` `np``.``ones``(``X_one_hot``.``shape``[``1``]``)``)`

`theta_hat` `=` `results``[``'``x``'``]`

`y_pred` `=` `X_one_hot` `@` `theta_hat`
`training_error` `=` `(``np``.``mean``(``anes_loss``(``100` `*` `(``y` `-` `y_pred``)``/` `y_pred``)``)``)`

`print``(``f``'``Training error:` `{``training_error``:``.2f``}``'``)`

Training error: 51.47

根据平均损失，这个模型比只包括周长和长度的先前模型表现更好。但是让我们试着简化这个模型，同时保持其准确性。为此，我们查看每个虚拟变量的系数，看看它们是否接近 0，以及彼此之间的情况。换句话说，我们想看到如果将系数包括在模型中会如何改变截距。系数的图表可以轻松地进行比较：

系数证实了我们在箱线图中看到的情况。驴子性别的系数接近零，这意味着知道性别并不真正改变体重预测。我们还看到，将超过 5 岁的驴子的年龄类别合并起来将简化模型而不会损失太多。最后，由于体况评分为 1.5 的驴子数量很少，其系数接近 2 的体况评分，我们倾向于将这两个类别合并。

鉴于这些发现，我们更新设计矩阵：

`def` `combine_bcs``(``X``)``:`
    `new_bcs_2` `=` `X``[``'``BCS_2.0``'``]` `+` `X``[``'``BCS_1.5``'``]`
    `return` `X``.``assign``(``*``*``{``'``BCS_2.0``'``:` `new_bcs_2``}``)``.``drop``(``columns``=``[``'``BCS_1.5``'``]``)`

`def` `combine_age_and_sex``(``X``)``:`
    `return` `X``.``drop``(``columns``=``[``'``Age_10-15``'``,` `'``Age_15-20``'``,` `'``Age_>20``'``,`
                           `'``Sex_gelding``'``,` `'``Sex_stallion``'``]``)`

`X_one_hot_simple` `=` `(`
    `X_one_hot``.``pipe``(``combine_bcs``)`
    `.``pipe``(``combine_age_and_sex``)`
`)`

然后我们拟合更简单的模型：

`results` `=` `minimize``(``training_loss``(``X_one_hot_simple``,` `y``)``,`
                   `np``.``ones``(``X_one_hot_simple``.``shape``[``1``]``)``)`
`theta_hat` `=` `results``[``'``x``'``]`
`y_pred` `=` `X_one_hot_simple` `@` `theta_hat`
`training_error` `=` `(``np``.``mean``(``anes_loss``(``100` `*` `(``y` `-` `y_pred``)``/` `y_pred``)``)``)`
`print``(``f``'``Training error:` `{``training_error``:``.2f``}``'``)`

Training error: 53.20

平均误差与更复杂的模型接近，因此我们决定采用这个更简单的模型。让我们显示系数并总结模型：

	var	theta_hat
0	intr	-175.25
1	长度	1.01
2	腹围	1.97
3	BCS_2.0	-6.33
4	BCS_2.5	-5.11
5	BCS_3.5	7.36
6	BCS_4.0	20.05
7	年龄 _2-5	-3.47
8	年龄 _<2	-6.49

我们的模型大致是：

Weight ≈ − 175 + Length + 2 Girth

在这个初步的近似之后，我们使用分类特征进行一些调整：

• BCS 2 或更低？减去 6.5 公斤。

• BCS 2.5？减去 5.1 公斤。

• BCS 3.5？增加 7.4 公斤。

• BCS 4？增加 20 公斤。

• 年龄小于 2 岁？减去 6.5 公斤。

• 年龄在 2 到 5 岁之间？减去 3.5 公斤。

这个模型似乎相当简单易行，因为在我们根据驴子的长度和腹围进行初始估计之后，我们根据几个是/否问题的答案添加或减少一些数字。让我们看看这个模型在预测测试集中驴子的体重方面表现如何。

模型评估

记住，在探索和建模剩余的 80%数据之前，我们将 20%的数据搁置。现在，我们已经准备好将我们从训练集中学到的应用到测试集中。也就是说，我们拿出我们拟合的模型，并用它来预测测试集中驴子的体重。为此，我们需要准备测试集。我们的模型使用驴子的腹围和长度，以及驴子年龄和体况评分的虚拟变量。我们将我们在训练集上的所有转换应用到我们的测试集上：

`y_test` `=` `test_set``[``'``Weight``'``]`

`X_test` `=` `(`
    `test_set``.``assign``(``intr``=``1``)`
    `[``[``'``intr``'``,` `'``Length``'``,` `'``Girth``'``,` `'``BCS``'``,` `'``Age``'``,` `'``Sex``'``]``]`
    `.``pipe``(``pd``.``get_dummies``,` `columns``=``[``'``BCS``'``,` `'``Age``'``,` `'``Sex``'``]``)`
    `.``drop``(``columns``=``[``'``BCS_3.0``'``,` `'``Age_5-10``'``,` `'``Sex_female``'``]``)`
    `.``pipe``(``combine_bcs``)`
    `.``pipe``(``combine_age_and_sex``)`
`)`

我们将我们对设计矩阵的所有操作汇总到我们在训练集建模中确定的最终版本中。现在我们准备使用我们在训练集上拟合的θ来为测试集中的那些驴子进行体重预测：

`y_pred_test` `=` `X_test` `@` `theta_hat`
`test_set_error` `=` `100` `*` `(``y_test` `-` `y_pred_test``)` `/` `y_pred_test`

然后我们可以绘制相对预测误差：

记住，正的相对误差意味着低估重量，这并不像高估重量那样严重。从这个残差图中，我们看到几乎所有的测试集重量都在预测的 10%之内，只有一个超过 10%的误差是在高估的一侧。鉴于我们的损失函数更严厉地惩罚了高估，这是合理的。

另一种散点图显示了实际值和预测值，同时标出了 10%误差线，提供了不同的视角：

对于较大重量的预测线，10%的误差线距离预测线更远。

我们已经实现了我们的目标！我们有一个使用易于获取的测量数据的模型，简单到可以在说明书上解释，并且可以在实际驴重量的预测中保持在 10%范围内。接下来，我们总结这个案例研究并反思我们的模型。

总结

在这个案例研究中，我们展示了建模的不同目的：描述、推理和预测。对于描述，我们寻求一个简单易懂的模型。我们手工制作了这个模型，从分析探索阶段的发现开始。我们每采取一项行动来包含一个特征在模型中，折叠类别或转换特征，都是我们在调查数据时做出的决策。

在对象模拟自然现象如驴的重量时，我们理想地会使用物理模型和统计模型。在这种情况下，物理模型是用圆柱体表示驴的形象。一个好奇的读者可能会指出，我们可以直接使用这个表示来估算驴（圆柱体）的重量，通过其长度和围长（因为围长是 2 π r ）：

w e i g h t ∝ g i r t h 2 × l e n g t h

这个物理模型表明，对数转换的重量大致上是围长和长度的线性函数：

log ⁡ ( w e i g h t ) ∝ 2 log ⁡ ( g i r t h ) + log ⁡ ( l e n g t h )

鉴于这个物理模型，您可能会想知道为什么我们没有在拟合模型时使用对数或平方变换。我们留给您去更详细地探讨这样的模型。但总的来说，如果测量值的范围较小，那么对数函数大致上是线性的。为了保持我们的模型简单，鉴于围长和重量之间的高相关性，我们选择不进行这些转换。

在这个建模过程中，我们进行了大量的数据挖掘。我们检查了所有可能的模型，这些模型由数值特征的线性组合构建，并检查了虚拟变量的系数，以决定是否合并类别。当我们使用这种迭代方法创建模型时，非常重要的是我们留出数据来评估模型。在新数据上评估模型可以让我们放心地知道我们选择的模型效果如何。我们留出的数据在建模时没有参与任何决策，因此它能很好地帮助我们了解模型在预测上的表现。

我们应该记住之前描述的数据范围及其潜在偏见。我们的模型在测试集上表现良好，但测试集和训练集来自同一数据收集过程。只要新数据的范围保持不变，我们期望我们的模型在实践中表现良好。

最后，这个案例研究展示了模型拟合通常是在简单与复杂之间、物理与统计模型之间取得平衡的过程。物理模型可以作为建模的良好起点，而统计模型则可以为物理模型提供信息。作为数据科学家，我们需要在分析的每一步都做出判断。建模既是艺术又是科学。

本案例研究及其前面几章集中讨论了拟合线性模型的问题。接下来，我们考虑一种不同类型的建模方法，用于解释或预测的响应变量是定性而非定量的情况。

第六部分：分类

第十九章：分类

本章继续探讨数据科学生命周期的第四阶段：拟合和评估模型以理解世界。到目前为止，我们已经描述了如何使用绝对误差拟合常数模型（第四章）以及使用平方误差拟合简单和多元线性模型（第十五章）。我们还拟合了带有不对称损失函数的线性模型（第十八章）和带有正则化损失的线性模型（第十六章）。在所有这些情况下，我们的目标是预测或解释数值结果的行为——公交等待时间、空气中的烟粒子和驴子的体重都是数值变量。

在本章中，我们扩展了对建模的视角。我们不再预测数值结果，而是构建模型来预测名义结果。这些模型使银行能够预测信用卡交易是否欺诈，医生能够将肿瘤分类为良性或恶性，以及您的电子邮件服务能够识别垃圾邮件并将其与常规邮件区分开来。这种类型的建模称为分类，在数据科学中广泛应用。

就像线性回归一样，我们制定一个模型，选择一个损失函数，通过最小化数据的平均损失来拟合模型，并评估拟合的模型。但与线性回归不同的是，我们的模型不是线性的，损失函数也不是平方误差，我们的评估比较不同类型的分类错误。尽管存在这些差异，模型拟合的整体结构在这种情况下仍然适用。回归和分类共同组成了监督学习的主要方法，即基于观察结果和协变量拟合模型的一般任务。

我们首先介绍一个示例，在本章中我们将一直使用它。

例子：受风损坏的树木

1999 年，一场风速超过 90 英里每小时的巨大风暴损坏了边界水道独木舟区野外（BWCAW）中数百万棵树木，该地区是美国东部最大的原始森林地带。为了了解树木对风害的敏感性，名叫罗伊·劳伦斯·里奇的研究人员对 BWCAW 进行了地面调查。在此研究后的几年里，其他研究人员利用这一数据集对风倒（即树木在强风中被连根拔起）进行了建模。

研究对象是 BWCAW 中的树木群落。访问框架是样线：穿过自然景观的直线。这些特定的样线从湖边开始，沿着地形的梯度直角行驶 250 至 400 米。沿着这些样线，调查员每隔 25 米停下来，检查一个 5 乘 5 米的小区。在每个小区，树木被计数，分类为倒伏或直立，以 6 英尺高处的直径测量，并记录它们的种类。

像这样的采样协议在研究自然资源时很常见。在 BWCAW 中，该地区 80%以上的土地距离湖泊不到 500 米，因此几乎覆盖了整个人口。该研究于 2000 年和 2001 年夏季进行，1999 年暴风雨和数据收集之间没有发生其他自然灾害。

收集了 3600 多棵树的测量数据，但在这个例子中，我们仅研究了黑云杉。有 650 多棵。我们读取了这些数据：

`trees` `=` `pd``.``read_csv``(``'``data/black_spruce.csv``'``)`
`trees`

	直径	暴风雨	状态
0	9.0	0.02	站立
1	11.0	0.03	站立
2	9.0	0.03	站立
...	...	...	...
656	9.0	0.94	倒下
657	17.0	0.94	倒下
658	8.0	0.98	倒下

659 rows × 3 columns

每行对应一个单独的树，具有以下属性：

直径

直径为厘米，测量高度在地面以上 6 英尺处的树木

暴风雨

暴风雨的严重程度（25 米宽区域内倒下的树木占比）

状态

树木“倒下”或“站立”

在我们转向建模之前，让我们进行一些探索性分析。首先，我们计算一些简单的摘要统计信息：

`trees``.``describe``(``)``[``3``:``]`

	直径	暴风雨
最小	5.0	0.02
25%	6.0	0.21
50%	8.0	0.36
75%	12.0	0.55
最大	32.0	0.98

基于四分位数，树直径的分布似乎向右倾斜。让我们用直方图比较站立和倒下树木的直径分布：

在暴风雨中倒下的树木直径分布以 12 厘米为中心，呈右偏态。相比之下，站立的树几乎都在 10 厘米以下，众数约为 6 厘米（研究中仅包括直径至少为 5 厘米的树木）。

还有一个要调查的特征是暴风雨的强度。我们将暴风雨强度与树木直径绘制成图，使用符号和标记颜色来区分站立的树木和倒下的树木。由于直径基本上是以厘米为单位测量的，所以许多树木具有相同的直径，因此我们通过为直径值添加一些噪声来减少过度绘制（参见第十一章）。我们还调整了标记颜色的不透明度，以显示图中的密集区域：

从这张图中可以看出，树木直径和暴风雨的强度与风倒有关：树木是被连根拔起还是留在原地。请注意，我们想要预测的风倒是一个名义变量。在下一节中，我们考虑了这如何影响预测问题。

建模和分类

我们希望创建一个解释树木易受风倒的模型。换句话说，我们需要为两级名义特征构建模型：倒下或站立。当响应变量是名义的时，这个建模任务称为分类。在这种情况下只有两个级别，所以这个任务更具体地称为二元分类。

一个常数模型

让我们首先考虑最简单的模型：一个始终预测一类的常数模型。我们使用C来表示常数模型的预测。对于我们的风倒数据集，这个模型将为每个输入预测C = 站立或C = 倒下。

在分类中，我们想追踪我们的模型多频繁地预测了正确的类别。现在，我们只是使用正确预测的计数。这有时被称为零一误差，因为损失函数有两个可能的值之一：当进行错误的预测时为 1，进行正确的预测时为 0。对于给定的观察结果y i和预测C，我们可以将这个损失函数表示为：

ℓ ( C , y ) = { 0 when  C  matches  y 1 when  C  is a mismatch for  y

当我们收集了数据y = [ y 1 , … , y n ]时，平均损失为：

L ( C , y ) = 1 n ∑ i = 1 n ℓ ( C , y ) = #  mismatches n

对于常数模型（见第四章），当C设置为最常见的类别时，模型可以最小化损失。

就黑云杉而言，我们有以下比例的站立和倒下的树木：

`trees``[``'``status``'``]``.``value_counts``(``)` `/` `len``(``trees``)`

status
standing    0.65
fallen      0.35
Name: count, dtype: float64

所以我们的预测是一棵树站立，而我们数据集的平均损失为0.35。

这就是说，这个预测并不特别有用或有见地。例如，在我们对树木数据集进行的 EDA 中，我们发现树木的大小与树木是否站立或倒下有关联。理想情况下，我们可以将这些信息纳入模型，但常数模型不允许我们这样做。让我们对如何将预测因子纳入我们的模型建立一些直觉。

检查尺寸和风倒之间的关系

我们想更仔细地研究树木尺寸与风倒的关系。为了方便起见，我们将名义风倒特征转换为 0-1 数字特征，其中 1 表示倒下的树木，0 表示站立的树木：

`trees``[``'``status_0_1``'``]` `=` `(``trees``[``'``status``'``]` `==` `'``fallen``'``)``.``astype``(``int``)`
`trees`

	直径	风暴	状态	状态 _0_1
0	9.0	0.02	站立	0
1	11.0	0.03	站立	0
2	9.0	0.03	站立	0
...	...	...	...	...
656	9.0	0.94	倒下	1
657	17.0	0.94	倒下	1
658	8.0	0.98	倒下	1

659 rows × 4 columns

这种表示在许多方面都是有用的。例如，status_0_1 的平均值是数据集中倒下的树木比例：

`pr_fallen` `=` `np``.``mean``(``trees``[``'``status_0_1``'``]``)`
`print``(``f``"``Proportion of fallen black spruce:` `{``pr_fallen``:``0.2f``}``"``)`

Proportion of fallen black spruce: 0.35

具备这种 0-1 特征，我们可以制作一张图来展示树木直径与风倒之间的关系。这类似于我们进行线性回归的过程，其中我们绘制结果变量与解释变量之间的散点图（参见 第十五章）。

在这里，我们将树木状态绘制为直径的函数，但是我们在状态中添加了一点随机噪声，以帮助我们查看每个直径处 0 和 1 值的密度。与之前一样，我们也会扰动直径值，并调整标记的不透明度以减少重叠绘制。我们还在倒下的树木比例处添加了一条水平线：

这个散点图显示较小的树更有可能直立，而较大的树更有可能倒下。请注意，树木的平均状态（0.35）基本上适合将一个恒定模型应用于响应变量。如果我们将树木直径视为一个解释特征，我们应该能够改进模型。

一个起点可能是计算不同直径树木的倒下比例。以下代码块将树木直径分成区间，并计算每个区间内倒下的树木比例：

`splits` `=` `[``4``,` `5``,` `6``,` `7``,` `8``,` `9``,` `10``,` `12``,` `14``,` `17``,` `20``,` `25``,` `32``]`
`tree_bins` `=` `(`
    `trees``[``"``status_0_1``"``]`
    `.``groupby``(``pd``.``cut``(``trees``[``"``diameter``"``]``,` `splits``)``)`
    `.``agg``(``[``"``mean``"``,` `"``count``"``]``)`
    `.``rename``(``columns``=``{``"``mean``"``:` `"``proportion``"``}``)`
    `.``assign``(``diameter``=``lambda` `df``:` `[``i``.``right` `for` `i` `in` `df``.``index``]``)`
`)`

我们可以将这些比例绘制成树木直径的函数图：

标记的大小反映了直径区间内树木的数量。我们可以利用这些比例来改进我们的模型。例如，对于直径为 6 厘米的树木，我们会将其分类为直立，而对于 20 厘米的树木，我们的分类则是倒下的。二元分类的一个自然起点是对观察到的比例进行建模，然后利用这些比例进行分类。接下来，我们为这些比例开发一个模型。

建模比例（和概率）

回顾一下，当我们建模时，我们需要选择三样东西：一个模型，一个损失函数，以及一种方法来最小化训练集上的平均损失。在前一节中，我们选择了一个恒定模型，0-1 损失，并进行了一些适合模型的证明。然而，这个恒定模型并没有包含预测变量。在本节中，我们通过引入一个称为逻辑模型的新模型来解决这个问题。

为了推动这些模型，注意到树木直径与倒下树木比例之间的关系似乎不是线性的。为了演示，让我们对这些数据拟合一个简单的线性模型，以显示它具有几个不良特征。使用 第十五章 中的技术，我们对树木状态与直径进行了线性模型的拟合：

`from` `sklearn``.``linear_model` `import` `LinearRegression`
`X` `=` `trees``[``[``'``diameter``'``]``]`
`y` `=` `trees``[``'``status_0_1``'``]`

`lin_reg` `=` `LinearRegression``(``)``.``fit``(``X``,` `y``)`

然后，我们将这条拟合线添加到比例散点图中：

显然，模型对比例的拟合效果并不理想。存在几个问题：

• 模型对大树给出大于 1 的比例。

• 模型没有捕捉到比例中的曲线特征。

• 极端点（例如直径为 30 厘米的树木）将拟合线向右移动，远离大部分数据。

为了解决这些问题，我们引入了逻辑模型。

逻辑模型

逻辑模型是最广泛使用的基础分类模型之一，是线性模型的简单扩展。逻辑函数，通常称为sigmoid 函数，定义如下：

logistic ( t ) = 1 1 + exp ⁡ ( − t )

警告

Sigmoid函数通常用σ ( t )表示。不幸的是，希腊字母σ在数据科学和统计学中有多种含义，如标准差、逻辑函数和置换。看到σ时，必须根据上下文理解其含义。

我们可以绘制逻辑函数以显示其 S 形状，并确认其输出在 0 到 1 之间。函数随着t单调增加，t的大值接近 1：

`def` `logistic``(``t``)``:`
    `return` `1.` `/` `(``1.` `+` `np``.``exp``(``-``t``)``)`

由于逻辑函数映射到 0 到 1 之间的区间，通常用于建模比例和概率。此外，我们可以将逻辑写成线性函数的形式，如θ 0 + θ 1 x：

σ ( θ 0 + θ 1 x ) = 1 1 + exp ⁡ ( − θ 0 − θ 1 x )

为了帮助你直观地理解该函数的形状，下图显示了我们变化θ 0和θ 1时的逻辑函数。

我们可以看到改变θ 1的大小会改变曲线的陡峭程度；离 0 越远，曲线越陡。改变θ 1的符号将曲线反映在竖直线x = 0周围。改变θ 0会使曲线左右移动。

逻辑函数可以看作是一种转换：将线性函数转换为非线性平滑曲线，其输出始终位于 0 到 1 之间。实际上，逻辑函数的输出具有更深层次的概率解释，接下来我们将描述它。

对数几率

记住，赔率是概率p的比率p / ( 1 − p )。例如，投掷一个公平硬币时，得到正面的赔率是 1；对于一个比尾部更有可能出现两倍的硬币（p = 2 / 3），得到正面的赔率是 2。逻辑模型也称为对数几率模型，因为逻辑函数与对数几率的线性函数重合。

我们可以在以下方程中看到这一点。为了展示这一点，我们将 Sigmoid 函数的分子和分母分别乘以 exp ⁡ ( t )：

σ ( t ) = 1 1 + exp ⁡ ( − t ) = exp ⁡ ( t ) 1 + exp ⁡ ( t ) ( 1 − σ ( t ) ) = 1 − exp ⁡ ( t ) 1 + exp ⁡ ( t ) = 1 1 + exp ⁡ ( t )

然后我们取对数几率并简化：

log ⁡ ( σ ( t ) 1 − σ ( t ) ) = log ⁡ ( exp ⁡ ( t ) ) = t

因此，对于 σ ( θ 0 + θ 1 x ) ，我们发现对数几率是 x 的线性函数：

log ⁡ ( σ ( θ 0 + θ 1 x ) 1 − σ ( θ 0 + θ 1 x ) ) = log ⁡ ( exp ⁡ ( θ 0 + θ 1 x ) ) = θ 0 + θ 1 x

以对数几率的形式表示 logistic 对于系数 θ 1 给出了一个有用的解释。假设解释变量增加 1，则几率的变化如下：

 odds  = exp ⁡ ( θ 0 + θ 1 ( x + 1 ) ) = exp ⁡ ( θ 1 ) × exp ⁡ ( θ 0 + θ 1 x )

我们看到几率增加或减少了 exp ⁡ ( θ 1 ) 倍。

注意

这里，log 函数是自然对数。由于自然对数在数据科学中是默认的，因此我们通常不必写成 ln。

接下来，让我们在我们的比例图中添加一个 logistic 曲线，以了解它对数据的拟合效果。

使用 Logistic 曲线

在下面的图中，我们在倒下的树木比例的图上添加了一个 logistic 曲线：

我们可以看到曲线相对比例很好。事实上，我们选择了这个特定的 logistic 通过将其拟合到数据。拟合的 logistic 回归是：

σ(-7.4 + 3.0x)

现在我们已经看到 logistic 曲线可以很好地建模概率，我们转向将 logistic 曲线拟合到数据的过程。在下一节中，我们继续我们建模的第二步：选择一个合适的损失函数。

Logistic 模型的损失函数

Logistic 模型给我们提供了概率（或经验比例），因此我们将损失函数写成 ℓ ( p , y ) ，其中 p 在 0 和 1 之间。响应采用两个值之一，因为我们的输出特征是二元分类。因此，任何损失函数都可以简化为：

ℓ ( p , y ) = { ℓ ( p , 0 ) if  y  is 0 ℓ ( p , 1 ) if  y  is 1

再次使用 0 和 1 来表示类别具有优势，因为我们可以方便地写出损失函数为：

ℓ ( p , y ) = y ℓ ( p , y ) + ( 1 − y ) ℓ ( p , 1 − y )

我们鼓励你通过考虑 y = 1 和 y = 0 两种情况来确认这种等价性。

Logistic 模型与对数损失配合得很好：

ℓ ( p , y ) = { − log ⁡ ( p ) if  y  is 1 − log ⁡ ( 1 − p ) if  y  is 0 = − y log ⁡ ( p ) − ( 1 − y ) log ⁡ ( 1 − p )

注意 log 损失在 0 和 1 处未定义，因为− log ⁡ ( p )当p趋向 0 时趋近于∞，类似地，− log ⁡ ( 1 − p )当p趋向 1 时也如此。我们在最小化过程中需要小心避免这些端点。下图显示了这两种损失函数的情况：

当y为 1（实线）时，在p接近 1 时损失较小，当y为 0（虚线）时，在p接近 0 时损失较小。

如果我们的目标是使用 log 损失对数据拟合一个常数，那么平均损失为：

L ( p , y ) = 1 n ∑ i [ − y i log ⁡ ( p ) − ( 1 − y i ) log ⁡ ( 1 − p ) ] = − n 1 n log ⁡ ( p ) − n 0 n log ⁡ ( 1 − p ) )

这里n 0和n 1分别是值为 0 和 1 的y i的数量。我们可以对p进行微分以找到最小值点。

∂ L ( p , y ) ∂ p = − n 1 n p + n 0 n ( 1 − p )

然后我们将导数设置为 0 并解出最小化值p ^：

0 = − n 1 n p ^ + n 0 n ( 1 − p ^ ) 0 = − p ^ ( 1 − p ^ ) n 1 p ^ + p ^ ( 1 − p ^ ) n 0 ( 1 − p ^ ) n 1 ( 1 − p ^ ) = n 0 p ^ p ^ = n 1 n

（最终的方程来自于注意到n 0 + n 1 = n。）

要基于逻辑函数拟合更复杂的模型，我们可以将σ ( θ 0 + θ 1 x )代入p。逻辑模型的损失变为：

ℓ ( σ ( θ 0 + θ 1 x ) , y ) = y ℓ ( σ ( θ 0 + θ 1 x ) , y ) + ( 1 − y ) ℓ ( σ ( θ 0 + θ 1 x ) , 1 − y ) = y log ⁡ ( σ ( θ 0 + θ 1 x ) ) + ( 1 − y ) log ⁡ ( σ ( θ 0 + θ 1 x ) )

对数据的损失进行平均，我们得到：

L ( θ 0 , θ 1 , x , y ) = 1 n ∑ i − y i log ⁡ ( σ ( θ 0 + θ 1 x i ) ) − ( 1 − y i ) log ⁡ ( 1 − σ ( θ 0 + θ 1 x i ) )

与平方损失不同，此损失函数没有封闭形式的解。我们使用像梯度下降这样的迭代方法（见第二十章）来最小化平均损失。这也是我们不在逻辑模型中使用平方误差损失的原因之一——平均平方误差是非凸的，这使得优化变得困难。凸性的概念在第二十章有更详细的讨论，图 20-4 提供了直观的图示。

注意

log 损失也称为逻辑损失和交叉熵损失。它的另一个名称是负对数似然。这个名称指的是使用似然性来拟合模型的技术，即我们的数据来自于某个概率分布的似然性。在这里我们不深入探讨这些替代方法的背景。

拟合逻辑模型（使用 log 损失）被称为逻辑回归。逻辑回归是广义线性模型的一个示例，它是带有非线性变换的线性模型。

我们可以使用 scikit-learn 来拟合逻辑模型。包的设计者使 API 与最小二乘法拟合线性模型非常相似（见第十五章）。首先，我们导入逻辑回归模块：

`from` `sklearn``.``linear_model` `import` `LogisticRegression`

然后，我们用结果 y，即树的状态，和协变量 X，即直径（我们已对其进行了对数变换），设置回归问题：

`trees``[``'``log_diam``'``]` `=` `np``.``log``(``trees``[``'``diameter``'``]``)` 
`X` `=` `trees``[``[``'``log_diam``'``]``]`
`y` `=` `trees``[``'``status_0_1``'``]`

然后，我们拟合逻辑回归，并检查直径的截距和系数：

`lr_model` `=` `LogisticRegression``(``)`
`lr_model``.``fit``(``X``,` `y``)`

`[``intercept``]` `=` `lr_model``.``intercept_`
`[``[``coef``]``]` `=` `lr_model``.``coef_`
`print``(``f``'``Intercept:` `{``intercept``:``.1f``}``'``)`
`print``(``f``'``Diameter coefficient:` `{``coef``:``.1f``}``'``)`

Intercept:           -7.4
Diameter coefficient: 3.0

在进行预测时，predict 函数返回预测的（最可能的）类别，而 predict_proba 返回预测的概率。对于直径为 6 的树，我们预计预测为 0（即 站立 ）的概率很高。我们来检查一下：

`diameter6` `=` `pd``.``DataFrame``(``{``'``log_diam``'``:` `[``np``.``log``(``6``)``]``}``)`
`[``pred_prof``]` `=` `lr_model``.``predict_proba``(``diameter6``)`
`print``(``f``'``Predicted probabilities:` `{``pred_prof``}``'``)`

Predicted probabilities: [0.87 0.13]

因此，模型预测直径为 6 的树以 站立 类别有 0.87 的概率，以 倒下 类别有 0.13 的概率。

现在我们已经用一个特征拟合了一个模型，我们可能想要看看是否包含另一个特征，比如风暴的强度，是否可以改善模型。为此，我们可以通过将一个特征添加到 X 中，并再次拟合模型来拟合多元逻辑回归。

请注意，逻辑回归拟合一个模型来预测概率——模型预测直径为 6 的树以 站立 类别有 0.87 的概率，以 倒下 类别有 0.13 的概率。由于概率可以是介于 0 和 1 之间的任何数，我们需要将概率转换回类别以执行分类。我们将在下一节中解决这个分类问题。

从概率到分类

我们在本章开始时介绍了一个二元分类问题，我们想要建模一个名义响应变量。到目前为止，我们已经使用逻辑回归来建模比例或概率，现在我们准备返回原始问题：我们使用预测的概率来分类记录。对于我们的例子，这意味着对于特定直径的树，我们使用逻辑回归中拟合的系数来估计其倒下的可能性。如果可能性很高，我们将树分类为倒下；否则，我们将其分类为站立。但是我们需要选择一个阈值来制定这个 决策规则。

sklearn 的逻辑回归模型的 predict 函数实现了基本的决策规则：如果预测的概率 p > 0.5 ，则预测 1 。否则，预测 0 。我们将这个决策规则以虚线叠加在模型预测之上：

在本节中，我们考虑一个更一般的决策规则。对于某些选择的 τ ，如果模型预测的概率 p > τ ，则预测 1 ，否则预测 0 。默认情况下，sklearn 设置 τ = 0.5 。让我们探讨当 τ 被设置为其他值时会发生什么。

选择适当的τ值取决于我们的目标。假设我们希望最大化准确率。分类器的准确率是正确预测的分数。我们可以计算不同阈值下的准确率，即不同的τ值：

`def` `threshold_predict``(``model``,` `X``,` `threshold``)``:`
    `return` `np``.``where``(``model``.``predict_proba``(``X``)``[``:``,` `1``]` `>` `threshold``,` `1.0``,` `0.0``)`

`def` `accuracy``(``threshold``,` `X``,` `y``)``:`
    `return` `np``.``mean``(``threshold_predict``(``lr_model``,` `X``,` `threshold``)` `==` `y``)`

`thresholds` `=` `np``.``linspace``(``0``,` `1``,` `200``)`
`accs` `=` `[``accuracy``(``t``,` `X``,` `y``)` `for` `t` `in` `thresholds``]`

要理解准确率如何随τ变化而变化，我们制作了一个图表：

注意，具有最高准确率的阈值并不完全在 0.5 处。在实践中，我们应该使用交叉验证来选择阈值（参见第十六章）。

最大化准确率的阈值可能不是 0.5，这有许多原因，但一个常见的原因是类别不平衡，其中一个类别比另一个类别频繁。类别不平衡可能导致模型将记录分类为更常见的类别。在极端情况下（如欺诈检测），当数据中只有很小一部分包含特定类别时，我们的模型可以通过始终预测频繁类别而不学习如何生成适合稀有类别的好分类器来实现高准确率。有一些管理类别不平衡的技术，例如：

• 对数据进行重新采样以减少或消除类别不平衡

• 调整损失函数以对较小类别施加更大的惩罚

在我们的示例中，类别不平衡并不那么严重，因此我们继续进行而不进行这些调整。

类别不平衡问题解释了为什么单靠准确率通常不是我们评判模型的方式。相反，我们希望区分不同类型的正确和错误分类。我们接下来描述这些内容。

混淆矩阵

在二元分类中可视化错误的一个方便方法是查看混淆矩阵。混淆矩阵比较模型预测与实际结果。在这种情况下存在两种类型的错误：

假阳性

当实际类别为 0（错误）但模型预测为 1（真实）

假阴性

当实际类别为 1（真实）但模型预测为 0（错误）

理想情况下，我们希望尽量减少两种错误，但我们经常需要平衡这两种来源之间的关系。

注意

“正面”和“负面”这些术语来自于疾病测试，其中指示存在疾病的测试被称为正面结果。这可能有点令人困惑，因为患病似乎并不是一件积极的事情。而y = 1表示“正面”案例。为了保持清晰，确认你对y = 1在你数据背景下的理解是个好主意。

scikit-learn有一个函数来计算和绘制混淆矩阵：

`from` `sklearn``.``metrics` `import` `confusion_matrix`
`mat` `=` `confusion_matrix``(``y``,` `lr_model``.``predict``(``X``)``)`
`mat`

array([[377,  49],
       [104, 129]])

理想情况下，我们希望在对角方格中看到所有计数 True negative 和 True positive。这意味着我们已经正确分类了所有内容。但这很少见，我们需要评估错误的规模。为此，比较率而不是计数更容易。接下来，我们描述不同的率以及何时可能更喜欢优先考虑其中之一。

精度与召回率

在某些情况下，错过阳性案例的成本可能会更高。例如，如果我们正在构建一个用于识别肿瘤的分类器，我们希望确保不会错过任何恶性肿瘤。相反，我们不太关心将良性肿瘤分类为恶性，因为病理学家仍然需要仔细查看以验证恶性分类。在这种情况下，我们希望在实际上为阳性的记录中具有很高的真阳性率。该率称为敏感度或召回率：

Recall = True Positives True Positives + False Negatives = True Positives Actually True

较高的召回率会冒着将假记录预测为真的风险（假阳性）。

另一方面，当将电子邮件分类为垃圾邮件（阳性）或非垃圾邮件（阴性）时，如果一封重要的电子邮件被放入垃圾邮件文件夹中，我们可能会感到烦恼。在这种情况下，我们希望有高的精度，即模型对于阳性预测的准确性：

Precision = True Positives True Positives + False Positives = True Positives Predicted True

较高精度的模型通常更有可能预测真实观察结果为负（更高的假阴性率）。

常见的分析比较不同阈值下的精度和召回率：

`from` `sklearn` `import` `metrics`
`precision``,` `recall``,` `threshold` `=` `(`
    `metrics``.``precision_recall_curve``(``y``,` `lr_model``.``predict_proba``(``X``)``[``:``,` `1``]``)``)`

`tpr_df` `=` `pd``.``DataFrame``(``{``"``threshold``"``:``threshold``,` 
                       `"``precision``"``:``precision``[``:``-``1``]``,` `"``recall``"``:` `recall``[``:``-``1``]``,` `}``)`

为了查看精度和召回率之间的关系，我们将它们都绘制在阈值τ上：

另一个评估分类器性能的常见图表是精度-召回率曲线，简称 PR 曲线。它绘制了每个阈值的精度-召回率对：

`fig` `=` `px``.``line``(``tpr_df``,` `x``=``"``recall``"``,` `y``=``"``precision``"``,`
              `labels``=``{``"``recall``"``:``"``Recall``"``,``"``precision``"``:``"``Precision``"``}``)`
`fig``.``update_layout``(``width``=``450``,` `height``=``250``,` `yaxis_range``=``[``0``,` `1``]``)`
`fig`

注意，曲线的右端反映了样本中的不平衡性。精度与样本中倒下树木的比例相匹配，为 0.35。为不同模型绘制多个 PR 曲线可以帮助比较模型。

使用精度和召回率使我们能够更好地控制哪种类型的错误更重要。例如，假设我们想确保至少 75%的倒下树木被分类为倒下。我们可以找到发生这种情况的阈值：

`fall75_ind` `=` `np``.``argmin``(``recall` `>``=` `0.75``)` `-` `1`

`fall75_threshold` `=` `threshold``[``fall75_ind``]`
`fall75_precision` `=` `precision``[``fall75_ind``]`
`fall75_recall` `=` `recall``[``fall75_ind``]`

Threshold: 0.33
Precision: 0.59
Recall:    0.81

我们发现约 41%（1 - 精度）的我们分类为倒下的树实际上是站立的。此外，我们发现低于此阈值的树木比例为：

`print``(``"``Proportion of samples below threshold:``"``,` 
      `f``"``{``np``.``mean``(``lr_model``.``predict_proba``(``X``)``[``:``,``1``]` `<` `fall75_threshold``)``:``0.2f``}``"``)`

Proportion of samples below threshold: 0.52

因此，我们已将 52%的样本分类为站立（负面）。特异性（也称为真负率）衡量分类器将属于负类的数据标记为负类的比例：

Specificity = True Negatives True Negatives + False Positives = True Negatives Predicted False

我们的阈值的特异性为：

`act_neg` `=` `(``y` `==` `0``)`
`true_neg` `=` `(``lr_model``.``predict_proba``(``X``)``[``:``,``1``]` `<` `fall75_threshold``)` `&` `act_neg`

Specificity: 0.70

换句话说，70%的被分类为站立的树实际上是站立的。

如我们所见，有几种使用 2x2 混淆矩阵的方法。理想情况下，我们希望准确率、精确率和召回率都很高。这种情况发生在大多数预测落在表格的对角线上，因此我们的预测几乎全部正确——真负类和真正类。不幸的是，在大多数情况下，我们的模型会有一定程度的错误。在我们的例子中，相同直径的树木包括倒下的和站立的混合，因此我们不能完美地根据它们的直径分类树木。在实践中，当数据科学家选择一个阈值时，他们需要考虑自己的背景来决定是优先考虑精确率、召回率还是特异性。

总结

在本章中，我们用一个解释变量拟合简单的逻辑回归，但是我们可以通过将更多特征添加到我们的设计矩阵中来轻松地包含模型中的其他变量。例如，如果某些预测变量是分类的，我们可以将它们作为独热编码特征包含进来。这些想法直接延续自第十五章。正则化技术（来自第十六章）也适用于逻辑回归。我们将在第二十一章的案例研究中整合所有这些建模技术——包括使用训练集-测试集分割来评估模型和交叉验证来选择阈值——以开发一个用于分类假新闻的模型。

逻辑回归是机器学习中的基石，因为它自然地扩展到更复杂的模型中。例如，逻辑回归是神经网络的基本组成部分之一。当响应变量有多于两个类别时，逻辑回归可以扩展为多项逻辑回归。适用于建模计数的逻辑回归的另一个扩展称为泊松回归。这些不同形式的回归与最大似然密切相关，其中响应的潜在模型分别为二项式、多项式或泊松分布，目标是优化参数的数据似然。这些模型家族也被称为广义线性模型。在所有这些场景中，不存在用于最小化损失的封闭形式解决方案，因此平均损失的优化依赖于数值方法，我们将在下一章中介绍。

第二十章：数值优化

在本书的这一部分，我们的建模过程应该感到很熟悉：我们定义一个模型，选择一个损失函数，并通过最小化训练数据上的平均损失来拟合模型。我们已经看到了几种最小化损失的技术。例如，我们在第十五章中使用了微积分和几何论证，找到了使用平方损失拟合线性模型的简单表达式。

但经验损失最小化并不总是那么简单。Lasso 回归在平均平方损失中加入L 1惩罚后，不再有闭合形式的解，而逻辑回归使用交叉熵损失来拟合非线性模型。在这些情况下，我们使用数值优化来拟合模型，系统地选择参数值以评估平均损失，以寻找最小化的值。

当我们在第四章介绍损失函数时，我们执行了一个简单的数值优化，以找到平均损失的最小化者。我们创建了一个θ值的网格，并在网格中的所有点上评估平均损失（见图 20-1）。具有最小平均损失的网格点我们认为是最佳拟合。不幸的是，这种网格搜索很快变得不切实际，原因如下：

• 对于具有许多特征的复杂模型，网格变得难以管理。对于仅具有四个特征和每个特征 100 个值的网格，我们必须评估100 4 = 100,000,000个网格点上的平均损失。

• 必须事先指定要搜索的参数值范围，以创建网格。当我们对范围没有很好的感觉时，我们需要从一个宽网格开始，可能需要在更窄的范围内重复网格搜索。

• 对于大量观测值，评估网格点上的平均损失可能会很慢。

图 20-1. 在网格点上搜索可能计算速度慢或不准确

在本章中，我们介绍利用损失函数的形状和平滑性寻找最小化参数值的数值优化技术。我们首先介绍梯度下降技术的基本思想，然后给出一个示例，描述使梯度下降起作用的损失函数的特性，最后提供了梯度下降的几个扩展。

梯度下降基础知识

梯度下降基于以下观念：对于许多损失函数，函数在参数的小邻域内大致是线性的。图 20-2 给出了这个基本思想的示意图。

图 20-2. 梯度下降技术是向最小化参数值的方向进行小增量移动的技术。

在图中，我们画出了损失曲线L在最小化值θ ^ 左侧某点的切线。注意到切线的斜率为负。在θ 右侧前进一个短距离θ + s （其中s 是一小量）给出切线上接近θ + s 处的损失，且该损失小于L ( θ ~ ) 。也就是说，由于斜率b 为负，并且切线在θ 附近近似损失函数，我们有：

L ( θ + s ) ≈ L ( θ ) + b × s < L ( θ )

因此，向θ 的右侧小步骤会减小损失。另一方面，在图 Figure 20-2 中θ ^ 的另一侧，切线是正的，向左侧小步骤也会减小损失。

当我们根据切线斜率的正负指示重复采取小步骤时，这导致平均损失值越来越小，最终使我们接近或达到最小化值θ ^ 。这就是梯度下降背后的基本思想。

更正式地说，为了最小化一般参数向量θ 的损失函数L ( θ ) ，梯度（一阶偏导数）决定了应该采取的步长和方向。如果我们将梯度∇ θ L ( θ ) 简单写作g ( θ ) ，那么梯度下降法指出每次增量或步长为− α g ( θ ) ，其中α 是某个小正数。然后，在新位置的平均损失是：

L ( θ + ( − α g ( θ ) ) ≈ L ( θ ) − α g ( θ ) T g ( θ ) < L ( θ )

注意 g ( θ ) 是一个 p × 1 向量，而 g ( θ ) T g ( θ ) 是正的。

梯度下降算法的步骤如下：

1. 选择一个起始值，称为 θ ( 0 )（一个常见的选择是 θ ( 0 ) = 0）。

2. 计算 θ ( t + 1 ) = θ ( t ) − α g ( θ ) 。

3. 重复步骤 2 直到 θ ( t + 1 ) 不再改变（或变化很小）为止。

数量 α 被称为学习率。设置 α 可能比较棘手。它需要足够小以避免超过最小值，但足够大以在合理步数内达到最小值（见 图 20-3）。有许多设置 α 的策略。例如，随着时间的推移减少 α 可能会很有用。当 α 在迭代之间变化时，我们使用符号 α ( t ) 表示学习率在搜索过程中变化。

图 20-3. 较小的学习率需要许多步骤才能收敛（左），较大的学习率可能会发散（右）；选择适当的学习率可以快速收敛到最小值（中）。

梯度下降算法简单而强大，因为我们可以用它来拟合许多类型的模型和许多类型的损失函数。它是拟合许多模型的计算工具首选，包括大数据集上的线性回归和逻辑回归。接下来，我们演示使用该算法来拟合巴士延误数据中的常量（来自 第四章）。

最小化 Huber 损失

Huber 损失 结合了绝对损失和平方损失，得到一个既可微（像平方损失）又对异常值不那么敏感（像绝对损失）的函数：

L ( θ , y ) = 1 n ∑ i = 1 n { 1 2 ( y i − θ ) 2 | y i − θ | ≤ γ γ ( | y i − θ | − 1 2 γ ) otherwise

由于 Huber 损失是可微的，我们可以使用梯度下降。我们首先找到平均 Huber 损失的梯度：

∇ θ L ( θ , y ) = 1 n ∑ i = 1 n { − ( y i − θ ) | y i − θ | ≤ γ − γ ⋅ sign ( y i − θ ) otherwise

我们创建了 huber_loss 和 grad_huber_loss 函数来计算平均损失及其梯度。我们编写这些函数时签名设计使我们能够指定参数以及我们平均的观察数据和损失函数的转折点：

`def` `huber_loss``(``theta``,` `dataset``,` `gamma``=``1``)``:`
    `d` `=` `np``.``abs``(``theta` `-` `dataset``)`
    `return` `np``.``mean``(`
        `np``.``where``(``d` `<``=` `gamma``,`
                 `(``theta` `-` `dataset``)``*``*``2` `/` `2.0``,`
                 `gamma` `*` `(``d` `-` `gamma` `/` `2.0``)``)`
    `)`

`def` `grad_huber_loss``(``theta``,` `dataset``,` `gamma``=``1``)``:`
    `d` `=` `np``.``abs``(``theta` `-` `dataset``)`
    `return` `np``.``mean``(`
        `np``.``where``(``d` `<``=` `gamma``,`
                 `-``(``dataset` `-` `theta``)``,`
                 `-``gamma` `*` `np``.``sign``(``dataset` `-` `theta``)``)`
    `)`

接下来，我们编写了梯度下降的简单实现。我们的函数签名包括损失函数、其梯度和要平均的数据。我们还提供学习率。

`def` `minimize``(``loss_fn``,` `grad_loss_fn``,` `dataset``,` `alpha``=``0.2``,` `progress``=``False``)``:`
    `'''`
 `Uses gradient descent to minimize loss_fn. Returns the minimizing value of`
 `theta_hat once theta_hat changes less than 0.001 between iterations.`
 `'''`
    `theta` `=` `0`
    `while` `True``:`
        `if` `progress``:`
            `print``(``f``'``theta:` `{``theta``:``.2f``}` `| loss:` `{``loss_fn``(``theta``,` `dataset``)``:``.3f``}``'``)`
        `gradient` `=` `grad_loss_fn``(``theta``,` `dataset``)`
        `new_theta` `=` `theta` `-` `alpha` `*` `gradient`

        `if` `abs``(``new_theta` `-` `theta``)` `<` `0.001``:`
            `return` `new_theta`

        `theta` `=` `new_theta`

请回想一下，公交车延误数据集包含超过 1,000 个测量值，即北行 C 线公交车在抵达西雅图第三大道和派克街站点时晚多少分钟：

`delays` `=` `pd``.``read_csv``(``'``data/seattle_bus_times_NC.csv``'``)`

在 第四章 中，我们为这些数据拟合了一个常数模型，以得到绝对损失和平方损失。我们发现绝对损失产生了数据的中位数，而平方损失产生了数据的均值：

`print``(``f``"``Mean:` `{``np``.``mean``(``delays``[``'``minutes_late``'``]``)``:``.3f``}``"``)`
`print``(``f``"``Median:` `{``np``.``median``(``delays``[``'``minutes_late``'``]``)``:``.3f``}``"``)`    

Mean:   1.920
Median: 0.742

现在我们使用梯度下降算法来找到最小化 Huber 损失的常数模型：

`%``%``time`
`theta_hat` `=` `minimize``(``huber_loss``,` `grad_huber_loss``,` `delays``[``'``minutes_late``'``]``)`
`print``(``f``'``Minimizing theta:` `{``theta_hat``:``.3f``}``'``)`
`print``(``)`

Minimizing theta: 0.701

CPU times: user 93 ms, sys: 4.24 ms, total: 97.3 ms
Wall time: 140 ms

Huber 损失的优化常数接近最小化绝对损失的值。这是由于 Huber 损失函数的形状决定的。它在尾部是线性的，因此不像绝对损失那样受到异常值的影响，也不像平方损失那样受到影响。

警告

我们编写了我们的 minimize 函数来演示算法背后的思想。在实践中，您应该使用经过充分测试的、数值上稳定的优化算法的实现。例如，scipy 包中有一个 minimize 方法，我们可以用它来找到平均损失的最小化器，甚至不需要计算梯度。这个算法可能比我们可能编写的任何一个算法都要快得多。事实上，我们在 第十八章 中使用它来创建我们自己的二次损失的非对称修改，特别是在我们希望损失对于最小值的一侧的错误更大而另一侧的影响较小的特殊情况下。

更一般地，我们通常在迭代之间θ ( t ) 不太改变时停止算法。在我们的函数中，我们停止当θ ( t + 1 ) − θ ( t ) 小于 0.001 时。在迭代次数较多时，例如 1,000 次，停止搜索也很常见。如果算法在 1,000 次迭代后仍未到达最小值，则可能是因为学习率过大，或者最小值可能存在于极限处± ∞。

当我们无法通过解析方法轻松求解最小值或者最小化计算成本很高时，梯度下降给了我们一个通用的最小化平均损失的方法。该算法依赖于平均损失函数的两个重要属性：它在θ 上既是凸的又是可微的。我们接下来讨论算法如何依赖这些属性。

凸和可微损失函数

正如其名字所示，梯度下降算法要求被最小化的函数是可微的。梯度∇ θ L ( θ ) ，使我们能够对θ 的小邻域内的平均损失进行线性近似。这个近似给出了我们的步长的方向（和大小），只要我们不超过最小值θ ^ ，我们就一定会最终到达它。嗯，只要损失函数也是凸的。

最小值的逐步搜索也依赖于损失函数是凸函数。左图中的函数是凸的，但右图中的函数不是。右图中的函数有一个局部最小值，根据算法开始的位置，它可能会收敛到这个局部最小值并完全错过真正的最小值。凸性质避免了这个问题。凸函数避免了局部最小值的问题。所以，通过适当的步长，梯度下降可以找到任何凸、可微函数的全局最优解θ。

图 20-4。对于非凸函数（右图），梯度下降可能会找到局部最小值而不是全局最小值，而凸函数（左图）不可能出现这种情况。

形式上，函数 f 如果对于任意的输入值 θ a 和 θ b，以及介于 0 和 1 之间的任意 q 都是凸的：

q f ( θ a ) + ( 1 − q ) f ( θ b ) ≥ f ( q θ a + ( 1 − q ) θ b )

这个不等式意味着连接函数的任意两个点的线段必须位于或位于函数本身之上。从启发式的角度来看，这意味着无论我们在梯度为负时向右走还是在梯度为正时向左走，只要我们采取足够小的步伐，我们就会朝向函数的最小值方向前进。

凸性的正式定义为我们提供了确定一个函数是否为凸函数的精确方式。我们可以利用这个定义来将平均损失函数 L ( θ ) 的凸性与损失函数 l ( θ ) 连接起来。迄今为止，在本章中，我们通过不提及数据来简化了 L ( θ ) 的表示。回顾一下：

L ( θ , X , y ) = 1 n ∑ i = 1 n l ( θ , x i , y i )

其中 X 是一个 n × p 设计矩阵，x i 是设计矩阵的第 i 行，对应于数据集中的第 i 个观测值。这意味着梯度可以表示为：

∇ θ L ( θ , X , y ) = 1 n ∑ i = 1 n ∇ θ l ( θ , x i , y i )

如果l ( θ , x i , y i )是关于θ的凸函数，则平均损失也是凸的。对于导数也是类似的：l ( θ , x i , y i )的导数被平均到数据上以评估L ( θ , X , y )的导数。我们将在练习中详细讨论凸性质。

现在，有了大量数据，计算θ ( t )可能会非常耗时，因为它涉及到所有( x i , y i )上的梯度∇ θ l的平均值。接下来，我们考虑梯度下降的变体，这些变体可以更快地计算，因为它们不会对所有数据进行平均。

梯度下降的变体

梯度下降的两个变体，随机梯度下降和小批量梯度下降，在计算平均损失的梯度时使用数据子集，并且对于具有大型数据集的优化问题很有用。第三个选择是牛顿法，它假设损失函数是两次可微的，并且使用损失函数的二次近似，而不是梯度下降中使用的线性近似。

回顾一下，梯度下降是根据梯度采取步骤的。在步骤t，我们从θ ( t )移动到：

θ ( t + 1 ) = θ ( t ) − α ⋅ ∇ θ L ( θ ( t ) , X , y )

由于 ∇ θ L ( θ , X , y ) 可以表示为损失函数 l 的平均梯度，我们有：

∇ θ L ( θ , X , y ) = 1 n ∑ i = 1 n ∇ θ l ( θ , x i , y i )

这种根据数据中每个点处损失的梯度的平均来表示平均损失梯度的方法说明了为什么这个算法也被称为批梯度下降。批梯度下降的两个变体使用较小数量的数据而不是完整的“批次”。第一个，随机梯度下降，在算法的每一步中只使用一个观察。

随机梯度下降

尽管批梯度下降通常能在相对较少的迭代中找到最优的 θ，但如果数据集包含许多观察结果，每次迭代可能需要很长时间来计算。为了克服这个困难，随机梯度下降通过单个、随机选择的数据点来近似整体梯度。由于此观察是随机选择的，我们期望使用在随机选择观察点处的梯度，平均而言会朝着正确的方向移动，从而最终收敛到最小化参数。

简而言之，为了进行随机梯度下降，我们将平均梯度替换为单个数据点处的梯度。因此，更新后的公式就是：

θ ( t + 1 ) = θ ( t ) − α ⋅ ∇ θ l ( θ ( t ) , x i , y i )

在这个公式中，i t h 观测值 ( x i , y i ) 是从数据中随机选择的。随机选择点对于随机梯度下降的成功至关重要。如果点不是随机选择的，算法可能会产生比批梯度下降更差的结果。

我们通常通过随机重新排列所有数据点并按照它们的重新排列顺序使用每个点，直到完成一整个数据的遍历来运行随机梯度下降。如果算法尚未收敛，那么我们会重新洗牌数据并再次遍历数据。每个迭代的随机梯度下降看一个数据点；每个完整的数据遍历称为epoch。

由于随机下降每次只检查一个数据点，有时会朝着极小化器θ ^的方向迈出步伐，但平均而言这些步骤是朝着正确的方向。由于该算法的更新速度比批量梯度下降快得多，因此它可以在批量梯度下降完成单次更新时朝着最优θ ^取得显著进展。

小批量梯度下降

正如其名称所示，小批量梯度下降 在批量梯度下降和随机梯度下降之间取得平衡，通过在每次迭代中随机选择更多的观测值来增加样本数。在小批量梯度下降中，我们对少量数据点的损失函数梯度进行平均，而不是单个点或所有点的梯度。我们让B表示从数据集中随机抽取的小批次数据点，并定义算法的下一步为：

θ ( t + 1 ) = θ ( t ) − α ⋅ 1 | B | ∑ i ∈ B ∇ θ l ( θ , x i , y i )

与随机梯度下降类似，我们通过随机洗牌数据来执行小批量梯度下降。然后我们将数据分割成连续的小批次，并按顺序迭代这些批次。每个 epoch 后，重新洗牌数据并选择新的小批次。

虽然我们已经区分了随机梯度下降和小批量梯度下降，随机梯度下降 有时被用作一个总称，包括任意大小的小批次的选择。

另一种常见的优化技术是牛顿法。

牛顿法

牛顿法利用二阶导数优化损失。其基本思想是在L ( θ ) 的小邻域内，用二次曲线而不是线性逼近来近似平均损失。对于一个小步长s，逼近如下所示：

L ( θ + s ) ≈ L ( θ ) + g ( θ ) T s + 1 2 s T H ( θ ) s

其中g ( θ ) = ∇ θ L ( θ )是梯度，H ( θ ) = ∇ θ 2 L ( θ )是L ( θ )的海森矩阵。更具体地说，H是θ中的二阶偏导数的p × p矩阵，具有i，j元素：

H i , j = ∂ 2 l ∂ θ i ∂ θ j

对L ( θ + s )的二次逼近在s = − [ H − 1 ( θ ) ] g ( θ ) 处具有最小值。（凸性意味着H是对称方阵，可以被反转。）然后算法中的一步从θ ( t )移动到：

θ ( t + 1 ) = θ ( t ) + 1 n ∑ i = 1 n − [ H − 1 ( θ ( t ) ] g ( θ ( t ) )

图 20-5 展示了牛顿法优化的思想。

图 20-5。牛顿法使用对曲线的局部二次逼近来朝着凸、两次可微函数的最小值迈出步伐

此技术在逼近准确且步长小的情况下会快速收敛。否则，牛顿法可能会发散，这通常发生在函数在某个维度上几乎平坦的情况下。当函数相对平坦时，导数接近于零，其倒数可能非常大。大步长可能会移动到离逼近准确点很远的θ处。（与梯度下降不同，没有学习率可以保持步长小。）

摘要

在本章中，我们介绍了几种利用损失函数的形状和平滑性进行数值优化的技术，以搜索最小化参数值。我们首先介绍了梯度下降，它依赖于损失函数的可微性。梯度下降，也称为批量梯度下降，通过迭代改善模型参数，直到模型达到最小损失。由于批量梯度下降在处理大数据集时计算复杂度高，我们通常改用随机梯度下降来拟合模型。

小批量梯度下降在运行在某些计算机上找到的图形处理单元（GPU）芯片时最为优化。由于这些硬件类型可以并行执行计算，使用小批量可以提高梯度的准确性，而不增加计算时间。根据 GPU 的内存大小，小批量大小通常设置在 10 到 100 个观测之间。

或者，如果损失函数是二次可微的，则牛顿法可以非常快速地收敛，尽管在迭代中计算一步较为昂贵。混合方法也很受欢迎，先用梯度下降（某种类型），然后切换算法至牛顿法。这种方法可以避免发散，并且比单独使用梯度下降更快。通常，在最优点附近，牛顿法使用的二阶近似更为合适且收敛速度快。

最后，另一个选项是自适应设置步长。此外，如果不同特征的规模不同或频率不同，则设置不同的学习率可能很重要。例如，单词计数在常见单词和罕见单词之间可能会有很大差异。

在第十九章介绍的逻辑回归模型是通过本章描述的数值优化方法拟合的。我们最后介绍了一个案例研究，使用逻辑回归来拟合一个具有数千个特征的复杂模型。

第二十一章：案例研究：检测假新闻

假新闻——为了欺骗他人而创造的虚假信息——是一个重要问题，因为它可能会伤害人们。例如，社交媒体上的帖子在 图 21-1 中自信地宣称手部消毒剂对冠状病毒无效。尽管事实不确，但它还是通过社交媒体传播开来：被分享了近 10 万次，很可能被数百万人看到。

图 21-1. 2020 年 3 月推特上流行的一条帖子错误地声称消毒剂不能杀死冠状病毒。

我们可能会想知道是否可以在不阅读故事的情况下自动检测假新闻。对于这个案例研究，我们遵循数据科学生命周期的步骤。我们首先细化我们的研究问题并获取新闻文章和标签的数据集。然后我们清理和转换数据。接下来，我们探索数据以理解其内容，并设计用于建模的特征。最后，我们使用逻辑回归构建模型，预测新闻文章是否真实或虚假，并评估其性能。

我们包括这个案例研究是因为它让我们重申数据科学中的几个重要概念。首先，自然语言数据经常出现，即使是基本技术也能进行有用的分析。其次，模型选择是数据分析的重要部分，在这个案例研究中，我们应用了交叉验证、偏差-方差权衡和正则化的学习成果。最后，即使在测试集上表现良好的模型，在实际应用中也可能存在固有限制，我们很快就会看到。

让我们首先细化我们的研究问题，并理解我们数据的范围。

问题和范围

我们最初的研究问题是：我们能自动检测假新闻吗？为了细化这个问题，我们考虑了用于建立检测假新闻模型的信息类型。如果我们手动分类了新闻故事，人们已阅读每个故事并确定其真假，那么我们的问题变成了：我们能否建立一个模型来准确预测新闻故事是否为假的，基于其内容？

为了解决这个问题，我们可以使用 FakeNewsNet 数据库，如 Shu et al 所述。该数据库包含来自新闻和社交媒体网站的内容，以及用户参与度等元数据。为简单起见，我们只查看数据集的政治新闻文章。该数据子集仅包括由 Politifact 进行事实检查的文章，Politifact 是一个声誉良好的非党派组织。数据集中的每篇文章都有基于 Politifact 评估的“真实”或“虚假”标签，我们将其用作基准真实性。

Politifact 使用非随机抽样方法选择文章进行事实核查。根据其网站，Politifact 的记者每天选择“最有新闻价值和重要性”的主张。Politifact 始于 2007 年，存储库发布于 2020 年，因此大多数文章发布于 2007 年到 2020 年之间。

总结这些信息，我们确定目标人群包括所有在线发布的政治新闻故事，时间跨度从 2007 年到 2020 年（我们也想列出这些故事的来源）。访问框架由 Politifact 确定，标识出当天最有新闻价值的主张。因此，这些数据的主要偏见来源包括：

覆盖偏见

新闻媒体仅限于 Politifact 监控的那些，这可能会忽略奥秘或短暂存在的网站。

选择偏见

数据仅限于 Politifact 认为足够有趣以进行事实核查的文章，这意味着文章可能偏向于广泛分享和具有争议性的文章。

测量偏见

故事是否应标记为“假”或“真”由一个组织（Politifact）决定，并反映了该组织在其事实核查方法中存在的偏见，无论是有意还是无意。

漂移

由于我们只有 2007 年到 2020 年间发布的文章，内容可能会有些漂移。话题在快速发展的新闻趋势中被推广和伪造。

我们在开始整理数据之前，会牢记这些数据的限制，以便将其整理成可分析的形式。

获取和整理数据

让我们使用FakeNewsNet 的 GitHub 页面将数据导入 Python。阅读存储库描述和代码后，我们发现该存储库实际上并不存储新闻文章本身。相反，运行存储库代码将直接从在线网页上抓取新闻文章（使用我们在第十四章中介绍的技术）。这带来了一个挑战：如果一篇文章不再在网上可用，那么它很可能会在我们的数据集中丢失。注意到这一点后，让我们继续下载数据。

注意

FakeNewsNet 代码突显了可重复研究中的一个挑战——在线数据集随时间变化，但如果在存储库中存储和共享这些数据可能会面临困难（甚至违法）。例如，FakeNewsNet 数据集的其他部分使用 Twitter 帖子，但如果创建者在其存储库中存储帖子副本则会违反 Twitter 的条款和服务。在处理从网络收集的数据时，建议记录数据收集日期并仔细阅读数据来源的条款和服务。

运行脚本下载 Politifact 数据大约需要一个小时。之后，我们将数据文件放入data/politifact文件夹中。Politifact 标记为假和真的文章分别位于data/politifact/fake和data/politifact/real文件夹中。让我们看一看其中一个标记为“真实”的文章：

`!`ls -l data/politifact/real `|` head -n `5`

total 0
drwxr-xr-x  2 sam  staff  64 Jul 14  2022 politifact100
drwxr-xr-x  3 sam  staff  96 Jul 14  2022 politifact1013
drwxr-xr-x  3 sam  staff  96 Jul 14  2022 politifact1014
drwxr-xr-x  2 sam  staff  64 Jul 14  2022 politifact10185
ls: stdout: Undefined error: 0

`!`ls -lh data/politifact/real/politifact1013/

total 16
-rw-r--r--  1 sam  staff   5.7K Jul 14  2022 news content.json

每篇文章的数据存储在名为 news content.json 的 JSON 文件中。让我们将一篇文章的 JSON 加载到 Python 字典中（参见 第十四章）：

`import` `json`
`from` `pathlib` `import` `Path`

`article_path` `=` `Path``(``'``data/politifact/real/politifact1013/news content.json``'``)`
`article_json` `=` `json``.``loads``(``article_path``.``read_text``(``)``)`

这里，我们将 article_json 中的键和值显示为表格：

	value
key
---	---
url	www.senate.gov/legislative…...
text	Roll Call Vote 111th Congress - 1st Session\n...
images	[statse.webtrendslive.com/dcs222dj3ow…...
top_img	www.senate.gov/resources/i…
keywords	[]
authors	[]
canonical_link
title	U.S. Senate: U.S. Senate Roll Call Votes 111th...
meta_data	{'viewport’: ‘width=device-width, initial-scal...
movies	[]
publish_date	None
source	www.senate.gov
summary

JSON 文件中有很多字段，但是对于这个分析，我们只关注几个与文章内容主要相关的字段：文章的标题、文本内容、URL 和发布日期。我们创建一个数据框，其中每一行代表一篇文章（新闻报道的粒度）。为此，我们将每个可用的 JSON 文件加载为 Python 字典，然后提取感兴趣的字段以存储为 pandas 的 DataFrame，命名为 df_raw：

`from` `pathlib` `import` `Path`

`def` `df_row``(``content_json``)``:`
    `return` `{`
        `'``url``'``:` `content_json``[``'``url``'``]``,`
        `'``text``'``:` `content_json``[``'``text``'``]``,`
        `'``title``'``:` `content_json``[``'``title``'``]``,`
        `'``publish_date``'``:` `content_json``[``'``publish_date``'``]``,`
    `}`

`def` `load_json``(``folder``,` `label``)``:`
    `filepath` `=` `folder` `/` `'``news content.json``'`
    `data` `=` `df_row``(``json``.``loads``(``filepath``.``read_text``(``)``)``)` `if` `filepath``.``exists``(``)` `else` `{``}`
    `return` `{`
        `*``*``data``,`
        `'``label``'``:` `label``,`
    `}`

`fakes` `=` `Path``(``'``data/politifact/fake``'``)`
`reals` `=` `Path``(``'``data/politifact/real``'``)`

`df_raw` `=` `pd``.``DataFrame``(``[``load_json``(``path``,` `'``fake``'``)` `for` `path` `in` `fakes``.``iterdir``(``)``]` `+`
                      `[``load_json``(``path``,` `'``real``'``)` `for` `path` `in` `reals``.``iterdir``(``)``]``)`

`df_raw``.``head``(``2``)`

	url	text	title	publish_date	label
0	dailybuzzlive.com/cannibals-arrested-florida/	Police in Vernal Heights, Florida, arrested 3-...	Cannibals Arrested in Florida Claim Eating Hum...	1.62e+09	fake
1	web.archive.org/web/2017122…...	WASHINGTON — Rod Jay Rosenstein, Deputy Attorn...	BREAKING: Trump fires Deputy Attorney General ...	1.45e+09	fake

探索这个数据框会揭示一些在开始分析之前我们想要解决的问题。例如：

• 一些文章无法下载。当出现这种情况时，url 列包含 NaN。

• 一些文章没有文本（例如只有视频内容的网页）。我们从数据框中删除这些文章。

• publish_date 列以 Unix 格式（自 Unix 纪元以来的秒数）存储时间戳，因此我们需要将它们转换为 pandas.Timestamp 对象。

• 我们对网页的基本 URL 感兴趣。然而，JSON 文件中的 source 字段与 url 列相比有许多缺失值，所以我们必须使用 url 列中的完整 URL 提取基本 URL。例如，从 dailybuzzlive.com/cannibals-arrested-florida/ 我们得到 dailybuzzlive.com。

• 一些文章是从存档网站（web.archive.org）下载的。当这种情况发生时，我们希望从原始的 URL 中提取实际的基本 URL，通过移除 web.archive.org 前缀。

• 我们希望将 title 和 text 列连接成一个名为 content 的单一列，其中包含文章的所有文本内容。

我们可以使用 pandas 函数和正则表达式来解决这些数据问题：

`import` `re`

`# [1], [2]`
`def` `drop_nans``(``df``)``:`
    `return` `df``[``~``(``df``[``'``url``'``]``.``isna``(``)` `|`
                `(``df``[``'``text``'``]``.``str``.``strip``(``)` `==` `'``'``)` `|` 
                `(``df``[``'``title``'``]``.``str``.``strip``(``)` `==` `'``'``)``)``]`

`# [3]`
`def` `parse_timestamps``(``df``)``:`
    `timestamp` `=` `pd``.``to_datetime``(``df``[``'``publish_date``'``]``,` `unit``=``'``s``'``,` `errors``=``'``coerce``'``)`
    `return` `df``.``assign``(``timestamp``=``timestamp``)`

`# [4], [5]`
`archive_prefix_re` `=` `re``.``compile``(``r``'``https://web.archive.org/web/``\``d+/``'``)`
`site_prefix_re` `=` `re``.``compile``(``r``'``(https?://)?(www``\``.)?``'``)`
`port_re` `=` `re``.``compile``(``r``'``:``\``d+``'``)`

`def` `url_basename``(``url``)``:`
    `if` `archive_prefix_re``.``match``(``url``)``:`
        `url` `=` `archive_prefix_re``.``sub``(``'``'``,` `url``)`
    `site` `=` `site_prefix_re``.``sub``(``'``'``,` `url``)``.``split``(``'``/``'``)``[``0``]`
    `return` `port_re``.``sub``(``'``'``,` `site``)`

`# [6]`
`def` `combine_content``(``df``)``:`
    `return` `df``.``assign``(``content``=``df``[``'``title``'``]` `+` `'` `'` `+` `df``[``'``text``'``]``)`

`def` `subset_df``(``df``)``:`
    `return` `df``[``[``'``timestamp``'``,` `'``baseurl``'``,` `'``content``'``,` `'``label``'``]``]`

`df` `=` `(``df_raw`
 `.``pipe``(``drop_nans``)`
 `.``reset_index``(``drop``=``True``)`
 `.``assign``(``baseurl``=``lambda` `df``:` `df``[``'``url``'``]``.``apply``(``url_basename``)``)`
 `.``pipe``(``parse_timestamps``)`
 `.``pipe``(``combine_content``)`
 `.``pipe``(``subset_df``)`
`)`

数据整理后，我们得到名为df的以下数据框架：

`df``.``head``(``2``)`

	timestamp	baseurl	content	label
0	2021-04-05 16:39:51	dailybuzzlive.com	佛罗里达州被捕的食人族声称吃...	假
1	2016-01-01 23:17:43	houstonchronicle-tv.com	突发新闻：特朗普解雇副检察...	假

现在我们已加载并清理了数据，可以进行探索性数据分析。

探索数据

我们正在探索的新闻文章数据集只是更大的 FakeNewsNet 数据集的一部分。因此，原始论文并未提供有关我们数据子集的详细信息。因此，为了更好地理解数据，我们必须自己进行探索。

在开始探索性数据分析之前，我们遵循标准做法，将数据分割为训练集和测试集。我们只使用训练集进行 EDA：

`from` `sklearn``.``model_selection` `import` `train_test_split`

`df``[``'``label``'``]` `=` `(``df``[``'``label``'``]` `==` `'``fake``'``)``.``astype``(``int``)`

`X_train``,` `X_test``,` `y_train``,` `y_test` `=` `train_test_split``(`
    `df``[``[``'``timestamp``'``,` `'``baseurl``'``,` `'``content``'``]``]``,` `df``[``'``label``'``]``,`
    `test_size``=``0.25``,` `random_state``=``42``,`
`)`

`X_train``.``head``(``2``)`

	timestamp	baseurl	content
164	2019-01-04 19:25:46	worldnewsdailyreport.com	中国月球车未发现美国...
28	2016-01-12 21:02:28	occupydemocrats.com	弗吉尼亚州共和党人要求学校检查...

让我们统计训练集中真假文章的数量：

`y_train``.``value_counts``(``)`

label
0    320
1    264
Name: count, dtype: int64

我们的训练集有 584 篇文章，实际文章比虚假文章多约 60 篇。接下来，我们检查这三个字段中是否存在缺失值：

`X_train``.``info``(``)`

<class 'pandas.core.frame.DataFrame'>
Index: 584 entries, 164 to 102
Data columns (total 3 columns):
 #   Column     Non-Null Count  Dtype         
---  ------     --------------  -----         
 0   timestamp  306 non-null    datetime64[ns]
 1   baseurl    584 non-null    object        
 2   content    584 non-null    object        
dtypes: datetime64ns, object(2)
memory usage: 18.2+ KB

时间戳几乎一半为空。如果我们在分析中使用它，这个特征将限制数据集。让我们仔细查看baseurl，它表示发布原始文章的网站。

探索出版商

要理解baseurl列，我们首先统计每个网站的文章数：

`X_train``[``'``baseurl``'``]``.``value_counts``(``)`

baseurl
whitehouse.gov               21
abcnews.go.com               20
nytimes.com                  17
                             ..
occupydemocrats.com           1
legis.state.ak.us             1
dailynewsforamericans.com     1
Name: count, Length: 337, dtype: int64

我们的训练集有 584 行，我们发现有 337 个独特的发布网站。这意味着数据集包含许多只有少数文章的出版物。每个网站发布的文章数量的直方图证实了这一点：

fig = px.histogram(X_train['baseurl'].value_counts(), width=450, height=250,
                   labels={"value": "Number of articles published at a URL"})

`fig``.``update_layout``(``showlegend``=``False``)`

此直方图显示，绝大多数网站（337 中的 261 个）在训练集中只有一篇文章，只有少数网站在训练集中有超过五篇文章。尽管如此，识别发布最多假或真文章的网站可能具有信息量。首先，我们找出发布最多假文章的网站：

接下来，我们列出发布最多真实文章的网站：

cnn.com 和 washingtonpost.com 出现在两个列表中。即使我们不知道这些网站的文章总数，我们可能会预期来自 yournewswire.com 的文章更有可能被标记为“假”，而来自 whitehouse.gov 的文章更有可能被标记为“真”。尽管如此，我们并不指望使用发布网站来预测文章的真实性会非常有效；数据集中大多数网站的文章数量实在太少。

接下来，让我们探索timestamp列，记录新闻文章的发布日期。

探索发布日期

将时间戳绘制在直方图上显示，大多数文章是在 2000 年之后发布的，尽管至少有一篇文章是在 1940 年之前发布的：

`fig` `=` `px``.``histogram``(`
    `X_train``[``"``timestamp``"``]``,`
    `labels``=``{``"``value``"``:` `"``Publication year``"``}``,` `width``=``550``,` `height``=``250``,`
`)`
`fig``.``update_layout``(``showlegend``=``False``)`

当我们更仔细地查看发布于 2000 年之前的新闻文章时，我们发现时间戳与文章的实际发布日期不符。这些日期问题很可能与网络爬虫从网页中收集到的不准确信息有关。我们可以放大直方图中 2000 年之后的区域：

`fig` `=` `px``.``histogram``(`
    `X_train``.``loc``[``X_train``[``"``timestamp``"``]` `>` `"``2000``"``,` `"``timestamp``"``]``,`
    `labels``=``{``"``value``"``:` `"``Publication year``"``}``,` `width``=``550``,` `height``=``250``,` 
`)`
`fig``.``update_layout``(``showlegend``=``False``)`

正如预期的那样，大多数文章是在 2007 年（Politifact 成立的年份）至 2020 年（FakeNewsNet 仓库发布的年份）之间发布的。但我们还发现，时间戳主要集中在 2016 年至 2018 年之间——这是有争议的 2016 年美国总统选举年及其后两年。这一发现进一步提示我们的分析局限性可能不适用于非选举年。

我们的主要目标是使用文本内容进行分类。接下来我们探索一些词频。

探索文章中的单词

我们想看看文章中使用的单词与文章是否被标记为“假”之间是否存在关系。一个简单的方法是查看像 军事 这样的单词，然后统计提到“军事”的文章中有多少被标记为“假”。对于 军事 来说，文章提到它的比例应该远高于或远低于数据集中的假文章比例 45%（数据集中假文章的比例：264/584）。

我们可以利用我们对政治话题的领域知识来挑选一些候选单词进行探索：

word_features = [
    # names of presidential candidates
    'trump', 'clinton',
    # congress words
    'state', 'vote', 'congress', 'shutdown',

    # other possibly useful words
    'military', 'princ', 'investig', 'antifa', 
    'joke', 'homeless', 'swamp', 'cnn', 'the'
]

然后，我们定义一个函数，为每个单词创建一个新特征，如果文章中出现该单词，则特征为True，否则为False：

`def` `make_word_features``(``df``,` `words``)``:`
    `features` `=` `{` `word``:` `df``[``'``content``'``]``.``str``.``contains``(``word``)` `for` `word` `in` `words` `}`
    `return` `pd``.``DataFrame``(``features``)`

这就像是单词存在的一种独热编码（参见第十五章）。我们可以使用这个函数进一步处理我们的数据，并创建一个包含每个选择的单词特征的新数据框架：

`df_words` `=` `make_word_features``(``X_train``,` `word_features``)`
`df_words``[``"``label``"``]` `=` `df``[``"``label``"``]`

`df_words``.``shape`

(584, 16)

`df_words``.``head``(``4``)`

	trump	clinton	state	vote	...	swamp	cnn	the	label
164	False	False	True	False	...	False	False	True	1
28	False	False	False	False	...	False	False	True	1
708	False	False	True	True	...	False	False	True	0
193	False	False	False	False	...	False	False	True	1

4 rows × 16 columns

现在我们可以找出这些文章中被标记为fake的比例。我们在以下图表中可视化了这些计算结果。在左图中，我们用虚线标记了整个训练集中fake文章的比例，这有助于我们理解每个单词特征的信息量—一个高信息量的单词将使得其点远离该线：

这个图表揭示了建模的一些有趣的考虑因素。例如，注意到单词antifa具有很高的预测性—所有提到单词antifa的文章都标记为fake。然而，antifa只出现在少数文章中。另一方面，单词the几乎出现在每篇文章中，但对于区分real和fake文章没有信息量，因为含有the的文章的比例与总体 fake 文章的比例相匹配。我们可能会更喜欢像vote这样的单词，它既有预测能力又出现在许多新闻文章中。

这个探索性分析使我们了解了我们的新闻文章发表的时间框架，数据中涵盖的广泛发布网站以及用于预测的候选词。接下来，我们为预测文章是真实还是虚假拟合模型。

建模

现在我们已经获得、清洗并探索了我们的数据，让我们拟合模型来预测文章是真实还是虚假。在本节中，我们使用逻辑回归因为我们面临二元分类问题。我们拟合了三种不同复杂度的模型。首先，我们拟合了一个仅使用单个手动选择的单词作为解释特征的模型。然后我们拟合了一个使用多个手动选择的单词的模型。最后，我们拟合了一个使用所有在训练集中的单词，并使用 tf-idf 转换向量化的模型（介绍见第十三章）。让我们从简单的单词模型开始。

单词模型

我们的探索性数据分析表明，单词vote与文章被标记为real或fake相关。为了验证这一点，我们使用一个二元特征拟合了逻辑回归模型：如果文章中出现单词vote则为1，否则为0。我们首先定义一个函数将文章内容转为小写：

`def` `lowercase``(``df``)``:`
    `return` `df``.``assign``(``content``=``df``[``'``content``'``]``.``str``.``lower``(``)``)`

对于我们的第一个分类器，我们只使用单词vote：

`one_word` `=` `[``'``vote``'``]`

我们可以将lowercase函数和来自我们的探索性数据分析的make_word_features函数链在一起成为一个scikit-learn管道。这提供了一种方便的方式来一次性地转换和拟合数据：

`from` `sklearn``.``pipeline` `import` `make_pipeline`
`from` `sklearn``.``linear_model` `import` `LogisticRegressionCV`
`from` `sklearn``.``preprocessing` `import` `FunctionTransformer`

`model1` `=` `make_pipeline``(`
    `FunctionTransformer``(``lowercase``)``,`
    `FunctionTransformer``(``make_word_features``,` `kw_args``=``{``'``words``'``:` `one_word``}``)``,`
    `LogisticRegressionCV``(``Cs``=``10``,` `solver``=``'``saga``'``,` `n_jobs``=``4``,` `max_iter``=``10000``)``,`
`)`

在使用时，前面的流水线将文章内容中的字符转换为小写，为每个感兴趣的单词创建一个二元特征的数据框，并使用L 2正则化在数据上拟合逻辑回归模型。另外，LogisticRegressionCV函数使用交叉验证（默认为五折）来选择最佳的正则化参数。（有关正则化和交叉验证的更多信息，请参见第十六章。）

让我们使用管道来拟合训练数据：

`%``%``time`

`model1``.``fit``(``X_train``,` `y_train``)`
`print``(``f``'``{``model1``.``score``(``X_train``,` `y_train``)``:``.1%``}` `accuracy on training set.``'``)`

64.9% accuracy on training set.
CPU times: user 110 ms, sys: 42.7 ms, total: 152 ms
Wall time: 144 ms

总体而言，单词分类器只能正确分类 65% 的文章。我们在训练集上绘制分类器的混淆矩阵，以查看它所犯的错误类型：

我们的模型经常将真实文章（0）错误地分类为虚假（1）。由于这个模型很简单，我们可以看一下这两种情况的概率：文章中是否包含vote这个词：

"vote" present: [[0.72 0.28]]
"vote" absent: [[0.48 0.52]]

当文章包含vote这个词时，模型会给出文章为真实的高概率，而当vote缺失时，概率略微倾向于文章为虚假。我们鼓励读者使用逻辑回归模型的定义和拟合系数来验证这一点：

`print``(``f``'``Intercept:` `{``log_reg``.``intercept_``[``0``]``:``.2f``}``'``)`
`[``[``coef``]``]` `=` `log_reg``.``coef_`
`print``(``f``'``"``vote``"` `Coefficient:` `{``coef``:``.2f``}``'``)`

Intercept: 0.08
"vote" Coefficient: -1.00

正如我们在第十九章中看到的那样，系数表示随着解释变量的变化而发生的几率变化的大小。对于像文章中的一个 0-1 变量这样的变量，这有着特别直观的含义。对于一个包含vote的文章，其为虚假的几率会减少一个因子exp ⁡ ( θ v o t e )，即：

`np``.``exp``(``coef``)` 

0.36836305405149367

注意

请记住，在这种建模场景中，标签0表示真实文章，标签1表示虚假文章。这可能看起来有点反直觉—我们说“真正的正例”是当模型正确预测一篇虚假文章为虚假时。在二元分类中，我们通常说“正面”的结果是指存在某种不寻常情况的结果。例如，测试结果为阳性的人可能会有这种疾病。

让我们通过引入额外的单词特征来使我们的模型变得更加复杂一些。

多单词模型

我们创建了一个模型，该模型使用了我们在训练集的探索性数据分析（EDA）中检查过的所有单词，除了the。让我们使用这 15 个特征来拟合一个模型：

`model2` `=` `make_pipeline``(`
    `FunctionTransformer``(``lowercase``)``,`
    `FunctionTransformer``(``make_word_features``,` `kw_args``=``{``'``words``'``:` `word_features``}``)``,`
    `LogisticRegressionCV``(``Cs``=``10``,` `solver``=``'``saga``'``,` `n_jobs``=``4``,` `max_iter``=``10000``)``,`
`)`

`%``%``time`

`model2``.``fit``(``X_train``,` `y_train``)`
`print``(``f``'``{``model2``.``score``(``X_train``,` `y_train``)``:``.1%``}` `accuracy on training set.``'``)`

74.8% accuracy on training set.
CPU times: user 1.54 s, sys: 59.1 ms, total: 1.6 s
Wall time: 637 ms

该模型比单词模型准确率高约 10 个百分点。从一个单词模型转换为一个 15 个单词模型仅获得 10 个百分点可能会有点令人惊讶。混淆矩阵有助于揭示出所犯错误的类型：

我们可以看到，这个分类器在准确分类真实文章方面做得更好。然而，在将虚假文章分类时，它会犯更多错误——有 59 篇虚假文章被错误分类为真实文章。在这种情况下，我们可能更关注将一篇文章误分类为虚假而实际上它是真实的。因此，我们希望有很高的精确度——正确预测为虚假文章的虚假文章比例：

`model1_precision` `=` `238` `/` `(``238` `+` `179``)`
`model2_precision` `=` `205` `/` `(``205` `+` `88``)`

`[``round``(``num``,` `2``)` `for` `num` `in` `[``model1_precision``,` `model2_precision``]``]`

[0.57, 0.7]

我们更大的模型中的精确度有所提高，但约有 30% 的被标记为虚假的文章实际上是真实的。让我们来看一下模型的系数：

通过观察它们的符号，我们可以快速解释系数。trump 和 investig 上的大正值表明模型预测包含这些词的新文章更有可能是虚假的。对于像 congress 和 vote 这样具有负权重的词来说情况相反。我们可以使用这些系数来比较文章是否包含特定词时的对数几率。

尽管这个更大的模型的表现比简单的单词模型更好，但我们不得不使用我们对新闻的知识手动挑选单词特征。如果我们漏掉了高度预测性的词怎么办？为了解决这个问题，我们可以使用 tf-idf 转换将所有文章中的所有单词合并起来。

使用 tf-idf 转换进行预测

对于第三个也是最后一个模型，我们使用了第十三章中的 term frequency-inverse document frequency (tf-idf) 转换来向量化训练集中所有文章的整个文本。回想一下，使用此转换，一篇文章被转换为一个向量，其中每个词的出现次数在任何 564 篇文章中都会出现。该向量由词在文章中出现的次数的归一化计数组成，除以该词的稀有度。tf-idf 对仅出现在少数文档中的词赋予更高的权重。这意味着我们的分类器用于预测的是训练集新闻文章中的所有单词。正如我们之前介绍 tf-idf 时所做的那样，首先我们移除停用词，然后对单词进行标记化，最后我们使用 scikit-learn 中的 TfidfVectorizer：

`tfidf` `=` `TfidfVectorizer``(``tokenizer``=``stemming_tokenizer``,` `token_pattern``=``None``)`

`from` `sklearn``.``compose` `import` `make_column_transformer`

`model3` `=` `make_pipeline``(`
    `FunctionTransformer``(``lowercase``)``,`
    `make_column_transformer``(``(``tfidf``,` `'``content``'``)``)``,`
    `LogisticRegressionCV``(``Cs``=``10``,`
                         `solver``=``'``saga``'``,`
                         `n_jobs``=``8``,`
                         `max_iter``=``1000``)``,`
    `verbose``=``True``,`
`)`

`%``%``time`

`model3``.``fit``(``X_train``,` `y_train``)`
`print``(``f``'``{``model3``.``score``(``X_train``,` `y_train``)``:``.1%``}` `accuracy on training set.``'``)`

[Pipeline]  (step 1 of 3) Processing functiontransformer, total=   0.0s
[Pipeline] . (step 2 of 3) Processing columntransformer, total=  14.5s
[Pipeline]  (step 3 of 3) Processing logisticregressioncv, total=   6.3s
100.0% accuracy on training set.
CPU times: user 50.2 s, sys: 508 ms, total: 50.7 s
Wall time: 34.2 s

我们发现这个模型在训练集上实现了 100% 的准确率。我们可以查看 tf-idf 转换器来更好地理解模型。让我们首先找出分类器使用的唯一标记的数量：

`tfidf` `=` `model3``.``named_steps``.``columntransformer``.``named_transformers_``.``tfidfvectorizer`
`n_unique_tokens` `=` `len``(``tfidf``.``vocabulary_``.``keys``(``)``)`
`print``(``f``'``{``n_unique_tokens``}` `tokens appeared across` `{``len``(``X_train``)``}` `examples.``'``)`

23800 tokens appeared across 584 examples.

这意味着我们的分类器有 23,812 个特征，比我们之前的模型大幅增加，之前的模型只有 15 个。由于我们无法显示那么多模型权重，我们显示了 10 个最负和 10 个最正的权重：

这些系数展示了该模型的一些怪癖。我们看到一些有影响力的特征对应于原始文本中的标点符号。目前尚不清楚我们是否应该清除模型中的标点符号。一方面，标点符号似乎没有单词传达的意义那么多。另一方面，似乎合理的是，例如，一篇文章中有很多感叹号可能有助于模型决定文章是真实还是假的。在这种情况下，我们决定保留标点符号，但是好奇的读者可以在去除标点符号后重复此分析，以查看生成的模型受到的影响。

最后，我们显示了所有三个模型的测试集误差：

	测试集误差
模型 1	0.61
模型 2	0.70
模型 3	0.88

正如我们所预料的那样，随着我们引入更多的特征，模型变得更加准确。使用 tf-idf 的模型比使用二进制手工选择的词特征的模型表现要好得多，但是它没有达到在训练集上获得的 100% 准确率。这说明了建模中的一种常见权衡：在给定足够的数据的情况下，更复杂的模型通常可以胜过更简单的模型，特别是在这种情况研究中，更简单的模型有太多的模型偏差而表现不佳的情况下。但是，复杂的模型可能更难解释。例如，我们的 tf-idf 模型有超过 20,000 个特征，这使得基本上不可能解释我们的模型如何做出决策。此外，与模型 2 相比，tf-idf 模型需要更长时间进行预测——它的速度慢了 100 倍。在决定使用哪种模型时，所有这些因素都需要考虑在内。

另外，我们需要注意我们的模型适用于什么。在这种情况下，我们的模型使用新闻文章的内容进行预测，这使得它们高度依赖于出现在训练集中的单词。然而，我们的模型可能不会在未来的新闻文章上表现得像在训练集中没有出现的单词那样好。例如，我们的模型使用 2016 年美国选举候选人的名字进行预测，但是它们不会知道要在 2020 或 2024 年纳入候选人的名字。为了在较长时间内使用我们的模型，我们需要解决这个漂移问题。

话虽如此，令人惊讶的是，一个逻辑回归模型在相对较少的特征工程（tf-idf）下也能表现良好。我们已经回答了我们最初的研究问题：我们的 tf-idf 模型在检测我们数据集中的假新闻方面表现出色，而且可能可以推广到训练数据覆盖的同一时间段内发布的其他新闻。

摘要

我们很快就要结束本章，从而结束这本书。我们从讨论数据科学生命周期开始这本书。让我们再次看看生命周期，在 图 21-2 中，以欣赏您所学到的一切。

图 21-2。数据科学生命周期的四个高级步骤，本书中我们详细探讨了每个步骤。

本案例研究逐步介绍了数据科学生命周期的每个阶段：

1. 许多数据分析从一个研究问题开始。本章中我们呈现的案例研究从询问我们是否可以创建自动检测假新闻模型开始。

2. 我们使用在线找到的代码将网页抓取到 JSON 文件中来获取数据。由于数据描述相对较少，我们需要清理数据以理解它。这包括创建新的特征来指示文章中特定词语的存在或缺失。

3. 我们的初步探索确定了可能对预测有用的单词。在拟合简单模型并探索它们的精确度和准确度后，我们进一步使用 tf-idf 转换文章，将每篇新闻文章转换为归一化的词向量。

4. 我们将向量化文本作为逻辑模型中的特征，并使用正则化和交叉验证拟合最终模型。最后，在测试集上找到拟合模型的准确度和精确度。

当我们像这样详细列出生命周期中的步骤时，步骤之间似乎流畅连接在一起。但现实是混乱的——正如图表所示，真实数据分析在各个步骤之间来回跳跃。例如，在我们的案例研究结束时，我们发现了可能促使我们重新访问生命周期早期阶段的数据清理问题。尽管我们的模型非常准确，但大部分训练数据来自 2016 年至 2018 年的时期，因此如果我们想要在该时间段之外的文章上使用它，就必须仔细评估模型的性能。

本质上，重要的是在数据分析的每个阶段牢记整个生命周期。作为数据科学家，你将被要求证明你的决策，这意味着你需要深入理解你的研究问题和数据。本书中的原则和技术将为你提供一套基础技能。在你的数据科学旅程中继续前进，我们建议你通过以下方式继续扩展你的技能：

• 重新审视本书的案例研究。首先复制我们的分析，然后深入探讨你对数据的疑问。

• 进行独立的数据分析。提出你感兴趣的研究问题，从网络中找到相关数据，并分析数据，看看数据与你的期望有多大匹配。这样做将使你对整个数据科学生命周期有第一手经验。

• 深入研究一个主题。我们在附加材料附录中提供了许多深入资源。选择你最感兴趣的资源，并深入了解。

世界需要像你这样能够利用数据得出结论的人，因此我们真诚地希望你能利用这些技能帮助他人制定有效的战略、打造更好的产品，并做出明智的决策。