想象一下预测骰子掷出 4 的结果(这是骰子的单数形式)。用纯粹的概率思维方法来解决这个问题,一个人简单地说骰子有六个面。我们假设每个面都是等可能的,所以掷出 4 的概率是 1/6,或 16.666%。
然而,一个狂热的统计学家可能会说:“不!我们需要掷骰子来获取数据。如果我们能掷出 30 次或更多次,而且我们掷得越多,就越好,那么我们才能有数据来确定掷出 4 的概率。” 如果我们假设骰子是公平的,这种方法可能看起来很愚蠢,但如果不是呢?如果是这样,收集数据是发现掷出 4 的概率的唯一方法。我们将在第三章中讨论假设检验。
假设你有一枚公平的硬币和一枚公平的六面骰子。你想找出硬币翻转为正面和骰子掷出 6 的概率。这是两个不同事件的两个单独概率,但我们想找出这两个事件同时发生的概率。这被称为联合概率。
这很容易理解,但为什么会这样呢?许多概率规则可以通过生成所有可能的事件组合来发现,这源自离散数学中的排列组合知识领域。对于这种情况,生成硬币和骰子之间的每种可能结果,将正面(H)和反面(T)与数字 1 到 6 配对。请注意,我在我们得到正面和一个 6 的结果周围放了星号“*”:
抛硬币和掷骰子时有 12 种可能的结果。我们感兴趣的只有一个,即“H6”,即得到一个正面和一个 6。因为只有一个结果符合我们的条件,而有 12 种可能的结果,所以得到一个正面和一个 6 的概率是 1/12。
非互斥事件又是什么呢,即可以同时发生的事件?让我们回到抛硬币和掷骰子的例子。得到正面或 6 的概率是多少?在你尝试将这些概率相加之前,让我们再次生成所有可能的结果,并突出显示我们感兴趣的结果:
为什么会这样?再次研究硬币翻转和骰子结果的组合,看看是否能找到一些可疑之处。注意当我们将概率相加时,我们在“H6”和“T6”中都重复计算了得到 6 的概率!如果这不清楚,请尝试找出得到正面或掷骰子得到 1 至 5 的概率:
嗯……稍微研究一下这些数字,问问咖啡是否真的是问题所在。再次注意,任何时候只有 0.5%的人口患有癌症。然而,65%的人口定期喝咖啡。如果咖啡会导致癌症,我们难道不应该有比 0.5%更高的癌症数字吗?难道不应该接近 65%吗?
这就是比例数字的狡猾之处。它们在没有任何给定上下文的情况下可能看起来很重要,而媒体头条肯定可以利用这一点来获取点击量:“新研究揭示 85%的癌症患者喝咖啡”可能会这样写。当然,这很荒谬,因为我们将一个常见属性(喝咖啡)与一个不常见属性(患癌症)联系起来。
人们很容易被条件概率搞混,因为条件的方向很重要,而且两个条件被混淆为相等。"在你是咖啡饮用者的情况下患癌症的概率"不同于"在你患癌症的情况下是咖啡饮用者的概率"。简单来说:很少有咖啡饮用者患癌症,但很多癌症患者喝咖啡。
如果你想在 Python 中计算这个问题,请查看示例 2-1。
示例 2-1. 在 Python 中使用贝叶斯定理
p_coffee_drinker = .65
p_cancer = .005
p_coffee_drinker_given_cancer = .85
p_cancer_given_coffee_drinker = p_coffee_drinker_given_cancer *
p_cancer / p_coffee_drinker
# prints 0.006538461538461539
print(p_cancer_given_coffee_drinker)
因此,某人在是咖啡饮用者的情况下患癌症的概率仅为 0.65%!这个数字与某人在患癌症的情况下是咖啡饮用者的概率(85%)非常不同。现在你明白为什么条件的方向很重要了吗?贝叶斯定理就是为了这个原因而有用。它还可以用于将多个条件概率链接在一起,根据新信息不断更新我们的信念。
如果你想更深入地探讨贝叶斯定理背后的直觉,请参考附录 A。现在只需知道它帮助我们翻转条件概率。接下来让我们讨论条件概率如何与联合和并集操作相互作用。
朴素贝叶斯
贝叶斯定理在一种常见的机器学习算法——朴素贝叶斯中扮演着核心角色。乔尔·格鲁斯在他的书《从零开始的数据科学》(O'Reilly)中对此进行了介绍。
联合和条件概率
让我们重新审视联合概率以及它们与条件概率的相互作用。我想找出某人是咖啡饮用者且患有癌症的概率。我应该将和相乘吗?还是应该使用代替?我应该使用哪一个?
如果我们已经确定我们的概率仅适用于患癌症的人群,那么使用 而不是 是不是更有意义?前者更具体,适用于已经建立的条件。因此,我们应该使用,因为 已经是我们联合概率的一部分。这意味着某人患癌症并且是咖啡饮用者的概率为 0.425%:
这种联合概率也适用于另一个方向。我可以通过将 和 相乘来找到某人是咖啡饮用者且患癌症的概率。正如你所看到的,我得到了相同的答案:
如果我们没有任何条件概率可用,那么我们能做的最好就是将 和 相乘,如下所示:
现在想想这个:如果事件 A 对事件 B 没有影响,那么这对条件概率P(B|A)意味着什么?这意味着P(B|A) = P(B),也就是说事件 A 的发生不会影响事件 B 发生的可能性。因此,我们可以更新我们的联合概率公式,无论这两个事件是否相关,都可以是:
最后让我们谈谈并集和条件概率。如果我想计算发生事件 A 或事件 B 的概率,但事件 A 可能会影响事件 B 的概率,我们更新我们的求和规则如下:
作为提醒,这也适用于互斥事件。如果事件 A 和事件 B 不能同时发生,则求和规则P(A|B) × P(B)将得到 0。
二项分布
在本章的其余部分,我们将学习两种概率分布:二项分布和贝塔分布。虽然我们在本书的其余部分不会使用它们,但它们本身是有用的工具,对于理解在一定数量的试验中事件发生的方式至关重要。它们也将是理解我们在第三章中大量使用的概率分布的好起点。让我们探讨一个可能在现实场景中发生的用例。
假设你正在研发一种新的涡轮喷气发动机,并进行了 10 次测试。结果显示有 8 次成功和 2 次失败:
你原本希望获得 90%的成功率,但根据这些数据,你得出结论你的测试只有 80%的成功率。每次测试都耗时且昂贵,因此你决定是时候回到起点重新设计了。
然而,你的一位工程师坚持认为应该进行更多的测试。“唯一确定的方法是进行更多的测试,”她辩称。“如果更多的测试产生了 90%或更高的成功率呢?毕竟,如果你抛硬币 10 次得到 8 次正面,这并不意味着硬币固定在 80%。”
你稍微考虑了工程师的论点,并意识到她有道理。即使是一个公平的硬币翻转也不会总是有相同的结果,尤其是只有 10 次翻转。你最有可能得到五次正面,但也可能得到三、四、六或七次正面。你甚至可能得到 10 次正面,尽管这极不可能。那么,如何确定在底层概率为 90%的情况下 80%成功的可能性?
一个在这里可能相关的工具是二项分布,它衡量了在n次试验中,k次成功可能发生的概率,给定了p的概率。
从视觉上看,二项分布看起来像图 2-1。
在这里,我们看到了在 10 次试验中每个柱子代表的k次成功的概率。这个二项分布假设了一个 90%的概率p,意味着成功发生的概率为 0.90(或 90%)。如果这是真的,那么我们在 10 次试验中获得 8 次成功的概率为 0.1937。在 10 次试验中获得 1 次成功的概率极低,为 0.000000008999,这就是为什么柱子几乎看不见。
我们还可以通过将八次或更少成功的柱子相加来计算八次或更少成功的概率。这将给我们八次或更少成功的概率为 0.2639。
![emds 0201]()
图 2-1. 一个二项分布
那么我们如何实现二项分布呢?我们可以相对容易地从头开始实现(如在附录 A 中分享的),或者我们可以使用像 SciPy 这样的库。示例 2-2 展示了我们如何使用 SciPy 的binom.pmf()
函数(PMF代表“概率质量函数”)来打印从 0 到 10 次成功的二项分布的所有 11 个概率。
示例 2-2. 使用 SciPy 进行二项分布
from scipy.stats import binom
n = 10
p = 0.9
for k in range(n + 1):
probability = binom.pmf(k, n, p)
print("{0} - {1}".format(k, probability))
# OUTPUT:
# 0 - 9.99999999999996e-11
# 1 - 8.999999999999996e-09
# 2 - 3.644999999999996e-07
# 3 - 8.748000000000003e-06
# 4 - 0.0001377809999999999
# 5 - 0.0014880347999999988
# 6 - 0.011160260999999996
# 7 - 0.05739562800000001
# 8 - 0.19371024449999993
# 9 - 0.38742048900000037
# 10 - 0.34867844010000004
正如你所看到的,我们提供n作为试验次数,p作为每次试验成功的概率,k作为我们想要查找概率的成功次数。我们迭代每个成功次数x与我们将看到该数量成功的相应概率。正如我们在输出中看到的,最可能的成功次数是九。
但是,如果我们将八个或更少成功的概率相加,我们会得到 0.2639。这意味着即使基础成功率为 90%,我们也有 26.39%的机会看到八个或更少的成功。所以也许工程师是对的:26.39%的机会并非没有,而且确实可能。
然而,在我们的模型中我们确实做了一个假设,接下来我们将讨论与 beta 分布相关的内容。
从头开始的二项分布
转到附录 A 了解如何在没有 scikit-learn 的情况下从头开始构建二项分布。
Beta 分布
我在使用二项分布的引擎测试模型时做了什么假设?我是否假设了一个参数为真,然后围绕它构建了整个模型?仔细思考并继续阅读。
我的二项分布可能存在问题的地方在于我假设了基础成功率为 90%。这并不是说我的模型毫无价值。我只是展示了如果基础成功率为 90%,那么在 10 次试验中看到 8 个或更少的成功的概率为 26.39%。所以工程师当然没有错,可能存在一个基础成功率为 90%的情况。
但让我们反过来思考这个问题:如果除了 90%之外还有其他基础成功率可以产生 8/10 的成功呢?我们能否用基础成功率 80%?70%?30%?来看到 8/10 的成功?当我们固定 8/10 的成功时,我们能探索概率的概率吗?
与其创建无数个二项分布来回答这个问题,我们可以使用一个工具。Beta 分布允许我们看到在给定alpha成功和beta失败的情况下事件发生的不同基础概率的可能性。
给出八次成功和两次失败的 beta 分布图表如图 2-2 所示。
![emds 0202]()
图 2-2. Beta 分布
Desmos 上的 Beta 分布
如果你想与 beta 分布互动,可以在这里找到 Desmos 图表。
请注意,x 轴表示从 0.0 到 1.0(0%到 100%)的所有基础成功率,y 轴表示在八次成功和两次失败的情况下给定该概率的可能性。换句话说,beta 分布允许我们看到在 8/10 成功的情况下概率的概率。把它看作是一种元概率,所以花点时间理解这个概念!
还要注意,贝塔分布是一个连续函数,这意味着它形成了一个十进制值的连续曲线(与二项分布中整洁且离散的整数相对)。这将使得贝塔分布的数学运算稍微困难,因为 y 轴上的给定密度值不是概率。我们通过曲线下的面积来找到概率。
贝塔分布是一种概率分布,这意味着整个曲线下的面积为 1.0,或者 100%。要找到概率,我们需要找到一个范围内的面积。例如,如果我们想评估 8/10 次成功将产生 90%或更高成功率的概率,我们需要找到 0.9 和 1.0 之间的面积,即 0.225,如图 2-3 中所阴影部分所示。
![emds 0203]()
图 2-3. 90%到 100%之间的面积,为 22.5%
就像我们使用二项分布一样,我们可以使用 SciPy 来实现贝塔分布。每个连续概率分布都有一个累积密度函数(CDF),用于计算给定 x 值之前的面积。假设我想计算到 90%(0.0 到 0.90)的面积,如图 2-4 中所阴影部分所示。
![emds 0204]()
图 2-4. 计算到 90%(0.0 到 0.90)的面积
使用 SciPy 的beta.cdf()
函数非常简单,我需要提供的唯一参数是 x 值,成功次数a和失败次数b,如示例 2-3 所示。
示例 2-3. 使用 SciPy 进行贝塔分布
from scipy.stats import beta
a = 8
b = 2
p = beta.cdf(.90, a, b)
# 0.7748409780000001
print(p)
根据我们的计算,底层成功概率为 90%或更低的概率为 77.48%。
我们如何计算成功概率为 90%或更高的概率,如图 2-5 中所阴影部分所示?
![emds 0205]()
图 2-5. 成功概率为 90%或更高的情况
我们的累积密度函数只计算边界左侧的面积,而不是右侧。想想我们的概率规则,对于概率分布,曲线下的总面积为 1.0。如果我们想找到一个事件的相反概率(大于 0.90 而不是小于 0.90),只需从 1.0 中减去小于 0.90 的概率,剩下的概率将捕获大于 0.90 的情况。图 2-6 说明了我们如何进行这种减法。
![emds 0206]()
图 2-6. 找到成功概率大于 90%的概率
示例 2-4 展示了我们如何在 Python 中计算这种减法操作。
示例 2-4. 在贝塔分布中进行减法以获得正确的面积
from scipy.stats import beta
a = 8
b = 2
p = 1.0 - beta.cdf(.90, a, b)
# 0.22515902199999993
print(p)
这意味着在 8/10 次成功的引擎测试中,成功率为 90% 或更高的概率只有 22.5%。但成功率低于 90% 的概率约为 77.5%。我们在这里的胜算不高,但如果我们感觉幸运,我们可以通过更多的测试来赌那 22.5% 的机会。如果我们的首席财务官为 26 次额外测试提供资金,结果是 30 次成功和 6 次失败,我们的贝塔分布将如 图 2-7 所示。
![emds 0207]()
图 2-7. 30 次成功和 6 次失败后的贝塔分布
注意我们的分布变得更窄,因此更有信心地认为成功率在一个更小的范围内。不幸的是,我们达到 90% 成功率最低的概率已经减少,从 22.5% 降至 13.16%,如 例子 2-5 所示。
例子 2-5. 具有更多试验的贝塔分布
from scipy.stats import beta
a = 30
b = 6
p = 1.0 - beta.cdf(.90, a, b)
# 0.13163577484183708
print(p)
此时,最好离开并停止做测试,除非你想继续赌博对抗那 13.16% 的机会,并希望峰值向右移动。
最后,我们如何计算中间的面积?如果我想找到成功率在 80% 和 90% 之间的概率,如 图 2-8 所示,该怎么办?
![emds 0208]()
图 2-8. 成功率在 80% 和 90% 之间的概率
仔细考虑一下你可能如何处理这个问题。如果我们像在 图 2-9 中那样,从 .80 后面的面积中减去 .90 后面的面积会怎样?
![emds 0209]()
图 2-9. 获取 .80 和 .90 之间的面积
这会给我们提供 .80 和 .90 之间的面积吗?是的,它会,并且会产生一个 .3386 或 33.86% 的概率。这是我们如何在 Python 中计算它的方法(例子 2-6)。
例子 2-6. 使用 SciPy 计算贝塔分布中间区域
from scipy.stats import beta
a = 8
b = 2
p = beta.cdf(.90, a, b) - beta.cdf(.80, a, b)
# 0.33863336200000016
print(p)
贝塔分布是一种迷人的工具,用于根据有限的观察来衡量事件发生与不发生的概率。它使我们能够推理关于概率的概率,并且我们可以在获得新数据时更新它。我们也可以将其用于假设检验,但我们将更加强调使用正态分布和 T 分布来进行这种目的,如 第三章 中所述。
从头开始的贝塔分布
要了解如何从头开始实现贝塔分布,请参考 附录 A。
结论
在本章中,我们涵盖了很多内容!我们不仅讨论了概率的基础知识、逻辑运算符和贝叶斯定理,还介绍了概率分布,包括二项式分布和贝塔分布。在下一章中,我们将讨论更著名的分布之一,正态分布,以及它与假设检验的关系。
如果你想了解更多关于贝叶斯概率和统计的知识,一本很棒的书是Bayesian Statistics the Fun Way,作者是 Will Kurt(No Starch Press)。还可以在 O’Reilly 平台上找到互动的Katacoda 场景。
练习
-
今天有 30%的降雨概率,你的雨伞订单准时到达的概率为 40%。你渴望今天在雨中散步,但没有雨伞你无法做到!
下雨的概率和你的雨伞到达的概率是多少?
-
今天有 30%的降雨概率,你的雨伞订单准时到达的概率为 40%。
只有在不下雨或者你的雨伞到达时,你才能出门办事。
不下雨或者你的雨伞到达的概率是多少?
-
今天有 30%的降雨概率,你的雨伞订单准时到达的概率为 40%。
然而,你发现如果下雨,你的雨伞准时到达的概率只有 20%。
下雨的概率和你的雨伞准时到达的概率是多少?
-
你从拉斯维加斯飞往达拉斯的航班上有 137 名乘客预订了座位。然而,这是拉斯维加斯的一个星期天早上,你估计每位乘客有 40%的可能性不会出现。
你正在努力计算要超售多少座位,以免飞机空飞。
至少有 50 名乘客不会出现的概率有多大?
-
你抛了一枚硬币 19 次,其中 15 次是正面,4 次是反面。
你认为这枚硬币有可能是公平的吗?为什么?
答案在附录 B 中。
第三章:描述性和推断性统计
统计学是收集和分析数据以发现有用发现或预测导致这些发现发生的原因的实践。概率在统计学中经常起着重要作用,因为我们使用数据来估计事件发生的可能性。
统计学可能并不总是受到赞誉,但它是许多数据驱动创新的核心。机器学习本身就是一种统计工具,寻找可能的假设来相关不同变量之间的关系。然而,即使对于专业统计学家来说,统计学也存在许多盲点。我们很容易陷入数据所说的内容,而忘记了要问数据来自何处。随着大数据、数据挖掘和机器学习加速推动统计算法的自动化,这些问题变得更加重要。因此,拥有扎实的统计学和假设检验基础非常重要,这样你就不会把这些自动化处理当作黑匣子。
在本节中,我们将介绍统计学和假设检验的基础知识。从描述性统计开始,我们将学习总结数据的常见方法。之后,我们将进入推断统计,试图根据样本揭示总体属性。
什么是数据?
定义“数据”可能看起来有些奇怪,因为我们都使用并认为理所当然。但我认为有必要这样做。如果你问任何人数据是什么,他们可能会回答类似于“你知道的...数据!就是...你知道的...信息!”而不会深入探讨。现在它似乎被宣传为至关重要的事物。不仅是真相的来源...还有智慧!这是人工智能的燃料,人们相信你拥有的数据越多,你就拥有的真相就越多。因此,你永远不可能拥有足够的数据。它将揭示重新定义你的业务策略所需的秘密,甚至可能创造人工通用智能。但让我提供一个关于数据是什么的实用观点。数据本身并不重要。数据的分析(以及它是如何产生的)是所有这些创新和解决方案的驱动力。
想象一下,如果你拿到一张一个家庭的照片。你能根据这张照片揭示这个家庭的故事吗?如果你有 20 张照片呢?200 张照片?2,000 张照片?你需要多少张照片才能了解他们的故事?你需要他们在不同情况下的照片吗?一个人和一起的照片?和亲戚朋友在一起的照片?在家里和工作中的照片?
数据就像照片一样;它提供了故事的快照。连续的现实和背景并没有完全捕捉到,也没有驱动这个故事的无限数量的变量。正如我们将讨论的,数据可能存在偏见。它可能存在缺口,缺少相关变量。理想情况下,我们希望有无限量的数据捕捉无限数量的变量,有如此之多的细节,我们几乎可以重新创造现实并构建替代现实!但这可能吗?目前,不可能。即使将全球最强大的超级计算机组合在一起,也无法接近以数据形式捕捉整个世界的全部内容。
因此,我们必须缩小范围,使我们的目标变得可行。父亲打高尔夫球的几张战略照片可以很容易地告诉我们他是否擅长高尔夫。但是仅凭照片来解读他的整个人生故事?那可能是不可能的。有很多东西是无法在快照中捕捉到的。这些实际问题也应该在处理数据项目时应用,因为数据实际上只是捕捉特定时间的快照,只捕捉到了它所针对的内容(就像相机一样)。我们需要保持我们的目标集中,这有助于收集相关和完整的数据。如果我们的目标过于宽泛和开放,我们可能会遇到虚假发现和不完整数据集的问题。这种实践,被称为数据挖掘,有其时机和地点,但必须小心进行。我们将在本章末重新讨论这个问题。
即使目标明确,我们仍然可能在数据方面遇到问题。让我们回到确定几张战略照片是否能告诉我们父亲是否擅长高尔夫的问题。也许如果你有一张他挥杆中的照片,你就能看出他的动作是否正确。或者也许如果你看到他在一个洞位上欢呼和击掌,你可以推断他得了一个好成绩。也许你只需拍一张他的记分卡的照片!但重要的是要注意所有这些情况都可能是伪造的或脱离了上下文。也许他在为别人欢呼,或者记分卡不是他的,甚至是伪造的。就像这些照片一样,数据并不能捕捉到背景或解释。这是一个非常重要的观点,因为数据提供线索,而不是真相。这些线索可以引导我们找到真相,也可能误导我们得出错误的结论。
这就是为什么对数据来源感到好奇是如此重要的技能。询问数据是如何创建的,由谁创建的,以及数据未捕捉到什么。很容易陷入数据所说的内容而忘记询问数据来自何处。更糟糕的是,有广泛的观点认为可以将数据填入机器学习算法中,并期望计算机解决所有问题。但正如谚语所说,“垃圾进,垃圾出”。难怪根据 VentureBeat 的数据,只有13%的机器学习项目成功。成功的机器学习项目对数据进行了思考和分析,以及产生数据的过程。
描述性统计与推断性统计
当你听到“统计学”这个词时,你会想到什么?是计算平均值、中位数、众数、图表、钟形曲线和其他用于描述数据的工具吗?这是统计学最常被理解的部分,称为描述性统计,我们用它来总结数据。毕竟,浏览百万条数据记录还是对其进行总结更有意义?我们将首先涵盖统计学的这一领域。
推断性统计试图揭示关于更大总体的属性,通常基于样本。它经常被误解,比描述性统计更难理解。通常我们对研究一个太大以至无法观察的群体感兴趣(例如,北美地区青少年的平均身高),我们必须仅使用该群体的少数成员来推断关于他们的结论。正如你所猜测的那样,这并不容易做到。毕竟,我们试图用一个可能不具代表性的样本来代表一个总体。我们将在探讨过程中探讨这些警告。
总体、样本和偏差
在我们深入研究描述性和推断性统计之前,将一些定义和与之相关的具体示例联系起来可能是一个好主意。
总体是我们想要研究的特定感兴趣的群体,比如“北美地区所有 65 岁以上的老年人”,“苏格兰所有金毛猎犬”,或者“洛斯阿尔托斯高中目前的高中二年级学生”。请注意我们对定义总体的边界。有些边界很宽泛,涵盖了广阔地理区域或年龄组的大群体。其他则非常具体和小,比如洛斯阿尔托斯高中的高中二年级学生。如何确定总体的定义取决于你对研究感兴趣的内容。
样本是总体的一个理想随机和无偏子集,我们用它来推断总体的属性。我们经常不得不研究样本,因为调查整个总体并不总是可能的。当然,如果人口小且易接触,那么获取一些人口是更容易的。但是测量北美地区所有 65 岁以上的老年人?那不太可能实际可行!
总体可以是抽象的!
需要注意的是,人口可以是理论的,而不是实际可触及的。在这些情况下,我们的人口更像是从抽象事物中抽取的样本。这里是我最喜欢的例子:我们对一个机场在下午 2 点到 3 点之间起飞的航班感兴趣,但在那个时间段内的航班数量不足以可靠地预测这些航班的延误情况。因此,我们可能将这个人口视为从所有在下午 2 点到 3 点之间起飞的理论航班中抽取的样本。
这类问题是为什么许多研究人员借助模拟生成数据的原因。模拟可能是有用的,但很少是准确的,因为模拟只捕捉了有限的变量,并且内置了假设。
如果我们要根据样本推断人口的属性,那么样本尽可能随机是很重要的,这样我们才不会扭曲我们的结论。这里举个例子。假设我是亚利桑那州立大学的一名大学生。我想找出美国大学生每周观看电视的平均小时数。我走出宿舍,开始随机询问路过的学生,几个小时后完成了数据收集。问题在哪里?
问题在于我们的学生样本可能存在偏见,这意味着它通过过度代表某一群体而牺牲其他群体来扭曲我们的发现。我的研究将人口定义为“美国的大学生”,而不是“亚利桑那州立大学的大学生”。我只是在一个特定大学对学生进行调查,代表整个美国的所有大学生!这真的公平吗?
不太可能全国各地的大学都具有相同的学生属性。如果亚利桑那州立大学的学生比其他大学的学生观看电视时间更长怎么办?使用他们来代表整个国家难道不会扭曲结果吗?也许这是可能的,因为在亚利桑那州的坦佩市通常太热了,所以看电视是一种常见的消遣(据说,我会知道;我在凤凰城住了很多年)。其他气候较温和的地方的大学生可能会进行更多的户外活动,看电视时间更少。
这只是一个可能的变量,说明用一个大学的学生样本来代表整个美国的大学生是一个不好的主意。理想情况下,我应该随机调查全国各地不同大学的大学生。这样我就有了更具代表性的样本。
然而,偏见并不总是地理性的。假设我竭尽全力在全美各地调查学生。我策划了一个社交媒体活动,在 Twitter 和 Facebook 上让各大学分享调查,这样他们的学生就会看到并填写。我收到了数百份关于全国学生电视习惯的回复,觉得我已经征服了偏见的野兽...或者我真的做到了吗?
如果那些足够在社交媒体上看到投票的学生也更有可能看更多电视呢?如果他们在社交媒体上花很多时间,他们可能不介意娱乐性的屏幕时间。很容易想象他们已经准备好在另一个标签上流媒体 Netflix 和 Hulu!这种特定类型的偏见,即特定群体更有可能加入样本的偏见,被称为自我选择偏差。
糟糕!你就是赢不了,是吗?如果你足够长时间地考虑,数据偏差似乎是不可避免的!而且通常确实如此。许多混杂变量,或者我们没有考虑到的因素,都会影响我们的研究。数据偏差问题昂贵且难以克服,而机器学习尤其容易受到影响。
克服这个问题的方法是真正随机地从整个人口中选择学生,他们不能自愿地加入或退出样本。这是减轻偏见的最有效方法,但正如你所想象的那样,这需要大量协调的资源。
好了,关于人口、样本和偏差的讨论就到此为止。让我们继续进行一些数学和描述性统计。只要记住,数学和计算机不会意识到数据中的偏差。这取决于你作为一名优秀的数据科学专业人员来检测!始终要问关于数据获取方式的问题,然后仔细审查该过程可能如何使数据产生偏差。
机器学习中的样本和偏差
这些与抽样和偏差有关的问题也延伸到机器学习领域。无论是线性回归、逻辑回归还是神经网络,都会使用数据样本来推断预测。如果数据存在偏差,那么它将引导机器学习算法做出有偏见的结论。
这方面有许多记录的案例。刑事司法一直是机器学习的一个棘手应用,因为它一再显示出在每个意义上都存在偏见,由于数据集中存在少数族群,导致对少数族群进行歧视。2017 年,沃尔沃测试了训练过捕捉鹿、麋鹿和驯鹿数据集的自动驾驶汽车。然而,它在澳大利亚没有驾驶数据,因此无法识别袋鼠,更不用说理解它们的跳跃动作了!这两个都是有偏见数据的例子。
描述性统计
描述性统计是大多数人熟悉的领域。我们将介绍一些基础知识,如均值、中位数和众数,然后是方差、标准差和正态分布。
均值和加权均值
均值是一组值的平均值。这个操作很简单:将值相加,然后除以值的数量。均值很有用,因为它显示了观察到的一组值的“重心”在哪里。
均值的计算方式对于人口和样本是相同的。示例 3-1 展示了八个值的样本以及如何在 Python 中计算它们的均值。
示例 3-1. 在 Python 中计算均值
# Number of pets each person owns
sample = [1, 3, 2, 5, 7, 0, 2, 3]
mean = sum(sample) / len(sample)
print(mean) # prints 2.875
正如你所看到的,我们对八个人关于他们拥有的宠物数量进行了调查。样本的总和为 23,样本中的项目数为 8,因此这给我们一个均值为 2.875,因为 23/8 = 2.875。
你将看到两个版本的均值:样本均值和总体均值如此表达:
回想一下求和符号表示将所有项目相加。n和N分别代表样本和总体大小,但在数学上它们表示相同的东西:项目的数量。对于称为样本均值(“x-bar”)和总体均值(“mu”)的命名也是如此。无论是还是都是相同的计算,只是根据我们处理的是样本还是总体而有不同的名称。
均值可能对你来说很熟悉,但关于均值有一些不太为人知的东西:均值实际上是一种称为加权均值的加权平均值。我们通常使用的均值给予每个值相同的重要性。但我们可以操纵均值,给每个项目赋予不同的权重:
当我们希望某些值对均值的贡献大于其他值时,这将会很有帮助。一个常见的例子是对学术考试进行加权以得出最终成绩。如果你有三次考试和一次期末考试,我们给予每次考试 20%的权重,期末考试 40%的权重,我们如何表达它在示例 3-2 中。
示例 3-2. 在 Python 中计算加权均值
# Three exams of .20 weight each and final exam of .40 weight
sample = [90, 80, 63, 87]
weights = [.20, .20, .20, .40]
weighted_mean = sum(s * w for s,w in zip(sample, weights)) / sum(weights)
print(weighted_mean) # prints 81.4
我们通过相应的乘法对每个考试分数进行加权,而不是通过值计数进行除法,而是通过权重总和进行除法。权重不必是百分比,因为用于权重的任何数字最终都将被比例化。在示例 3-3 中,我们对每个考试赋予“1”的权重,但对最终考试赋予“2”的权重,使其比考试的权重大一倍。我们仍然会得到相同的答案 81.4,因为这些值仍然会被比例化。
示例 3-3. 在 Python 中计算加权平均值
# Three exams of .20 weight each and final exam of .40 weight
sample = [90, 80, 63, 87]
weights = [1.0, 1.0, 1.0, 2.0]
weighted_mean = sum(s * w for s,w in zip(sample, weights)) / sum(weights)
print(weighted_mean) # prints 81.4
中位数
中位数是一组有序值中间最靠近的值。您按顺序排列值,中位数将是中间的值。如果您有偶数个值,您将平均两个中间值。我们可以看到在示例 3-4 中,我们样本中拥有的宠物数量的中位数为 7:
0, 1, 5, *7*, 9, 10, 14
示例 3-4. 在 Python 中计算中位数
# Number of pets each person owns
sample = [0, 1, 5, 7, 9, 10, 14]
def median(values):
ordered = sorted(values)
print(ordered)
n = len(ordered)
mid = int(n / 2) - 1 if n % 2 == 0 else int(n/2)
if n % 2 == 0:
return (ordered[mid] + ordered[mid+1]) / 2.0
else:
return ordered[mid]
print(median(sample)) # prints 7
当数据被异常值或与其他值相比极端大和小的值所扭曲时,中位数可以成为均值的有用替代。这里有一个有趣的轶事来理解为什么。1986 年,北卡罗来纳大学教堂山分校地理学毕业生的年起薪平均值为$250,000。其他大学平均为$22,000。哇,UNC-CH 一定有一个了不起的地理学项目!
但实际上,为什么北卡罗来纳大学的地理学项目如此赚钱呢?嗯...迈克尔·乔丹是他们的毕业生之一。确实,有史以来最著名的 NBA 球员之一,确实是从北卡罗来纳大学获得地理学学位的。然而,他开始他的职业生涯是打篮球,而不是学习地图。显然,这是一个造成了巨大异常值的混杂变量,它严重扭曲了收入平均值。
这就是为什么在受异常值影响较大的情况下(例如与收入相关的数据)中,中位数可能比均值更可取。它对异常值不太敏感,并严格根据它们的相对顺序将数据划分到中间,而不是准确地落在数字线上。当您的中位数与均值非常不同时,这意味着您的数据集具有异常值。
众数
众数是出现频率最高的一组值。当您的数据重复时,想要找出哪些值出现最频繁时,它就变得很有用。
当没有任何值出现超过一次时,就没有众数。当两个值出现的频率相等时,数据集被认为是双峰的。在示例 3-5 中,我们计算了我们的宠物数据集的众数,果然我们看到这是双峰的,因为 2 和 3 都出现最多(并且相等)。
示例 3-5. 在 Python 中计算众数
# Number of pets each person owns
from collections import defaultdict
sample = [1, 3, 2, 5, 7, 0, 2, 3]
def mode(values):
counts = defaultdict(lambda: 0)
for s in values:
counts[s] += 1
max_count = max(counts.values())
modes = [v for v in set(values) if counts[v] == max_count]
return modes
print(mode(sample)) # [2, 3]
在实际应用中,除非您的数据重复,否则众数不经常使用。这通常在整数、类别和其他离散变量中遇到。
方差和标准差
当我们开始讨论方差和标准差时,这就变得有趣起来了。方差和标准差让人们感到困惑的一点是,对于样本和总体有一些计算差异。我们将尽力清楚地解释这些差异。
总体方差和标准差
在描述数据时,我们经常对测量平均值和每个数据点之间的差异感兴趣。这让我们对数据的“分布”有了一种感觉。
假设我对研究我的工作人员拥有的宠物数量感兴趣(请注意,我将这定义为我的总体,而不是样本)。我的工作人员有七个人。
我取得他们拥有的所有宠物数量的平均值,我得到 6.571。让我们从每个值中减去这个平均值。这将展示给我们每个值距离平均值有多远,如表 3-1 所示。
表 3-1. 我的员工拥有的宠物数量
Value |
Mean |
Difference |
0 |
6.571 |
-6.571 |
1 |
6.571 |
-5.571 |
5 |
6.571 |
-1.571 |
7 |
6.571 |
0.429 |
9 |
6.571 |
2.429 |
10 |
6.571 |
3.429 |
14 |
6.571 |
7.429 |
让我们在数轴上用“X”表示平均值来可视化这一点,见图 3-1。
![emds 0301]()
图 3-1. 可视化我们数据的分布,其中“X”是平均值
嗯...现在考虑为什么这些信息有用。这些差异让我们了解数据的分布情况以及值距离平均值有多远。有没有一种方法可以将这些差异合并成一个数字,快速描述数据的分布情况?
你可能会想要取这些差异的平均值,但当它们相加时,负数和正数会互相抵消。我们可以对绝对值求和(去掉负号并使所有值为正)。一个更好的方法是在求和之前对这些差异进行平方。这不仅消除了负值(因为平方一个负数会使其变为正数),而且放大了较大的差异,并且在数学上更容易处理(绝对值的导数不直观)。之后,对平方差异取平均值。这将给我们方差,一个衡量数据分布范围的指标。
这里是一个数学公式,展示如何计算方差:
更正式地,这是总体的方差:
在 Python 中计算我们宠物示例的总体方差在示例 3-6 中展示。
示例 3-6. 在 Python 中计算方差
# Number of pets each person owns
data = [0, 1, 5, 7, 9, 10, 14]
def variance(values):
mean = sum(values) / len(values)
_variance = sum((v - mean) ** 2 for v in values) / len(values)
return _variance
print(variance(data)) # prints 21.387755102040813
因此,我的办公室员工拥有的宠物数量的方差为 21.387755。好的,但这到底意味着什么?可以合理地得出结论,更高的方差意味着更广泛的分布,但我们如何将其与我们的数据联系起来?这个数字比我们的任何观察结果都要大,因为我们进行了大量的平方和求和,将其放在完全不同的度量标准上。那么我们如何将其压缩回原来的尺度?
平方的反面是平方根,所以让我们取方差的平方根,得到标准差。这是将方差缩放为以“宠物数量”表示的数字,这使得它更有意义:
要在 Python 中实现,我们可以重用variance()
函数并对其结果进行sqrt()
。现在我们有了一个std_dev()
函数,如示例 3-7 所示。
示例 3-7. 在 Python 中计算标准差
from math import sqrt
# Number of pets each person owns
data = [0, 1, 5, 7, 9, 10, 14]
def variance(values):
mean = sum(values) / len(values)
_variance = sum((v - mean) ** 2 for v in values) / len(values)
return _variance
def std_dev(values):
return sqrt(variance(values))
print(std_dev(data)) # prints 4.624689730353898
运行示例 3-7 中的代码,你会看到我们的标准差约为 4.62 只宠物。因此,我们可以用我们开始的尺度来表示我们的扩散,这使得我们的方差更容易解释。我们将在第五章中看到标准差的一些重要应用。
为什么是平方?
关于方差,如果中的指数让你感到困扰,那是因为它提示你对其进行平方根运算以获得标准差。这是一个小小的提醒,告诉你正在处理需要进行平方根运算的平方值。
样本方差和标准差
在前一节中,我们讨论了总体的方差和标准差。然而,当我们为样本计算这两个公式时,我们需要应用一个重要的调整:
你注意到了区别吗?当我们对平方差值求平均时,我们除以n-1 而不是总项目数n。为什么要这样做?我们这样做是为了减少样本中的任何偏差,不低估基于我们的样本的总体方差。通过在除数中计算少一个项目的值,我们增加了方差,因此捕捉到了样本中更大的不确定性。
如果我们的宠物数据是一个样本,而不是一个总体,我们应相应地进行调整。在示例 3-8 中,我修改了之前的variance()
和std_dev()
Python 代码,可选择提供一个参数is_sample
,如果为True
,则从方差的除数中减去 1。
示例 3-8. 计算样本的标准差
from math import sqrt
# Number of pets each person owns
data = [0, 1, 5, 7, 9, 10, 14]
def variance(values, is_sample: bool = False):
mean = sum(values) / len(values)
_variance = sum((v - mean) ** 2 for v in values) /
(len(values) - (1 if is_sample else 0))
return _variance
def std_dev(values, is_sample: bool = False):
return sqrt(variance(values, is_sample))
print("VARIANCE = {}".format(variance(data, is_sample=True))) # 24.95238095238095
print("STD DEV = {}".format(std_dev(data, is_sample=True))) # 4.99523582550223
请注意,在示例 3-8 中,我的方差和标准差与之前将它们视为总体而不是样本的示例相比有所增加。回想一下,在示例 3-7 中,将其视为总体时,标准差约为 4.62。但是在这里将其视为样本(通过从方差分母中减去 1),我们得到的值约为 4.99。这是正确的,因为样本可能存在偏差,不完全代表总体。因此,我们增加方差(从而增加标准差)以增加对值分布范围的估计。较大的方差/标准差显示出较大范围的不确定性。
就像均值(样本为 ,总体为 ),你经常会看到一些符号表示方差和标准差。样本和均值的标准差分别由s和指定。这里再次是样本和总体标准差的公式:
方差将是这两个公式的平方,撤销平方根。因此,样本和总体的方差分别为和:
再次,平方有助于暗示应该取平方根以获得标准差。
正态分布
在上一章中我们提到了概率分布,特别是二项分布和贝塔分布。然而,最著名的分布是正态分布。正态分布,也称为高斯分布,是一种对称的钟形分布,大部分质量集中在均值附近,其传播程度由标准差定义。两侧的“尾巴”随着远离均值而变得更细。
图 3-2 是金毛犬体重的正态分布。请注意,大部分质量集中在 64.43 磅的均值附近。
![emds 0302]()
图 3-2. 一个正态分布
发现正态分布
正态分布在自然界、工程、科学和其他领域中经常出现。我们如何发现它呢?假设我们对 50 只成年金毛犬的体重进行抽样,并将它们绘制在数轴上,如图 3-3 所示。
![emds 0303]()
图 3-3. 50 只金毛犬体重的样本
请注意,我们在中心附近有更多的值,但随着向左或向右移动,我们看到的值越来越少。根据我们的样本,看起来我们很不可能看到体重为 57 或 71 磅的金毛犬。但体重为 64 或 65 磅?是的,这显然很可能。
有没有更好的方法来可视化这种可能性,以查看我们更有可能从人群中抽样到哪些金毛犬的体重?我们可以尝试创建一个直方图,它根据相等长度的数值范围将值分组(或“箱”),然后使用柱状图显示每个范围内的值的数量。在图 3-4 中,我们创建了一个将值分组为 0.5 磅范围的直方图。
这个直方图并没有显示出数据的任何有意义的形状。原因是因为我们的箱太小了。我们没有足够大或无限的数据量来有意义地在每个箱中有足够的点。因此,我们必须使我们的箱更大。让我们使每个箱的长度为三磅,如图 3-5 所示。
![emds 0304]()
图 3-4. 一张金毛犬体重的直方图
![emds 0305]()
图 3-5. 一个更有意义的直方图
现在我们有了进展!正如你所看到的,如果我们将箱的大小调整得恰到好处(在本例中,每个箱的范围为三磅),我们开始在数据中看到一个有意义的钟形。这不是一个完美的钟形,因为我们的样本永远不会完全代表人群,但这很可能是我们的样本遵循正态分布的证据。如果我们使用适当的箱大小拟合一个直方图,并将其缩放为面积为 1.0(这是概率分布所需的),我们会看到一个粗略的钟形曲线代表我们的样本。让我们在图 3-6 中将其与我们的原始数据点一起展示。
看着这个钟形曲线,我们可以合理地期望金毛寻回犬的体重最有可能在 64.43(均值)左右,但在 55 或 73 时不太可能。比这更极端的情况变得非常不太可能。
![emds 0306]()
图 3-6。拟合到数据点的正态分布
正态分布的性质
正态分布具有几个重要的属性,使其有用:
概率密度函数(PDF)
标准偏差在正态分布中起着重要作用,因为它定义了它的“扩散程度”。实际上,它是与均值一起的参数之一。创建正态分布的概率密度函数(PDF)如下:
哇,这真是一个冗长的话题,不是吗?我们甚至在第一章中看到我们的朋友欧拉数和一些疯狂的指数。这是我们如何在 Python 中表达它的示例 3-9。
示例 3-9。Python 中的正态分布函数
# normal distribution, returns likelihood
def normal_pdf(x: float, mean: float, std_dev: float) -> float:
return (1.0 / (2.0 * math.pi * std_dev ** 2) ** 0.5) *
math.exp(-1.0 * ((x - mean) ** 2 / (2.0 * std_dev ** 2)))
这个公式中有很多要解开的内容,但重要的是它接受均值和标准偏差作为参数,以及一个 x 值,这样你就可以查找给定值的可能性。
就像第二章中的贝塔分布一样,正态分布是连续的。这意味着要获取一个概率,我们需要积分一系列x值以找到一个区域。
实际上,我们将使用 SciPy 来为我们进行这些计算。
累积分布函数(CDF)
对于正态分布,垂直轴不是概率,而是数据的可能性。要找到概率,我们需要查看一个给定范围,然后找到该范围下曲线的面积。假设我想找到金毛寻回犬体重在 62 到 66 磅之间的概率。图 3-7 显示了我们要找到面积的范围。
![emds 0307]()
图 3-7。测量 62 到 66 磅之间概率的 CDF
我们已经在第二章中用贝塔分布完成了这个任务,就像贝塔分布一样,这里也有一个累积密度函数(CDF)。让我们遵循这种方法。
正如我们在上一章中学到的,CDF 提供了给定分布的给定 x 值的面积直到。让我们看看我们的金毛寻回犬正态分布的 CDF 是什么样子,并将其与 PDF 放在图 3-8 中进行参考。
![emds 0308]()
图 3-8. PDF 与其 CDF 并排
注意两个图之间的关系。CDF 是一个 S 形曲线(称为 Sigmoid 曲线),它在 PDF 中投影到该范围的面积。观察图 3-9,当我们捕捉从负无穷到 64.43(均值)的面积时,我们的 CDF 显示的值恰好为.5 或 50%!
![emds 0309]()
图 3-9. 金毛寻回犬体重的 PDF 和 CDF,测量均值的概率
这个区域的面积为 0.5 或 50%,这是由于我们正态分布的对称性,我们可以期望钟形曲线的另一侧也有 50%的面积。
要在 Python 中使用 SciPy 计算到 64.43 的面积,使用norm.cdf()
函数,如示例 3-10 所示。
示例 3-10. Python 中的正态分布 CDF
from scipy.stats import norm
mean = 64.43
std_dev = 2.99
x = norm.cdf(64.43, mean, std_dev)
print(x) # prints 0.5
就像我们在第二章中所做的那样,我们可以通过减去面积来演绎地找到中间范围的面积。如果我们想找到 62 到 66 磅之间观察到金毛寻回犬的概率,我们将计算到 66 的面积并减去到 62 的面积,如图 3-10 所示。
![emds 0310]()
图 3-10. 寻找中间范围的概率
在 Python 中使用 SciPy 进行这个操作就像在示例 3-11 中所示的两个 CDF 操作相减一样简单。
示例 3-11. 使用 CDF 获取中间范围概率
from scipy.stats import norm
mean = 64.43
std_dev = 2.99
x = norm.cdf(66, mean, std_dev) - norm.cdf(62, mean, std_dev)
print(x) # prints 0.4920450147062894
你应该发现在 62 到 66 磅之间观察到金毛寻回犬的概率为 0.4920,或者大约为 49.2%。
逆 CDF
当我们在本章后面开始进行假设检验时,我们会遇到需要查找 CDF 上的一个区域然后返回相应 x 值的情况。当然,这是对 CDF 的反向使用,所以我们需要使用逆 CDF,它会翻转轴,如图 3-11 所示。
这样,我们现在可以查找一个概率然后返回相应的 x 值,在 SciPy 中我们会使用norm.ppf()
函数。例如,我想找到 95%的金毛寻回犬体重在以下。当我使用示例 3-12 中的逆 CDF 时,这很容易做到。
![emds 0311]()
图 3-11. 逆 CDF,也称为 PPF 或分位数函数
示例 3-12. 在 Python 中使用逆 CDF(称为ppf()
)
from scipy.stats import norm
x = norm.ppf(.95, loc=64.43, scale=2.99)
print(x) # 69.3481123445849
我发现 95%的金毛寻回犬体重为 69.348 磅或更少。
你也可以使用逆 CDF 生成遵循正态分布的随机数。如果我想创建一个模拟,生成一千个真实的金毛寻回犬体重,我只需生成一个介于 0.0 和 1.0 之间的随机值,传递给逆 CDF,然后返回体重值,如示例 3-13 所示。
示例 3-13. 从正态分布生成随机数
import random
from scipy.stats import norm
for i in range(0,1000):
random_p = random.uniform(0.0, 1.0)
random_weight = norm.ppf(random_p, loc=64.43, scale=2.99)
print(random_weight)
当然,NumPy 和其他库可以为您生成分布的随机值,但这突出了逆 CDF 很方便的一个用例。
从头开始实现 CDF 和逆 CDF
要学习如何在 Python 中从头开始实现 CDF 和逆 CDF,请参考附录 A。
Z 分数
将正态分布重新缩放,使均值为 0,标准差为 1 是很常见的,这被称为标准正态分布。这样可以轻松比较一个正态分布的扩展与另一个正态分布,即使它们具有不同的均值和方差。
特别重要的是,标准正态分布将所有 x 值表示为标准差,称为Z 分数。将 x 值转换为 Z 分数使用基本的缩放公式:
这是一个例子。我们有两个来自两个不同社区的房屋。社区 A 的平均房屋价值为$140,000,标准差为$3,000。社区 B 的平均房屋价值为$800,000,标准差为$10,000。
现在我们有两个来自每个社区的房屋。社区 A 的房屋价值为$150,000,社区 B 的房屋价值为$815,000。哪个房屋相对于其社区的平均房屋更昂贵?
如果我们将这两个值用标准差表示,我们可以相对于每个社区的均值进行比较。使用 Z 分数公式:
因此,A 社区的房屋实际上相对于其社区要昂贵得多,因为它们的 Z 分数分别为和 1.5。
下面是我们如何将来自具有均值和标准差的给定分布的 x 值转换为 Z 分数,反之亦然,如示例 3-14 所示。
示例 3-14. 将 Z 分数转换为 x 值,反之亦然
def z_score(x, mean, std):
return (x - mean) / std
def z_to_x(z, mean, std):
return (z * std) + mean
mean = 140000
std_dev = 3000
x = 150000
# Convert to Z-score and then back to X
z = z_score(x, mean, std_dev)
back_to_x = z_to_x(z, mean, std_dev)
print("Z-Score: {}".format(z)) # Z-Score: 3.333
print("Back to X: {}".format(back_to_x)) # Back to X: 150000.0
z_score()
函数将接受一个 x 值,并根据均值和标准差将其按标准偏差进行缩放。z_to_x()
函数将接受一个 Z 分数,并将其转换回 x 值。研究这两个函数,你可以看到它们的代数关系,一个解决 Z 分数,另一个解决 x 值。然后,我们将一个 x 值为 8.0 转换为的 Z 分数,然后将该 Z 分数转换回 x 值。
推断统计学
我们迄今为止所涵盖的描述性统计学是常见的。然而,当我们涉及推断统计学时,样本和总体之间的抽象关系得以充分发挥。这些抽象的微妙之处不是你想要匆忙通过的,而是要花时间仔细吸收。正如前面所述,我们作为人类有偏见,很快就会得出结论。成为一名优秀的数据科学专业人员需要你抑制那种原始的欲望,并考虑其他解释可能存在的可能性。推测没有任何解释,某一发现只是巧合和随机的也是可以接受的(也许甚至是开明的)。
首先让我们从为所有推断统计学奠定基础的定理开始。
中心极限定理
正态分布之所以有用的原因之一是因为它在自然界中经常出现,比如成年金毛犬的体重。然而,在自然人口之外的更迷人的背景下,它也会出现。当我们开始从人口中抽取足够大的样本时,即使该人口不遵循正态分布,正态分布仍然会出现。
假设我正在测量一个真正随机且均匀的人口。0.0 到 1.0 之间的任何值都是同等可能的,没有任何值有任何偏好。但当我们从这个人口中抽取越来越大的样本,计算每个样本的平均值,然后将它们绘制成直方图时,会发生一些有趣的事情。在示例 3-15 中运行这段 Python 代码,并观察图 3-12 中的图表。
示例 3-15. 在 Python 中探索中心极限定理
# Samples of the uniform distribution will average out to a normal distribution.
import random
import plotly.express as px
sample_size = 31
sample_count = 1000
# Central limit theorem, 1000 samples each with 31
# random numbers between 0.0 and 1.0
x_values = [(sum([random.uniform(0.0, 1.0) for i in range(sample_size)]) / \
sample_size)
for _ in range(sample_count)]
y_values = [1 for _ in range(sample_count)]
px.histogram(x=x_values, y = y_values, nbins=20).show()
![emds 0312]()
图 3-12. 取样本均值(每个大小为 31)并绘制它们
等等,当作为 31 个一组的均匀随机数取样然后求平均时,为什么会大致形成正态分布?任何数字都是等可能的,对吧?分布不应该是平坦的而不是钟形曲线吗?
发生了什么呢?样本中的单个数字本身不会产生正态分布。分布将是平坦的,其中任何数字都是等可能的(称为均匀分布)。但当我们将它们作为样本分组并求平均时,它们形成一个正态分布。
这是因为中心极限定理,它表明当我们对一个人口取足够大的样本,计算每个样本的均值,并将它们作为一个分布绘制时,会发生有趣的事情:
-
样本均值的平均值等于总体均值。
-
如果总体是正态的,那么样本均值将是正态的。
-
如果总体不是正态分布,但样本量大于 30,样本均值仍然会大致形成正态分布。
-
样本均值的标准差等于总体标准差除以n的平方根:
为什么上述所有内容很重要?这些行为使我们能够根据样本推断出关于人口的有用信息,即使对于非正态分布的人口也是如此。如果你修改前面的代码并尝试较小的样本量,比如 1 或 2,你将看不到正态分布出现。但当你接近 31 或更多时,你会看到我们收敛到一个正态分布,如图 3-13 所示。
![emds 0313]()
图 3-13。较大的样本量接近正态分布
三十一是统计学中的教科书数字,因为这时我们的样本分布通常会收敛到总体分布,特别是当我们测量样本均值或其他参数时。当你的样本少于 31 个时,这时你必须依赖 T 分布而不是正态分布,随着样本量的减小,正态分布的尾部会变得越来越厚。我们稍后会简要讨论这一点,但首先让我们假设在谈论置信区间和检验时,我们的样本至少有 31 个。
置信区间
你可能听过“置信区间”这个术语,这经常让统计新手和学生感到困惑。置信区间是一个范围计算,显示我们有多大信心相信样本均值(或其他参数)落在总体均值的范围内。
基于 31 只金毛猎犬的样本,样本均值为 64.408,样本标准差为 2.05,我有 95%的信心总体均值在 63.686 和 65.1296 之间。 我怎么知道这个?让我告诉你,如果你感到困惑,回到这段话,记住我们试图实现什么。我有意突出了它!
我首先选择一个置信水平(LOC),这将包含所需概率的总体均值范围。我希望有 95%的信心,我的样本均值落在我将计算的总体均值范围内。这就是我的 LOC。我们可以利用中心极限定理并推断出总体均值的这个范围是什么。首先,我需要关键 z 值,这是在标准正态分布中给出 95%概率的对称范围,如图 3-14 中所示。
我们如何计算包含 0.95 面积的对称范围?这作为一个概念比作为一个计算更容易理解。你可能本能地想使用 CDF,但随后你可能会意识到这里有更多的变动部分。
首先,您需要利用逆 CDF。从逻辑上讲,为了获得中心对称区域的 95%面积,我们将切掉剩余 5%面积的尾部。将剩余的 5%面积分成两半,每个尾部将给我们 2.5%的面积。因此,我们想要查找 x 值的面积是 0.025 和 0.975,如图 3-15 所示。
![emds 0314]()
图 3-14。标准正态分布中心的 95%对称概率
![emds 0315]()
图 3-15。我们想要给出面积为 0.025 和 0.975 的 x 值
我们可以查找面积为 0.025 和面积为 0.975 的 x 值,这将给出包含 95%面积的中心范围。然后我们将返回包含这个区域的相应下限和上限 z 值。请记住,我们在这里使用标准正态分布,因此它们除了是正/负之外是相同的。让我们按照 Python 中示例 3-16 中所示的方式计算这个。
示例 3-16。检索关键的 z 值
from scipy.stats import norm
def critical_z_value(p):
norm_dist = norm(loc=0.0, scale=1.0)
left_tail_area = (1.0 - p) / 2.0
upper_area = 1.0 - ((1.0 - p) / 2.0)
return norm_dist.ppf(left_tail_area), norm_dist.ppf(upper_area)
print(critical_z_value(p=.95))
# (-1.959963984540054, 1.959963984540054)
好的,所以我们得到了±1.95996,这是我们在标准正态分布中心捕获 95%概率的关键 z 值。接下来,我将利用中心极限定理来产生误差边界(E),这是样本均值周围包含在该置信水平下的总体均值的范围。回想一下,我们的 31 只金毛猎犬样本的均值为 64.408,标准差为 2.05。计算这个误差边界的公式是:
如果我们将这个误差边界应用于样本均值,最终我们就得到了置信区间!
这里是如何在 Python 中从头到尾计算这个置信区间的,详见示例 3-17。
示例 3-17. 在 Python 中计算置信区间
from math import sqrt
from scipy.stats import norm
def critical_z_value(p):
norm_dist = norm(loc=0.0, scale=1.0)
left_tail_area = (1.0 - p) / 2.0
upper_area = 1.0 - ((1.0 - p) / 2.0)
return norm_dist.ppf(left_tail_area), norm_dist.ppf(upper_area)
def confidence_interval(p, sample_mean, sample_std, n):
# Sample size must be greater than 30
lower, upper = critical_z_value(p)
lower_ci = lower * (sample_std / sqrt(n))
upper_ci = upper * (sample_std / sqrt(n))
return sample_mean + lower_ci, sample_mean + upper_ci
print(confidence_interval(p=.95, sample_mean=64.408, sample_std=2.05, n=31))
# (63.68635915701992, 65.12964084298008)
因此,解释这个的方式是“基于我的 31 只金毛犬体重样本,样本均值为 64.408,样本标准差为 2.05,我有 95%的信心总体均值在 63.686 和 65.1296 之间。” 这就是我们描述置信区间的方式。
这里还有一件有趣的事情要注意,就是在我们的误差边界公式中,n 变大时,我们的置信区间变得更窄!这是有道理的,因为如果我们有一个更大的样本,我们对于总体均值落在更小范围内的信心就更大,这就是为什么它被称为置信区间。
这里需要注意的一个警告是,为了使这个工作起效,我们的样本大小必须至少为 31 个项目。这回到了中心极限定理。如果我们想将置信区间应用于更小的样本,我们需要使用方差更高的分布(更多不确定性的更胖尾巴)。这就是 T 分布的用途,我们将在本章末讨论它。
在第五章中,我们将继续使用线性回归的置信区间。
理解 P 值
当我们说某事是统计显著时,我们指的是什么?我们经常听到这个词被随意地使用,但在数学上它意味着什么?从技术上讲,这与一个叫做 p 值的东西有关,这对许多人来说是一个难以理解的概念。但我认为当你将 p 值的概念追溯到它的发明时,它会更有意义。虽然这只是一个不完美的例子,但它传达了一些重要的思想。
1925 年,数学家罗纳德·费舍尔在一个聚会上。他的同事穆里尔·布里斯托尔声称她能通过品尝来区分是先倒茶还是先倒牛奶。对这一说法感到好奇,罗纳德当场设置了一个实验。
他准备了八杯茶。四杯先倒牛奶,另外四杯先倒茶。然后他把它们呈现给他的鉴赏家同事,并要求她识别每一杯的倒入顺序。令人惊讶的是,她全部正确识别了,这种情况发生的概率是 70 分之 1,或者 0.01428571。
这个 1.4%的概率就是我们所谓的p 值,即某事发生的概率是偶然而非因为一个假设的解释。不深入组合数学的兔子洞,穆里尔完全猜对杯子的概率是 1.4%。这到底告诉了你什么?
当我们设计一个实验,无论是确定有机甜甜圈是否导致体重增加,还是住在电力线附近是否导致癌症,我们总是要考虑到随机运气起了作用的可能性。就像穆里尔仅仅通过猜测正确识别茶杯的概率为 1.4%一样,总是有一种可能性,即随机性给了我们一个好手牌,就像老丨虎丨机一样。这有助于我们构建我们的零假设(H[0]),即所研究的变量对实验没有影响,任何积极的结果只是随机运气。备择假设(H[1]) 提出所研究的变量(称为控制变量) 导致了积极的结果。
传统上,统计显著性的阈值是 5%或更低的 p 值,即 0.05。由于 0.014 小于 0.05,这意味着我们可以拒绝零假设,即穆里尔是随机猜测的。然后我们可以提出备择假设,即穆里尔有特殊能力判断是先倒茶还是牛奶。
现在,这个茶会的例子没有捕捉到的一点是,当我们计算 p 值时,我们捕捉到了所有该事件或更罕见事件的概率。当我们深入研究下一个使用正态分布的例子时,我们将解决这个问题。
假设检验
过去的研究表明,感冒的平均恢复时间为 18 天,标准偏差为 1.5 天,并且符合正态分布。
这意味着在 15 到 21 天之间恢复的概率约为 95%,如图 3-16 和示例 3-18 所示。
![emds 0316]()
图 3-16。在 15 到 21 天之间有 95%的恢复机会。
示例 3-18。计算在 15 到 21 天之间恢复的概率
from scipy.stats import norm
# Cold has 18 day mean recovery, 1.5 std dev
mean = 18
std_dev = 1.5
# 95% probability recovery time takes between 15 and 21 days.
x = norm.cdf(21, mean, std_dev) - norm.cdf(15, mean, std_dev)
print(x) # 0.9544997361036416
我们可以从剩下的 5%概率推断出,恢复需要超过 21 天的概率为 2.5%,而少于 15 天的概率也为 2.5%。记住这一点,因为这将在后面至关重要!这驱动了我们的 p 值。
现在假设一个实验性新药物被给予了一个由 40 人组成的测试组,他们从感冒中恢复需要平均 16 天,如图 3-17 所示。
![emds 0317]()
图 3-17。服用药物的一组人恢复需要 16 天。
这种药物有影响吗?如果你思考足够长时间,你可能会意识到我们所问的是:这种药物是否显示出统计上显著的结果?或者这种药物没有起作用,16 天的恢复只是与测试组的巧合?第一个问题构成了我们的备择假设,而第二个问题构成了我们的零假设。
我们可以通过两种方式来计算这个问题:单尾检验和双尾检验。我们将从单尾检验开始。
单尾检验
当我们进行单尾检验时,通常使用不等式来构建零假设和备择假设。我们假设围绕着总体均值,并说它要么大于/等于 18(零假设 ),要么小于 18(备择假设 ):
要拒绝我们的零假设,我们需要表明服用药物的患者的样本均值不太可能是巧合的。由于传统上认为小于或等于.05 的 p 值在统计上是显著的,我们将使用这个作为我们的阈值(图 3-17)。当我们在 Python 中使用逆 CDF 计算时,如图 3-18 和示例 3-19 所示,我们发现大约 15.53 是给我们左尾部分.05 面积的恢复天数。
![emds 0318]()
图 3-18. 获取其后面 5%面积的 x 值
示例 3-19. 获取其后面 5%面积的 x 值的 Python 代码
from scipy.stats import norm
# Cold has 18 day mean recovery, 1.5 std dev
mean = 18
std_dev = 1.5
# What x-value has 5% of area behind it?
x = norm.ppf(.05, mean, std_dev)
print(x) # 15.53271955957279
因此,如果我们在样本组中实现平均 15.53 天或更少的恢复时间,我们的药物被认为在统计上足够显著地显示了影响。然而,我们的恢复时间的样本均值实际上是 16 天,不在这个零假设拒绝区域内。因此,如图 3-19 所示,统计显著性测试失败了。
![emds 0319]()
图 3-19. 我们未能证明我们的药物测试结果在统计上是显著的
到达那个 16 天标记的面积是我们的 p 值,为.0912,并且我们在 Python 中计算如示例 3-20 所示。
示例 3-20. 计算单尾 p 值
from scipy.stats import norm
# Cold has 18 day mean recovery, 1.5 std dev
mean = 18
std_dev = 1.5
# Probability of 16 or less days
p_value = norm.cdf(16, mean, std_dev)
print(p_value) # 0.09121121972586788
由于.0912 的 p 值大于我们的统计显著性阈值.05,我们不认为药物试验是成功的,并且未能拒绝我们的零假设。
双尾检验
我们之前进行的测试称为单尾检验,因为它只在一个尾部寻找统计显著性。然而,通常更安全和更好的做法是使用双尾检验。我们将详细说明原因,但首先让我们计算它。
要进行双尾检验,我们将零假设和备择假设构建在“相等”和“不相等”的结构中。在我们的药物测试中,我们会说零假设的平均恢复时间为 18 天。但我们的备择假设是平均恢复时间不是 18 天,这要归功于新药物:
这有一个重要的含义。我们构建我们的备择假设不是测试药物是否改善感冒恢复时间,而是是否有任何影响。这包括测试它是否增加了感冒的持续时间。这有帮助吗?记住这个想法。
自然地,这意味着我们将 p 值的统计显著性阈值扩展到两个尾部,而不仅仅是一个。如果我们正在测试 5% 的统计显著性,那么我们将其分割,并将每个尾部分配 2.5%。如果我们的药物的平均恢复时间落在任一区域内,我们的测试将成功,并拒绝零假设(图 3-20)。
![emds 0320]()
图 3-20. 双尾检验
下尾和上尾的 x 值分别为 15.06 和 20.93,这意味着如果我们低于或高于这些值,我们将拒绝零假设。这两个值是使用图 3-21 和示例 3-21 中显示的逆 CDF 计算的。记住,为了得到上尾,我们取 .95 然后加上显著性阈值的 .025 部分,得到 .975。
![emds 0321]()
图 3-21. 计算正态分布中心 95% 的区域
示例 3-21. 计算 5% 统计显著性范围
from scipy.stats import norm
# Cold has 18 day mean recovery, 1.5 std dev
mean = 18
std_dev = 1.5
# What x-value has 2.5% of area behind it?
x1 = norm.ppf(.025, mean, std_dev)
# What x-value has 97.5% of area behind it
x2 = norm.ppf(.975, mean, std_dev)
print(x1) # 15.060054023189918
print(x2) # 20.93994597681008
药物测试组的样本均值为 16,16 既不小于 15.06 也不大于 20.9399。因此,就像单尾检验一样,我们仍然未能拒绝零假设。我们的药物仍然没有显示出任何统计显著性,也没有任何影响,如图 3-22 所示。
![emds 0322]()
图 3-22. 双尾检验未能证明统计显著性
但是 p 值是什么?这就是双尾检验变得有趣的地方。我们的 p 值将捕捉不仅仅是 16 左侧的区域,还将捕捉右尾的对称等效区域。由于 16 比均值低 4 天,我们还将捕捉高于 20 的区域,即比均值高 4 天(图 3-23)。这是捕捉事件或更罕见事件的概率,位于钟形曲线的两侧。
![emds 0323]()
图 3-23. p 值添加对称侧以获得统计显著性
当我们将这两个区域相加时,我们得到一个 p 值为 .1824。这远大于 .05,因此绝对不符合我们的 p 值阈值为 .05(示例 3-22)。
示例 3-22. 计算双尾 p 值
from scipy.stats import norm
# Cold has 18 day mean recovery, 1.5 std dev
mean = 18
std_dev = 1.5
# Probability of 16 or less days
p1 = norm.cdf(16, mean, std_dev)
# Probability of 20 or more days
p2 = 1.0 - norm.cdf(20, mean, std_dev)
# P-value of both tails
p_value = p1 + p2
print(p_value) # 0.18242243945173575
那么为什么在双尾检验中还要添加对称区域的相反侧?这可能不是最直观的概念,但首先记住我们如何构建我们的假设:
如果我们在“等于 18”与“不等于 18”的情况下进行测试,我们必须捕捉两侧任何等于或小于的概率。毕竟,我们试图证明显著性,这包括任何同样或更不可能发生的事情。我们在只使用“大于/小于”逻辑的单尾检验中没有这种特殊考虑。但当我们处理“等于/不等于”时,我们的兴趣范围向两个方向延伸。
那么双尾检验有什么实际影响?它如何影响我们是否拒绝零假设?问问自己:哪一个设定了更高的阈值?你会注意到,即使我们的目标是表明我们可能减少了某些东西(使用药物减少感冒恢复时间),重新构建我们的假设以显示任何影响(更大或更小)会创建一个更高的显著性阈值。如果我们的显著性阈值是小于或等于 0.05 的 p 值,我们的单尾检验在 p 值 0.0912 时更接近接受,而双尾检验则是大约两倍的 p 值 0.182。
这意味着双尾检验更难拒绝零假设,并需要更强的证据来通过测试。还要考虑这一点:如果我们的药物可能加重感冒并使其持续时间更长呢?捕捉这种可能性并考虑该方向的变化可能会有所帮助。这就是为什么在大多数情况下双尾检验更可取。它们往往更可靠,不会偏向于某一个方向的假设。
我们将在第五章和第六章再次使用假设检验和 p 值。
T 分布:处理小样本
让我们简要讨论如何处理 30 个或更少的较小样本;当我们在第五章进行线性回归时,我们将需要这个。无论是计算置信区间还是进行假设检验,如果样本中有 30 个或更少的项目,我们会选择使用 T 分布而不是正态分布。T 分布类似于正态分布,但尾部更胖,以反映更多的方差和不确定性。图 3-24 显示了一个正态分布(虚线)和一个自由度为 1 的 T 分布(实线)。
![emds 0324]()
图 3-24。T 分布与正态分布并列;请注意更胖的尾部。
样本大小越小,T 分布的尾部就越胖。但有趣的是,在接近 31 个项目后,T 分布几乎无法与正态分布区分开来,这很好地反映了中心极限定理的思想。
示例 3-23 展示了如何找到 95%置信度的临界 t 值。当样本量为 30 或更少时,您可以在置信区间和假设检验中使用这个值。概念上与关键 z 值相同,但我们使用 T 分布而不是正态分布来反映更大的不确定性。样本量越小,范围越大,反映了更大的不确定性。
示例 3-23. 用 T 分布获取临界值范围
from scipy.stats import t
# get critical value range for 95% confidence
# with a sample size of 25
n = 25
lower = t.ppf(.025, df=n-1)
upper = t.ppf(.975, df=n-1)
print(lower, upper)
-2.063898561628021 2.0638985616280205
请注意,df
是“自由度”参数,正如前面所述,它应该比样本量少一个。
大数据考虑和得克萨斯神枪手谬误
在我们结束本章之前,最后一个想法。正如我们所讨论的,随机性在验证我们的发现中起着如此重要的作用,我们总是要考虑到它的可能性。不幸的是,随着大数据、机器学习和其他数据挖掘工具的出现,科学方法突然变成了一种倒退的实践。这可能是危险的;请允许我演示一下,从加里·史密斯的书标准偏差(Overlook Press)中借鉴一个例子。
假设我从一副公平的牌中抽取四张牌。这里没有游戏或目标,只是抽取四张牌并观察它们。我得到两个 10,一个 3 和一个 2。我说:“这很有趣,我得到两个 10,一个 3 和一个 2。这有意义吗?接下来我抽的四张牌也会是两个连续的数字和一对吗?这里的潜在模型是什么?”
看到我做了什么吗?我拿了完全随机的东西,不仅寻找了模式,而且试图从中建立一个预测模型。这里微妙发生的是,我从未把得到这四张具有特定模式的牌作为我的目标。我是在它们发生之后观察到它们。
这正是数据挖掘每天都会遇到的问题:在随机事件中找到巧合的模式。随着大量数据和快速算法寻找模式,很容易找到看起来有意义但实际上只是随机巧合的东西。
这也类似于我朝墙壁开枪。然后我在弹孔周围画一个靶子,然后邀请朋友过来炫耀我的惊人射击技巧。荒谬,对吧?嗯,很多数据科学家每天都在比喻性地做这件事,这就是得克萨斯神枪手谬误。他们开始行动而没有目标,偶然发现了一些罕见的东西,然后指出他们发现的东西在某种程度上具有预测价值。
问题在于大数定律表明罕见事件很可能会发生;我们只是不知道是哪些。当我们遇到罕见事件时,我们会突出显示甚至推测可能导致它们的原因。问题在于:某个特定人赢得彩票的概率极低,但总会有人赢得彩票。当有人赢得彩票时,我们为什么会感到惊讶?除非有人预测了赢家,否则除了一个随机的人走运之外,没有发生任何有意义的事情。
这也适用于相关性,我们将在第五章中学习。在具有数千个变量的庞大数据集中,很容易找到具有 0.05 p 值的统计显著性发现吗?当然!我会找到成千上万个这样的发现。我甚至会展示证据表明尼古拉斯·凯奇电影的数量与一年中溺水的人数相关。
因此,为了防止得到德克萨斯神枪手谬误并成为大数据谬论的受害者,请尝试使用结构化的假设检验并为此目标收集数据。如果您使用数据挖掘,请尝试获取新鲜数据以查看您的发现是否仍然成立。最后,始终考虑事情可能是巧合的可能性;如果没有常识性的解释,那么很可能是巧合。
我们学会了在收集数据之前进行假设,但数据挖掘是先收集数据,然后进行假设。具有讽刺意味的是,我们在开始时更客观,因为我们会寻找数据来有意识地证明和反驳我们的假设。
结论
我们在本章学到了很多,你应该为走到这一步感到自豪。这可能是本书中较难的主题之一!我们不仅学习了从平均值到正态分布的描述性统计,还解决了置信区间和假设检验的问题。
希望你能稍微不同地看待数据。数据是某种东西的快照,而不是对现实的完全捕捉。单独的数据并不是很有用,我们需要背景、好奇心和分析数据来源的地方,然后我们才能对其进行有意义的洞察。我们讨论了如何描述数据以及根据样本推断更大人口的属性。最后,我们解决了一些数据挖掘中可能遇到的谬论,以及如何通过新鲜数据和常识来纠正这些谬论。
如果您需要回顾本章内容,不要感到难过,因为需要消化的内容很多。如果您想在数据科学和机器学习职业中取得成功,进入假设检验的思维模式也很重要。很少有从业者将统计和假设检验概念与机器学习联系起来,这是不幸的。
理解性和可解释性是机器学习的下一个前沿,因此在阅读本书的其余部分和职业生涯中继续学习和整合这些想法。
练习
-
你购买了一卷 1.75 毫米的线材用于你的 3D 打印机。你想测量线材直径与 1.75 毫米的实际接近程度。你使用卡钳工具在卷轴上抽取了五次直径值:
1.78, 1.75, 1.72, 1.74, 1.77
计算这组数值的均值和标准差。
-
一家制造商表示 Z-Phone 智能手机的平均使用寿命为 42 个月,标准差为 8 个月。假设服从正态分布,一个随机的 Z-Phone 使用寿命在 20 到 30 个月之间的概率是多少?
-
我怀疑我的 3D 打印机线材的平均直径不是广告中宣传的 1.75 毫米。我用我的工具抽取了 34 个测量值。样本均值为 1.715588,样本标准差为 0.029252。
我整个线轴的均值的 99% 置信区间是多少?
-
你的营销部门开始了一项新的广告活动,并想知道它是否影响了销售。过去销售额平均为每天 $10,345,标准差为 $552。新的广告活动持续了 45 天,平均销售额为 $11,641。
这次广告活动是否影响了销售?为什么?(使用双尾检验以获得更可靠的显著性。)
答案在附录 B 中。
第四章:线性代数
稍微转换一下思路,让我们远离概率和统计,进入线性代数领域。有时人们会混淆线性代数和基本代数,可能认为它与使用代数函数y = mx + b绘制线条有关。这就是为什么线性代数可能应该被称为“向量代数”或“矩阵代数”,因为它更加抽象。线性系统发挥作用,但以一种更加形而上的方式。
那么,线性代数到底是什么?嗯,线性代数关注线性系统,但通过向量空间和矩阵来表示它们。如果你不知道什么是向量或矩阵,不用担心!我们将深入定义和探索它们。线性代数对于许多应用领域的数学、统计学、运筹学、数据科学和机器学习都至关重要。当你在这些领域中处理数据时,你正在使用线性代数,也许你甚至不知道。
你可以暂时不学习线性代数,使用机器学习和统计库为你完成所有工作。但是,如果你想理解这些黑匣子背后的直觉,并更有效地处理数据,理解线性代数的基础是不可避免的。线性代数是一个庞大的主题,可以填满厚厚的教科书,所以当然我们不能在这本书的一章中完全掌握它。然而,我们可以学到足够多,以便更加熟练地应用它并有效地在数据科学领域中导航。在本书的剩余章节中,包括第五章和第七章,也将有机会应用它。
什么是向量?
简单来说,向量是空间中具有特定方向和长度的箭头,通常代表一段数据。它是线性代数的中心构建模块,包括矩阵和线性变换。在其基本形式中,它没有位置的概念,所以始终想象它的尾部从笛卡尔平面的原点(0,0)开始。
图 4-1 展示了一个向量,它在水平方向移动三步,在垂直方向移动两步。
![emds 0401]()
图 4-1。一个简单的向量
再次强调,向量的目的是直观地表示一段数据。如果你有一条房屋面积为 18,000 平方英尺,估值为 26 万美元的数据记录,我们可以将其表示为一个向量[18000, 2600000],在水平方向移动 18000 步,在垂直方向移动 260000 步。
我们数学上声明一个向量如下:
我们可以使用简单的 Python 集合,如 Python 列表,在示例 4-1 中所示声明一个向量。
示例 4-1. 使用列表在 Python 中声明向量
v = [3, 2]
print(v)
然而,当我们开始对向量进行数学计算,特别是在执行诸如机器学习之类的任务时,我们应该使用 NumPy 库,因为它比纯 Python 更高效。您还可以使用 SymPy 执行线性代数运算,在本章中,当小数变得不方便时,我们偶尔会使用它。然而,在实践中,您可能主要使用 NumPy,因此我们主要会坚持使用它。
要声明一个向量,您可以使用 NumPy 的array()
函数,然后可以像示例 4-2 中所示传递一组数字给它。
示例 4-2. 使用 NumPy 在 Python 中声明向量
import numpy as np
v = np.array([3, 2])
print(v)
Python 速度慢,但其数值库不慢
Python 是一个计算速度较慢的语言平台,因为它不像 Java、C#、C 等编译为较低级别的机器代码和字节码。它在运行时动态解释。然而,Python 的数值和科学库并不慢。像 NumPy 这样的库通常是用低级语言如 C 和 C++编写的,因此它们在计算上是高效的。Python 实际上充当“胶水代码”,为您的任务集成这些库。
向量有无数的实际应用。在物理学中,向量通常被认为是一个方向和大小。在数学中,它是 XY 平面上的一个方向和比例,有点像运动。在计算机科学中,它是存储数据的一组数字。作为数据科学专业人员,我们将最熟悉计算机科学的上下文。然而,重要的是我们永远不要忘记视觉方面,这样我们就不会把向量看作是神秘的数字网格。没有视觉理解,几乎不可能掌握许多基本的线性代数概念,如线性相关性和行列式。
这里有一些更多的向量示例。在图 4-2 中,请注意一些向量在 X 和 Y 轴上具有负方向。具有负方向的向量在我们后面合并时会产生影响,基本上是相减而不是相加。
![emds 0402]()
图 4-2. 不同向量的抽样
还要注意,向量可以存在于超过两个维度。接下来我们声明一个沿着 x、y 和 z 轴的三维向量:
要创建这个向量,我们在 x 方向走了四步,在 y 方向走了一步,在 z 方向走了两步。这在图 4-3 中有可视化展示。请注意,我们不再在二维网格上显示向量,而是在一个三维空间中,有三个轴:x、y 和 z。
![emds 0403]()
图 4-3. 一个三维向量
当然,我们可以使用三个数值在 Python 中表示这个三维向量,就像在示例 4-3 中声明的那样。
示例 4-3. 在 Python 中使用 NumPy 声明一个三维向量
import numpy as np
v = np.array([4, 1, 2])
print(v)
像许多数学模型一样,可视化超过三维是具有挑战性的,这是我们在本书中不会花费精力去做的事情。但从数字上来看,仍然很简单。示例 4-4 展示了我们如何在 Python 中数学上声明一个五维向量。
示例 4-4. 在 Python 中使用 NumPy 声明一个五维向量
import numpy as np
v = np.array([6, 1, 5, 8, 3])
print(v)
添加和组合向量
单独看,向量并不是非常有趣。它表示一个方向和大小,有点像在空间中的移动。但当你开始组合向量,也就是向量相加时,事情就变得有趣起来。我们实际上将两个向量的运动合并成一个单一的向量。
假设我们有两个向量 和 如图 4-4 所示。我们如何将这两个向量相加呢?
![emds 0404]()
图 4-4. 将两个向量相加
我们将在稍后讨论为什么添加向量是有用的。但如果我们想要结合这两个向量,包括它们的方向和大小,那会是什么样子呢?从数字上来看,这很简单。你只需将各自的 x 值相加,然后将 y 值相加,得到一个新的向量,如示例 4-5 所示。
示例 4-5. 使用 NumPy 在 Python 中将两个向量相加
from numpy import array
v = array([3,2])
w = array([2,-1])
# sum the vectors
v_plus_w = v + w
# display summed vector
print(v_plus_w) # [5, 1]
但这在视觉上意味着什么呢?要将这两个向量视觉上相加,将一个向量连接到另一个向量的末端,然后走到最后一个向量的顶端(图 4-5)。你最终停留的位置就是一个新向量,是这两个向量相加的结果。
![emds 0405]()
图 4-5. 将两个向量相加得到一个新向量
如图 4-5 所示,当我们走到最后一个向量 的末端时,我们得到一个新向量 [5, 1]。这个新向量是将 和 相加的结果。在实践中,这可以简单地将数据相加。如果我们在一个区域中总结房屋价值和其平方英尺,我们将以这种方式将多个向量相加成一个单一向量。
请注意,无论我们是先将 加上 还是反之,都不影响结果,这意味着它是可交换的,操作顺序不重要。如果我们先走 再走 ,我们最终得到的向量 [5, 1] 与 图 4-6 中可视化的相同。
![emds 0406]()
图 4-6. 向量的加法是可交换的
缩放向量
缩放是增加或减小向量的长度。你可以通过乘以一个称为标量的单个值来增加/减小向量。图 4-7 是向量被放大 2 倍的示例。
![emds 0407]()
图 4-7. 向量的缩放
从数学上讲,你将向量的每个元素乘以标量值:
在 Python 中执行这个缩放操作就像将向量乘以标量一样简单,就像在示例 4-6 中编写的那样。
示例 4-6. 在 Python 中使用 NumPy 缩放数字
from numpy import array
v = array([3,1])
# scale the vector
scaled_v = 2.0 * v
# display scaled vector
print(scaled_v) # [6 2]
在这里,图 4-8 中的被缩小了一半。
![emds 0408]()
图 4-8. 将向量缩小一半
这里需要注意的一个重要细节是,缩放向量不会改变其方向,只会改变其大小。但是有一个轻微的例外,如图 4-9 中所示。当你用一个负数乘以一个向量时,它会翻转向量的方向。
![emds 0409]()
图 4-9. 负标量会翻转向量方向
仔细思考一下,通过负数进行缩放实际上并没有改变方向,因为它仍然存在于同一条线上。这引出了一个称为线性相关的关键概念。
跨度和线性相关性
这两个操作,相加两个向量并对它们进行缩放,带来了一个简单但强大的想法。通过这两个操作,我们可以组合两个向量并对它们进行缩放,以创建我们想要的任何结果向量。图 4-10 展示了取两个向量和,进行缩放和组合的六个示例。这些向量和,固定在两个不同方向上,可以被缩放和相加以创建任何新向量。
![emds 0410]()
图 4-10。缩放两个相加的向量允许我们创建任何新向量
再次,和在方向上固定,除了使用负标量进行翻转,但我们可以使用缩放自由地创建由组成的任何向量。
整个可能向量空间称为span,在大多数情况下,我们的 span 可以通过缩放和求和这两个向量来创建无限数量的向量。当我们有两个不同方向的向量时,它们是线性独立的,并且具有这种无限的 span。
但在哪种情况下我们受限于可以创建的向量?思考一下并继续阅读。
当两个向量存在于相同方向,或存在于同一条线上时会发生什么?这些向量的组合也被限制在同一条线上,将我们的 span 限制在那条线上。无论如何缩放,结果的和向量也被困在同一条线上。这使它们成为线性相关,如图 4-11 所示。
![emds 0411]()
图 4-11。线性相关向量
这里的 span 被困在与由两个向量组成的同一条线上。因为这两个向量存在于相同的基础线上,我们无法通过缩放灵活地创建任何新向量。
在三维或更高维空间中,当我们有一组线性相关的向量时,我们经常会被困在较低维度的平面上。这里有一个例子,即使我们有三维向量如图 4-12 所述,我们也被困在二维平面上。
![emds 0412]()
图 4-12。三维空间中的线性相关性;请注意我们的张成空间被限制在一个平面上
后面我们将学习一个简单的工具叫做行列式来检查线性相关性,但我们为什么关心两个向量是线性相关还是线性无关呢?当它们线性相关时,很多问题会变得困难或无法解决。例如,当我们在本章后面学习方程组时,一个线性相关的方程组会导致变量消失,使问题无法解决。但如果你有线性无关性,那么从两个或更多向量中创建任何你需要的向量的灵活性将变得无价,以便解决问题!
线性变换
将具有固定方向的两个向量相加,但按比例缩放它们以获得不同的组合向量的概念非常重要。这个组合向量,除了线性相关的情况外,可以指向任何方向并具有我们选择的任何长度。这为线性变换建立了一种直观,我们可以使用一个向量以类似函数的方式来转换另一个向量。
基向量
想象我们有两个简单的向量和(“i-hat”和“j-hat”)。这些被称为基向量,用于描述对其他向量的变换。它们通常长度为 1,并且指向垂直正方向,如图 4-13 中可视化的。
![emds 0413]()
图 4-13。基向量和
将基向量视为构建或转换任何向量的构建块。我们的基向量用一个 2×2 矩阵表示,其中第一列是,第二列是:
矩阵是一组向量(如,),可以有多行多列,是打包数据的便捷方式。我们可以使用和通过缩放和相加来创建任何我们想要的向量。让我们从长度为 1 开始,并在图 4-14 中展示结果向量。
![emds 0414]()
图 4-14. 从基向量创建向量
我希望向量落在[3, 2]处。如果我们将拉伸 3 倍,拉伸 2 倍,那么会发生什么?首先,我们分别按照下面的方式对它们进行缩放:
如果我们在这两个方向上拉伸空间,那么这对 会有什么影响呢?嗯,它将会随着 和 一起拉伸。这被称为线性变换,我们通过跟踪基向量的移动来进行向量的拉伸、挤压、剪切或旋转。在这种情况下(图 4-15),缩放 和 已经沿着我们的向量 拉伸了空间。
![emds 0415]()
图 4-15. 一个线性变换
但是 会落在哪里呢?很容易看出它会落在这里,即 [3, 2]。记住向量 是由 和 相加而成的。因此,我们只需将拉伸的 和 相加,就能看出向量 落在哪里:
一般来说,通过线性变换,你可以实现四种运动,如图 4-16 所示。
![emds 0416]()
图 4-16。线性变换可以实现四种运动
这四种线性变换是线性代数的核心部分。缩放向量将拉伸或挤压它。旋转将转动向量空间,而反转将翻转向量空间,使和交换位置。
值得注意的是,你不能有非线性的变换,导致曲线或波浪状的变换不再遵守直线。这就是为什么我们称其为线性代数,而不是非线性代数!
矩阵向量乘法
这将引出线性代数中的下一个重要概念。在变换后跟踪和落在哪里的概念很重要,因为这不仅允许我们创建向量,还可以变换现有向量。如果你想要真正的线性代数启示,想想为什么创建向量和变换向量实际上是相同的事情。这完全取决于相对性,考虑到你的基向量在变换前后都是起点。
给定作为矩阵打包的基向量和,变换向量的公式是:
是第一列[a, c], 是第二列[b, d]。我们将这两个基向量打包成一个矩阵,再次表示为一个在两个或更多维度中以数字网格形式表达的向量集合。通过应用基向量对向量进行转换被称为矩阵向量乘法。这一开始可能看起来有些牵强,但这个公式是对缩放和添加和的快捷方式,就像我们之前对两个向量进行加法,并将该转换应用于任何向量。
因此,实际上,矩阵确实是表示为基向量的转换。
要在 Python 中使用 NumPy 执行这种转换,我们需要将基向量声明为矩阵,然后使用dot()
运算符将其应用于向量(参见示例 4-7)。dot()
运算符将执行我们刚刚描述的矩阵和向量之间的缩放和加法。这被称为点积,我们将在本章中探讨它。
示例 4-7. NumPy 中的矩阵向量乘法
from numpy import array
# compose basis matrix with i-hat and j-hat
basis = array(
[[3, 0],
[0, 2]]
)
# declare vector v
v = array([1,1])
# create new vector
# by transforming v with dot product
new_v = basis.dot(v)
print(new_v) # [3, 2]
当考虑基向量时,我更喜欢将基向量分解然后将它们组合成一个矩阵。只需注意你需要转置,或者交换列和行。这是因为 NumPy 的array()
函数会执行我们不希望的相反方向,将每个向量填充为一行而不是一列。在 NumPy 中的转置示例可在示例 4-8 中看到。
示例 4-8. 分离基向量并将它们应用为一个转换
from numpy import array
# Declare i-hat and j-hat
i_hat = array([2, 0])
j_hat = array([0, 3])
# compose basis matrix using i-hat and j-hat
# also need to transpose rows into columns
basis = array([i_hat, j_hat]).transpose()
# declare vector v
v = array([1,1])
# create new vector
# by transforming v with dot product
new_v = basis.dot(v)
print(new_v) # [2, 3]
这里有另一个例子。让我们从向量为[2, 1]开始,和分别从[1, 0]和[0, 1]开始。然后我们将和转换为[2, 0]和[0, 3]。向量会发生什么?通过手工使用我们的公式进行数学计算,我们得到如下结果:
示例 4-9 展示了这个解决方案在 Python 中的应用。
示例 4-9。使用 NumPy 转换向量
from numpy import array
# Declare i-hat and j-hat
i_hat = array([2, 0])
j_hat = array([0, 3])
# compose basis matrix using i-hat and j-hat
# also need to transpose rows into columns
basis = array([i_hat, j_hat]).transpose()
# declare vector v 0
v = array([2,1])
# create new vector
# by transforming v with dot product
new_v = basis.dot(v)
print(new_v) # [4, 3]
向量现在落在[4, 3]处。图 4-17 展示了这个变换的样子。
![emds 0417]()
图 4-17。一个拉伸的线性变换
这是一个将事情提升到一个新水平的例子。让我们取值为[2, 1]的向量。和起始于[1, 0]和[0, 1],但随后被变换并落在[2, 3]和[2, -1]。会发生什么?让我们在图 4-18 和示例 4-10 中看看。
![emds 0418]()
图 4-18。一个进行旋转、剪切和空间翻转的线性变换
示例 4-10。一个更复杂的变换
from numpy import array
# Declare i-hat and j-hat
i_hat = array([2, 3])
j_hat = array([2, -1])
# compose basis matrix using i-hat and j-hat
# also need to transpose rows into columns
basis = array([i_hat, j_hat]).transpose()
# declare vector v 0
v = array([2,1])
# create new vector
# by transforming v with dot product
new_v = basis.dot(v)
print(new_v) # [6, 5]
这里发生了很多事情。我们不仅缩放了和,还拉长了向量。我们实际上还剪切、旋转和翻转了空间。当和在顺时针方向上交换位置时,你会知道空间被翻转了,我们将在本章后面学习如何通过行列式检测这一点。
矩阵乘法
我们学会了如何将一个向量和一个矩阵相乘,但是将两个矩阵相乘到底实现了什么?将矩阵乘法看作是对向量空间应用多个变换。每个变换就像一个函数,我们首先应用最内层的变换,然后依次向外应用每个后续变换。
这里是我们如何对任意向量(值为[x, y])应用旋转和剪切:
我们实际上可以通过使用这个公式将这两个变换合并,将一个变换应用到最后一个上。你需要将第一个矩阵的每一行与第二个矩阵的每一列相乘并相加,按照“上下!上下!”的模式进行:
因此,我们实际上可以通过使用这个公式将这两个单独的变换(旋转和剪切)合并为一个单一的变换:
要在 Python 中使用 NumPy 执行这个操作,你可以简单地使用 matmul()
或 @
运算符来组合这两个矩阵(示例 4-11)。然后我们将转身并将这个合并的变换应用于一个向量 [1, 2]。
示例 4-11. 组合两个变换
from numpy import array
# Transformation 1
i_hat1 = array([0, 1])
j_hat1 = array([-1, 0])
transform1 = array([i_hat1, j_hat1]).transpose()
# Transformation 2
i_hat2 = array([1, 0])
j_hat2 = array([1, 1])
transform2 = array([i_hat2, j_hat2]).transpose()
# Combine Transformations
combined = transform2 @ transform1
# Test
print("COMBINED MATRIX:\n {}".format(combined))
v = array([1, 2])
print(combined.dot(v)) # [-1, 1]
使用 dot()
与 matmul()
以及 @
一般来说,你应该更倾向于使用 matmul()
和它的简写 @
来组合矩阵,而不是 NumPy 中的 dot()
运算符。前者通常对于高维矩阵和元素广播有更好的策略。
如果你喜欢深入研究这些实现细节,这个 StackOverflow 问题是一个很好的起点。
请注意,我们也可以将每个变换分别应用于向量 ,并且仍然会得到相同的结果。如果你用这三行替换最后一行,分别应用每个变换,你仍然会得到那个新向量 [-1, 1]:
rotated = transform1.dot(v)
sheered = transform2.dot(rotated)
print(sheered) # [-1, 1]
请注意,应用每个变换的顺序很重要!如果我们在 变换 2
上应用 变换 1
,我们会得到一个不同的结果 [-2, 3],如 示例 4-12 中计算的那样。因此矩阵点积不是可交换的,这意味着你不能改变顺序并期望得到相同的结果!
示例 4-12. 反向应用变换
from numpy import array
# Transformation 1
i_hat1 = array([0, 1])
j_hat1 = array([-1, 0])
transform1 = array([i_hat1, j_hat1]).transpose()
# Transformation 2
i_hat2 = array([1, 0])
j_hat2 = array([1, 1])
transform2 = array([i_hat2, j_hat2]).transpose()
# Combine Transformations, apply sheer first and then rotation
combined = transform1 @ transform2
# Test
print("COMBINED MATRIX:\n {}".format(combined))
v = array([1, 2])
print(combined.dot(v)) # [-2, 3]
将每个变换看作一个函数,并且我们从最内层到最外层应用它们,就像嵌套函数调用一样。
行列式
当我们进行线性变换时,有时会“扩展”或“挤压”空间,这种情况发生的程度可能会有所帮助。从向量空间中的采样区域中取出一个样本区域,看看在我们对 和 进行缩放后会发生什么?
![emds 0420]()
图 4-20. 一个行列式测量线性变换如何缩放一个区域
请注意,面积增加了 6.0 倍,这个因子被称为行列式。行列式描述了在向量空间中采样区域随线性变换的尺度变化,这可以提供有关变换的有用信息。
示例 4-13 展示了如何在 Python 中计算这个行列式。
示例 4-13. 计算行列式
from numpy.linalg import det
from numpy import array
i_hat = array([3, 0])
j_hat = array([0, 2])
basis = array([i_hat, j_hat]).transpose()
determinant = det(basis)
print(determinant) # prints 6.0
简单的剪切和旋转不会影响行列式,因为面积不会改变。图 4-21 和 示例 4-14 展示了一个简单的剪切,行列式仍然保持为 1.0,表明它没有改变。
![emds 0421]()
图 4-21. 简单的剪切不会改变行列式
示例 4-14. 一个剪切的行列式
from numpy.linalg import det
from numpy import array
i_hat = array([1, 0])
j_hat = array([1, 1])
basis = array([i_hat, j_hat]).transpose()
determinant = det(basis)
print(determinant) # prints 1.0
但是缩放会增加或减少行列式,因为这会增加/减少采样区域。当方向翻转时( , 顺时针交换位置),那么行列式将为负。图 4-22 和 示例 4-15 说明了一个行列式展示了一个不仅缩放而且翻转了向量空间方向的变换。
![emds 0422]()
图 4-22. 在翻转空间上的行列式是负的
示例 4-15. 一个负的行列式
from numpy.linalg import det
from numpy import array
i_hat = array([-2, 1])
j_hat = array([1, 2])
basis = array([i_hat, j_hat]).transpose()
determinant = det(basis)
print(determinant) # prints -5.0
因为这个行列式是负的,我们很快就看到方向已经翻转。但行列式告诉你的最关键的信息是变换是否线性相关。如果行列式为 0,那意味着所有的空间都被压缩到一个较低的维度。
在图 4-23 中,我们看到两个线性相关的变换,其中一个二维空间被压缩成一维,一个三维空间被压缩成两维。在这两种情况下,面积和体积分别为 0!
![emds 0423]()
图 4-23。二维和三维中的线性依赖
示例 4-16 展示了前述 2D 示例的代码,将整个二维空间压缩成一个一维数轴。
示例 4-16。行列式为零
from numpy.linalg import det
from numpy import array
i_hat = array([-2, 1])
j_hat = array([3, -1.5])
basis = array([i_hat, j_hat]).transpose()
determinant = det(basis)
print(determinant) # prints 0.0
因此,测试行列式为 0 对于确定一个变换是否具有线性依赖性非常有帮助。当你遇到这种情况时,你可能会发现手头上有一个困难或无法解决的问题。
特殊类型的矩阵
我们应该涵盖一些值得注意的矩阵情况。
方阵
方阵是行数和列数相等的矩阵:
它们主要用于表示线性变换,并且是许多操作的要求,如特征分解。
身份矩阵
身份矩阵是一个对角线为 1 的方阵,而其他值为 0:
身份矩阵有什么了不起的地方?当你有一个身份矩阵时,实际上是撤销了一个变换并找到了起始基向量。这在下一节解方程组中将起到重要作用。
逆矩阵
逆矩阵是一个可以撤销另一个矩阵变换的矩阵。假设我有矩阵A:
矩阵A的逆被称为。我们将在下一节学习如何使用 Sympy 或 NumPy 计算逆,但这就是矩阵A的逆:
当我们在和A之间进行矩阵乘法时,我们得到一个身份矩阵。我们将在下一节关于方程组的 NumPy 和 Sympy 中看到这一点。
对角矩阵
与身份矩阵类似的是对角矩阵,它在对角线上有非零值,而其余值为 0。对角矩阵在某些计算中是可取的,因为它们代表应用于向量空间的简单标量。它出现在一些线性代数操作中。
三角矩阵
与对角矩阵类似的是三角矩阵,它在一个三角形值的对角线前有非零值,而其余值为 0。
三角矩阵在许多数值分析任务中是理想的,因为它们通常更容易在方程组中解决。它们也出现在某些分解任务中,比如LU 分解。
稀疏矩阵
有时,你会遇到大部分元素为零且只有很少非零元素的矩阵。这些被称为稀疏矩阵。从纯数学的角度来看,它们并不是特别有趣。但从计算的角度来看,它们提供了创造效率的机会。如果一个矩阵大部分为 0,稀疏矩阵的实现将不会浪费空间存储大量的 0,而是只跟踪非零的单元格。
当你有大型稀疏矩阵时,你可能会明确使用稀疏函数来创建你的矩阵。
方程组和逆矩阵
线性代数的基本用例之一是解方程组。学习逆矩阵也是一个很好的应用。假设你被提供以下方程,需要解出x、y和z:
你可以尝试手动尝试不同的代数运算来分离三个变量,但如果你想让计算机解决它,你需要将这个问题用矩阵表示,如下所示。将系数提取到矩阵A中,方程右侧的值提取到矩阵B中,未知变量提取到矩阵X中:
线性方程组的函数是AX = B。我们需要用另一个矩阵X转换矩阵A,结果得到矩阵B:
我们需要“撤销”A,这样我们就可以隔离X并获得x、y和z的值。撤销A的方法是取A的逆,用表示,并通过矩阵乘法应用于A。我们可以用代数方式表达这一点:
要计算矩阵A的逆,我们可能会使用计算机,而不是手动使用高斯消元法寻找解决方案,这在本书中我们不会涉及。这是矩阵A的逆:
当我们将与A相乘时,将创建一个单位矩阵,一个除了对角线上为 1 之外全为零的矩阵。单位矩阵是线性代数中相当于乘以 1 的概念,意味着它基本上没有影响,将有效地隔离x、y和z的值:
要在 Python 中看到这个单位矩阵的运作,您将需要使用 SymPy 而不是 NumPy。 NumPy 中的浮点小数不会使单位矩阵显得那么明显,但在示例 4-17 中以符号方式进行,我们将看到一个清晰的符号输出。请注意,在 SymPy 中进行矩阵乘法时,我们使用星号 *** 而不是 @。
示例 4-17. 使用 SymPy 研究逆矩阵和单位矩阵
from sympy import *
# 4x + 2y + 4z = 44
# 5x + 3y + 7z = 56
# 9x + 3y + 6z = 72
A = Matrix([
[4, 2, 4],
[5, 3, 7],
[9, 3, 6]
])
# dot product between A and its inverse
# will produce identity function
inverse = A.inv()
identity = inverse * A
# prints Matrix([[-1/2, 0, 1/3], [11/2, -2, -4/3], [-2, 1, 1/3]])
print("INVERSE: {}".format(inverse))
# prints Matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
print("IDENTITY: {}".format(identity))
在实践中,浮点精度的缺失不会对我们的答案产生太大影响,因此使用 NumPy 来解决 x 应该是可以的。示例 4-18 展示了使用 NumPy 的解决方案。
示例 4-18. 使用 NumPy 解决一组方程
from numpy import array
from numpy.linalg import inv
# 4x + 2y + 4z = 44
# 5x + 3y + 7z = 56
# 9x + 3y + 6z = 72
A = array([
[4, 2, 4],
[5, 3, 7],
[9, 3, 6]
])
B = array([
44,
56,
72
])
X = inv(A).dot(B)
print(X) # [ 2\. 34\. -8.]
因此 x = 2,y = 34,z = –8. 示例 4-19 展示了 SymPy 中的完整解决方案,作为 NumPy 的替代方案。
示例 4-19. 使用 SymPy 解决一组方程
from sympy import *
# 4x + 2y + 4z = 44
# 5x + 3y + 7z = 56
# 9x + 3y + 6z = 72
A = Matrix([
[4, 2, 4],
[5, 3, 7],
[9, 3, 6]
])
B = Matrix([
44,
56,
72
])
X = A.inv() * B
print(X) # Matrix([[2], [34], [-8]])
这里是数学符号表示的解决方案:
希望这给了你逆矩阵的直觉以及它们如何用于解方程组的用途。
线性规划中的方程组
这种解方程组的方法也用于线性规划,其中不等式定义约束条件,目标是最小化/最大化。
PatrickJMT 在线性规划方面有很多优质视频。我们也在附录 A 中简要介绍了它。
在实际操作中,你很少需要手动计算逆矩阵,可以让计算机为你完成。但如果你有需求或好奇,你会想了解高斯消元法。PatrickJMT 在 YouTube 上的视频展示了许多关于高斯消元法的视频。
特征向量和特征值
矩阵分解是将矩阵分解为其基本组件,类似于分解数字(例如,10 可以分解为 2 × 5)。
矩阵分解对于诸如找到逆矩阵、计算行列式以及线性回归等任务非常有帮助。根据你的任务,有许多分解矩阵的方法。在第五章中,我们将使用一种矩阵分解技术,QR 分解,来执行线性回归。
但在本章中,让我们专注于一种常见方法,称为特征分解,这经常用于机器学习和主成分分析。在这个层面上,我们没有精力深入研究每个应用。现在,只需知道特征分解有助于将矩阵分解为在不同机器学习任务中更易处理的组件。还要注意它只适用于方阵。
在特征分解中,有两个组成部分:用 lambda 表示的特征值 和用 v 表示的特征向量,如 图 4-24 所示。
![emds 0424]()
图 4-24. 特征向量和特征值
如果我们有一个方阵 A,它有以下特征值方程:
如果 A 是原始矩阵,它由特征向量 和特征值 组成。对于父矩阵的每个维度,都有一个特征向量和特征值,并非所有矩阵都可以分解为特征向量和特征值。有时甚至会出现复数(虚数)。
示例 4-20 是我们如何在 NumPy 中为给定矩阵 计算特征向量和特征值的方法。
示例 4-20. 在 NumPy 中执行特征分解
from numpy import array, diag
from numpy.linalg import eig, inv
A = array([
[1, 2],
[4, 5]
])
eigenvals, eigenvecs = eig(A)
print("EIGENVALUES")
print(eigenvals)
print("\nEIGENVECTORS")
print(eigenvecs)
"""
EIGENVALUES
[-0.46410162 6.46410162]
EIGENVECTORS
[[-0.80689822 -0.34372377]
[ 0.59069049 -0.9390708 ]]
"""
那么我们如何从特征向量和特征值重建矩阵 A 呢?回想一下这个公式:
我们需要对重构 A 的公式进行一些调整:
在这个新公式中, 是特征向量, 是对角形式的特征值, 是 的逆矩阵。对角形式意味着向量被填充成一个零矩阵,并且占据对角线,类似于单位矩阵的模式。
示例 4-21 在 Python 中将示例完整地展示,从分解矩阵开始,然后重新组合它。
示例 4-21. 在 NumPy 中分解和重组矩阵
from numpy import array, diag
from numpy.linalg import eig, inv
A = array([
[1, 2],
[4, 5]
])
eigenvals, eigenvecs = eig(A)
print("EIGENVALUES")
print(eigenvals)
print("\nEIGENVECTORS")
print(eigenvecs)
print("\nREBUILD MATRIX")
Q = eigenvecs
R = inv(Q)
L = diag(eigenvals)
B = Q @ L @ R
print(B)
"""
EIGENVALUES
[-0.46410162 6.46410162]
EIGENVECTORS
[[-0.80689822 -0.34372377]
[ 0.59069049 -0.9390708 ]]
REBUILD MATRIX
[[1\. 2.]
[4\. 5.]]
"""
正如你所看到的,我们重建的矩阵就是我们开始的那个矩阵。
结论
线性代数可能令人困惑,充满了神秘和值得思考的想法。你可能会发现整个主题就像一个大的兔子洞,你是对的!然而,如果你想要有一个长期成功的数据科学职业,继续对它感到好奇是个好主意。它是统计计算、机器学习和其他应用数据科学领域的基础。最终,它是计算机科学的基础。你当然可以一段时间不了解它,但在某个时候你会遇到理解上的限制。
你可能会想知道这些理念如何实际应用,因为它们可能感觉很理论。不用担心;我们将在本书中看到一些实际应用。但是理论和几何解释对于在处理数据时具有直觉是很重要的,通过直观地理解线性变换,你就能为以后可能遇到的更高级概念做好准备。
如果你想更多地了解线性规划,没有比3Blue1Brown 的 YouTube 播放列表“线性代数的本质”更好的地方了。PatrickJMT 的线性代数视频也很有帮助。
如果你想更加熟悉 NumPy,推荐阅读 Wes McKinney 的 O'Reilly 图书《Python 数据分析》(第二版)。它并不太关注线性代数,但提供了关于在数据集上使用 NumPy、Pandas 和 Python 的实用指导。
练习
-
向量的值为[1, 2],然后发生了一个变换。落在[2, 0],落在[0, 1.5]。会落在哪里?
-
向量的值为[1, 2],然后发生了一个变换。落在[-2, 1],落在[1, -2]。会落在哪里?
-
一个变换落在[1, 0],落在[2, 2]。这个变换的行列式是多少?
-
一个线性变换中是否可以进行两个或更多个线性变换?为什么?
-
解方程组以求解x、y和z:
-
以下矩阵是否线性相关?为什么?
答案在附录 B 中。
第五章:线性回归
数据分析中最实用的技术之一是通过观察数据点拟合一条直线,以展示两个或更多变量之间的关系。回归试图将一个函数拟合到观察数据中,以对新数据进行预测。线性回归将一条直线拟合到观察数据中,试图展示变量之间的线性关系,并对尚未观察到的新数据进行预测。
看一张线性回归的图片可能比阅读描述更有意义。在图 5-1 中有一个线性回归的例子。
线性回归是数据科学和统计学的中流砥柱,不仅应用了我们在前几章学到的概念,还为后续主题如神经网络(第七章)和逻辑回归(第六章)奠定了新的基础。这种相对简单的技术已经存在了两百多年,当代被称为一种机器学习形式。
机器学习从业者通常采用不同的验证方法,从数据的训练-测试分割开始。统计学家更有可能使用像预测区间和相关性这样的指标来进行统计显著性分析。我们将涵盖这两种思维方式,以便读者能够弥合这两个学科之间日益扩大的鸿沟,从而最好地装备自己。
![emds 0501]()
图 5-1 线性回归的示例,将一条直线拟合到观察数据中
一个基本的线性回归
我想研究狗的年龄与它看兽医的次数之间的关系。在一个虚构的样本中,我们有 10 只随机的狗。我喜欢用简单的数据集(真实或其他)来理解复杂的技术,这样我们就能了解技术的优势和局限性,而不会被复杂的数据搞混。让我们将这个数据集绘制成图 5-2 所示。
![emds 0502]()
图 5-2 绘制了 10 只狗的样本,显示它们的年龄和看兽医的次数
我们可以清楚地看到这里存在着线性相关性,意味着当这些变量中的一个增加/减少时,另一个也以大致相同的比例增加/减少。我们可以在图 5-3 中画一条线来展示这样的相关性。
![emds 0503]()
图 5-3 拟合我们数据的一条线
我将在本章后面展示如何计算这条拟合线。我们还将探讨如何计算这条拟合线的质量。现在,让我们专注于执行线性回归的好处。它使我们能够对我们以前没有见过的数据进行预测。我的样本中没有一只 8.5 岁的狗,但我可以看着这条线估计这只狗一生中会有 21 次兽医就诊。我只需看看当x = 8.5 时,y = 21.218,如图 5-4 所示。另一个好处是我们可以分析可能存在关系的变量,并假设相关的变量之间是因果关系。
现在线性回归的缺点是什么?我不能期望每个结果都会完全落在那条线上。毕竟,现实世界的数据是嘈杂的,从不完美,也不会遵循一条直线。它可能根本不会遵循一条直线!在那条线周围会有误差,点会在线的上方或下方。当我们谈论 p 值、统计显著性和预测区间时,我们将在数学上涵盖这一点,这些内容描述了我们的线性回归有多可靠。另一个问题是我们不应该使用线性回归来预测超出我们拥有数据范围之外的情况,这意味着我们不应该在x < 0 和x > 10 的情况下进行预测,因为我们没有这些值之外的数据。
![emds 0504]()
图 5-4。使用线性回归进行预测,看到一个 8.5 岁的狗预测将有约 21.2 次兽医就诊
不要忘记抽样偏差!
我们应该质疑这些数据以及它们是如何抽样的,以便检测偏见。这是在单个兽医诊所吗?多个随机诊所?通过使用兽医数据是否存在自我选择偏见,只调查拜访兽医的狗?如果这些狗是在相同的地理位置抽样的,那会不会影响数据?也许在炎热的沙漠气候中的狗更容易因中暑和被蛇咬而去看兽医,这会使我们样本中的兽医就诊次数增加。
正如在第三章中讨论的那样,将数据视为真理的神谕已经变得时髦。然而,数据只是从一个总体中抽取的样本,我们需要对我们的样本有多好地代表性进行判断。对数据的来源同样感兴趣(如果不是更多),而不仅仅是数据所说的内容。
使用 SciPy 进行基本线性回归
在本章中,我们有很多关于线性回归要学习的内容,但让我们从一些代码开始执行我们已经了解的内容。
有很多平台可以执行线性回归,从 Excel 到 Python 和 R。但在本书中,我们将坚持使用 Python,从 scikit-learn 开始为我们完成工作。我将在本章后面展示如何“从头开始”构建线性回归,以便我们掌握像梯度下降和最小二乘这样的重要概念。
示例 5-1 是我们如何使用 scikit-learn 对这 10 只狗进行基本的、未经验证的线性回归的样本。我们使用 Pandas 获取这些数据,将其转换为 NumPy 数组,使用 scikit-learn 进行线性回归,并使用 Plotly 在图表中显示它。
示例 5-1. 使用 scikit-learn 进行线性回归
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
# Import points
df = pd.read_csv('https://bit.ly/3goOAnt', delimiter=",")
# Extract input variables (all rows, all columns but last column)
X = df.values[:, :-1]
# Extract output column (all rows, last column)
Y = df.values[:, -1]
# Fit a line to the points
fit = LinearRegression().fit(X, Y)
# m = 1.7867224, b = -16.51923513
m = fit.coef_.flatten()
b = fit.intercept_.flatten()
print("m = {0}".format(m))
print("b = {0}".format(b))
# show in chart
plt.plot(X, Y, 'o') # scatterplot
plt.plot(X, m*X+b) # line
plt.show()
首先,我们从GitHub 上的这个 CSV导入数据。我们使用 Pandas 将两列分离为X和Y数据集。然后,我们将LinearRegression
模型拟合到输入的X数据和输出的Y数据。然后我们可以得到描述我们拟合线性函数的m和b系数。
在图中,您将确实看到一条拟合线穿过这些点,如图 5-5 所示。
![emds 0505]()
图 5-5. SciPy 将拟合一条回归线到您的数据
是什么决定了最佳拟合线到这些点?让我们接下来讨论这个问题。
残差和平方误差
统计工具如 scikit-learn 如何得出适合这些点的线?这归结为机器学习训练中的两个基本问题:
-
什么定义了“最佳拟合”?
-
我们如何得到那个“最佳拟合”呢?
第一个问题有一个相当确定的答案:我们最小化平方,或更具体地说是平方残差的和。让我们来详细解释一下。画出任何一条穿过点的线。残差是线和点之间的数值差异,如图 5-6 所示。
![emds 0506]()
图 5-6. 残差是线和点之间的差异
线上方的点将具有正残差,而线下方的点将具有负残差。换句话说,它是预测 y 值(来自线)与实际 y 值(来自数据)之间的减去差异。残差的另一个名称是误差,因为它反映了我们的线在预测数据方面有多么错误。
让我们计算这 10 个点与线y = 1.93939x + 4.73333 之间的差异,以及示例 5-2 中每个点的残差和示例 5-3。
示例 5-2. 计算给定线和数据的残差
import pandas as pd
# Import points
points = pd.read_csv('https://bit.ly/3goOAnt', delimiter=",").itertuples()
# Test with a given line
m = 1.93939
b = 4.73333
# Calculate the residuals
for p in points:
y_actual = p.y
y_predict = m*p.x + b
residual = y_actual - y_predict
print(residual)
示例 5-3. 每个点的残差
-1.67272
1.3878900000000005
-0.5515000000000008
2.5091099999999997
-0.4302799999999998
-1.3696699999999993f
0.6909400000000012
-2.2484499999999983
2.812160000000002
-1.1272299999999973
如果我们要通过我们的 10 个数据点拟合一条直线,我们很可能希望尽可能地减小这些残差,使线和点之间的间隙尽可能小。但是我们如何衡量“总体”呢?最好的方法是采用平方和,简单地对每个残差进行平方,或者将每个残差相乘,然后将它们求和。我们取每个实际的 y 值,并从中减去从线上取得的预测 y 值,然后对所有这些差异进行平方和。
一个直观的思考方式如图 5-7 所示,我们在每个残差上叠加一个正方形,每条边的长度都是残差。我们将所有这些正方形的面积相加,稍后我们将学习如何通过确定最佳m和b来找到我们可以实现的最小和。
![emds 0507]()
图 5-7. 可视化平方和,即所有正方形的面积之和,其中每个正方形的边长等于残差
让我们修改我们在示例 5-4 中的代码来找到平方和。
示例 5-4. 计算给定直线和数据的平方和
import pandas as pd
# Import points
points = pd.read_csv("https://bit.ly/2KF29Bd").itertuples()
# Test with a given line
m = 1.93939
b = 4.73333
sum_of_squares = 0.0
# calculate sum of squares
for p in points:
y_actual = p.y
y_predict = m*p.x + b
residual_squared = (y_predict - y_actual)**2
sum_of_squares += residual_squared
print("sum of squares = {}".format(sum_of_squares))
# sum of squares = 28.096969704500005
下一个问题是:如何找到能产生最小平方和的m和b值,而不使用像 scikit-learn 这样的库?让我们接着看。
寻找最佳拟合直线
现在我们有一种方法来衡量给定直线与数据点的质量:平方和。我们能够使这个数字越低,拟合就越好。那么如何找到能产生最小平方和的正确m和b值呢?
我们可以采用几种搜索算法,试图找到解决给定问题的正确值集。你可以尝试蛮力方法,随机生成m和b值数百万次,并选择产生最小平方和的值。这种方法效果不佳,因为即使找到一个体面的近似值也需要无尽的时间。我们需要一些更有指导性的东西。我将为你整理五种技术:闭式方程、矩阵求逆、矩阵分解、梯度下降和随机梯度下降。还有其他搜索算法,比如爬山算法(在附录 A 中有介绍),但我们将坚持使用常见的方法。
闭式方程
一些读者可能会问是否有一个公式(称为闭式方程)通过精确计算来拟合线性回归。答案是肯定的,但仅适用于只有一个输入变量的简单线性回归。对于具有多个输入变量和大量数据的许多机器学习问题,这种奢侈是不存在的。我们可以使用线性代数技术进行扩展,我们很快将讨论这一点。我们还将借此机会学习诸如随机梯度下降之类的搜索算法。
对于只有一个输入和一个输出变量的简单线性回归,以下是计算m和b的闭式方程。示例 5-5 展示了如何在 Python 中进行这些计算。
示例 5-5. 计算简单线性回归的m和b
import pandas as pd
# Load the data
points = list(pd.read_csv('https://bit.ly/2KF29Bd', delimiter=",").itertuples())
n = len(points)
m = (n*sum(p.x*p.y for p in points) - sum(p.x for p in points) *
sum(p.y for p in points)) / (n*sum(p.x**2 for p in points) -
sum(p.x for p in points)**2)
b = (sum(p.y for p in points) / n) - m * sum(p.x for p in points) / n
print(m, b)
# 1.9393939393939394 4.7333333333333325
这些用于计算m和b的方程式是从微积分中推导出来的,如果您有兴趣发现公式的来源,我们稍后在本章中将使用 SymPy 进行一些微积分工作。目前,您可以插入数据点数n,并迭代 x 和 y 值来执行刚才描述的操作。
今后,我们将学习更适用于处理大量数据的现代技术。闭式方程式往往不适用于大规模应用。
计算复杂度
闭式方程式不适用于较大数据集的原因是由于计算机科学中称为计算复杂度的概念,它衡量算法在问题规模增长时所需的时间。这可能值得熟悉一下;以下是关于这个主题的两个很棒的 YouTube 视频:
逆矩阵技术
今后,我有时会用不同的名称交替使用系数m和b,分别为和,这是您在专业世界中经常看到的惯例,所以现在可能是一个毕业的好时机。
虽然我们在第四章中专门致力于线性代数,但在您刚接触数学和数据科学时,应用它可能有点令人不知所措。这就是为什么本书中的大多数示例将使用纯 Python 或 scikit-learn。但是,在合适的情况下,我会加入线性代数,以展示线性代数的实用性。如果您觉得这一部分令人不知所措,请随时继续阅读本章的其余部分,稍后再回来。
我们可以使用转置和逆矩阵,我们在第四章中介绍过,来拟合线性回归。接下来,我们根据输入变量值矩阵和输出变量值向量计算系数向量。在不深入微积分和线性代数证明的兔子洞中,这是公式:
你会注意到在矩阵上执行了转置和逆操作,并与矩阵乘法结合。这是我们在 NumPy 中执行此操作的方式,在例子 5-6 中得到我们的系数m和b。
例 5-6. 使用逆矩阵和转置矩阵拟合线性回归
import pandas as pd
from numpy.linalg import inv
import numpy as np
# Import points
df = pd.read_csv('https://bit.ly/3goOAnt', delimiter=",")
# Extract input variables (all rows, all columns but last column)
X = df.values[:, :-1].flatten()
# Add placeholder "1" column to generate intercept
X_1 = np.vstack([X, np.ones(len(X))]).T
# Extract output column (all rows, last column)
Y = df.values[:, -1]
# Calculate coefficents for slope and intercept
b = inv(X_1.transpose() @ X_1) @ (X_1.transpose() @ Y)
print(b) # [1.93939394, 4.73333333]
# Predict against the y-values
y_predict = X_1.dot(b)
这并不直观,但请注意我们必须在X列旁边堆叠一个“列”为 1 的列。原因是这将生成截距系数。由于这一列全为 1,它实际上生成了截距而不仅仅是一个斜率。
当你有大量数据和大量维度时,计算机可能开始崩溃并产生不稳定的结果。这是矩阵分解的一个用例,我们在线性代数的第四章中学到了。在这种特定情况下,我们取我们的矩阵X,附加一个额外的列 1 来生成截距就像以前一样,然后将其分解为两个组件矩阵Q和R:
避免更多的微积分兔子洞,这里是我们如何使用Q和R来在矩阵形式b中找到 beta 系数值:
而例子 5-7 展示了我们如何在 Python 中使用 NumPy 使用前述QR分解公式执行线性回归。
例 5-7. 使用 QR 分解执行线性回归
import pandas as pd
from numpy.linalg import qr, inv
import numpy as np
# Import points
df = pd.read_csv('https://bit.ly/3goOAnt', delimiter=",")
# Extract input variables (all rows, all columns but last column)
X = df.values[:, :-1].flatten()
# Add placeholder "1" column to generate intercept
X_1 = np.vstack([X, np.ones(len(X))]).transpose()
# Extract output column (all rows, last column)
Y = df.values[:, -1]
# calculate coefficents for slope and intercept
# using QR decomposition
Q, R = qr(X_1)
b = inv(R).dot(Q.transpose()).dot(Y)
print(b) # [1.93939394, 4.73333333]
通常,QR分解是许多科学库用于线性回归的方法,因为它更容易处理大量数据,并且更稳定。我所说的稳定是什么意思?数值稳定性是算法保持错误最小化的能力,而不是在近似中放大错误。请记住,计算机只能工作到某个小数位数,并且必须进行近似,因此我们的算法不应随着这些近似中的复合错误而恶化变得重要。
感到不知所措吗?
如果你觉得这些线性代数示例中的线性回归让人不知所措,不要担心!我只是想提供一个线性代数实际用例的曝光。接下来,我们将专注于其他你可以使用的技术。
梯度下降
梯度下降是一种优化技术,利用导数和迭代来最小化/最大化一组参数以达到目标。要了解梯度下降,让我们进行一个快速的思想实验,然后在一个简单的例子中应用它。
关于梯度下降的思想实验
想象一下,你在一个山脉中夜晚拿着手电筒。你试图到达山脉的最低点。在你迈出一步之前,你可以看到你周围的斜坡。你朝着斜坡明显向下的方向迈步。对于更大的斜坡,你迈出更大的步伐,对于更小的斜坡,你迈出更小的步伐。最终,你会发现自己在一个斜率为 0 的低点,一个值为 0。听起来不错,对吧?这种使用手电筒的方法被称为梯度下降,我们朝着斜坡向下的方向迈步。
在机器学习中,我们经常将我们将遇到的所有可能的平方损失总和视为多山的地形。我们想要最小化我们的损失,并且我们通过导航损失地形来实现这一点。为了解决这个问题,梯度下降有一个吸引人的特点:偏导数就像是那盏手电筒,让我们能够看到每个参数(在这种情况下是m和b,或者和)的斜率。我们朝着m和b的斜率向下的方向迈步。对于更大的斜率,我们迈出更大的步伐,对于更小的斜率,我们迈出更小的步伐。我们可以通过取斜率的一部分来简单地计算这一步的长度。这一部分被称为我们的学习率。学习率越高,它运行得越快,但精度会受到影响。但学习率越低,训练所需的时间就越长,需要更多的迭代。
决定学习率就像在选择蚂蚁、人类或巨人来踏下斜坡。蚂蚁(小学习率)会迈出微小的步伐,花费不可接受的长时间才能到达底部,但会准确无误地到达。巨人(大学习率)可能会一直跨过最小值,以至于无论走多少步都可能永远无法到达。人类(适度学习率)可能具有最平衡的步幅,在速度和准确性之间找到正确的平衡,以到达最小值。
先学会走再学会跑
对于函数,让我们找到产生该函数最低点的 x 值。虽然我们可以通过代数方法解决这个问题,但让我们使用梯度下降来做。
这是我们试图做的可视化效果。如图 5-8 所示,我们希望“步进”x朝向斜率为 0 的最小值。
![emds 0508]()
图 5-8. 朝向斜率接近 0 的局部最小值迈进
在示例 5-8 中,函数f(x)
及其对x的导数为dx_f(x)
。回想一下,我们在第一章中讨论了如何使用 SymPy 计算导数。找到导数后,我们继续执行梯度下降。
示例 5-8. 使用梯度下降找到抛物线的最小值
import random
def f(x):
return (x - 3) ** 2 + 4
def dx_f(x):
return 2*(x - 3)
# The learning rate
L = 0.001
# The number of iterations to perform gradient descent
iterations = 100_000
# start at a random x
x = random.randint(-15,15)
for i in range(iterations):
# get slope
d_x = dx_f(x)
# update x by subtracting the (learning rate) * (slope)
x -= L * d_x
print(x, f(x)) # prints 2.999999999999889 4.0
如果我们绘制函数(如图 5-8 所示),我们应该看到函数的最低点明显在x = 3 处,前面的代码应该非常接近这个点。学习率用于在每次迭代中取斜率的一部分并从 x 值中减去它。较大的斜率将导致较大的步长,而较小的斜率将导致较小的步长。经过足够的迭代,x将最终到达函数的最低点(或足够接近),其中斜率为 0。
梯度下降和线性回归
现在你可能想知道我们如何将其用于线性回归。嗯,这个想法是一样的,只是我们的“变量”是m和b(或和)而不是x。原因在于:在简单线性回归中,我们已经知道 x 和 y 值,因为这些值作为训练数据提供。我们需要解决的“变量”实际上是参数m和b,因此我们可以找到最佳拟合线,然后接受一个x变量来预测一个新的 y 值。
我们如何计算m和b的斜率?我们需要这两者的偏导数。我们要对哪个函数求导?记住我们试图最小化损失,这将是平方和。因此,我们需要找到我们的平方和函数对m和b的导数。
我像示例 5-9 中所示实现了这两个m和b的偏导数。我们很快将学习如何在 SymPy 中执行此操作。然后我执行梯度下降来找到m和b:100,000 次迭代,学习率为 0.001 就足够了。请注意,您将学习率设得越小,速度就越慢,需要的迭代次数就越多。但如果设得太高,它将运行得很快,但近似度较差。当有人说机器学习算法正在“学习”或“训练”时,实际上就是在拟合这样一个回归。
示例 5-9. 执行线性回归的梯度下降
import pandas as pd
# Import points from CSV
points = list(pd.read_csv("https://bit.ly/2KF29Bd").itertuples())
# Building the model
m = 0.0
b = 0.0
# The learning Rate
L = .001
# The number of iterations
iterations = 100_000
n = float(len(points)) # Number of elements in X
# Perform Gradient Descent
for i in range(iterations):
# slope with respect to m
D_m = sum(2 * p.x * ((m * p.x + b) - p.y) for p in points)
# slope with respect to b
D_b = sum(2 * ((m * p.x + b) - p.y) for p in points)
# update m and b
m -= L * D_m
b -= L * D_b
print("y = {0}x + {1}".format(m, b))
# y = 1.9393939393939548x + 4.733333333333227
嗯,不错!那个近似值接近我们的闭式方程解。但有什么问题吗?仅仅因为我们通过最小化平方和找到了“最佳拟合直线”,这并不意味着我们的线性回归就很好。最小化平方和是否保证了一个很好的模型来进行预测?并不完全是这样。现在我向你展示了如何拟合线性回归,让我们退一步,重新审视全局,确定给定的线性回归是否是首选的预测方式。但在我们这样做之前,这里有一个展示 SymPy 解决方案的更多绕路。
使用 SymPy 进行线性回归的梯度下降
如果你想要得到这两个关于平方和函数的导数的 SymPy 代码,分别为 m 和 b,这里是 示例 5-10 中的代码。
示例 5-10. 计算 m 和 b 的偏导数
from sympy import *
m, b, i, n = symbols('m b i n')
x, y = symbols('x y', cls=Function)
sum_of_squares = Sum((m*x(i) + b - y(i)) ** 2, (i, 0, n))
d_m = diff(sum_of_squares, m)
d_b = diff(sum_of_squares, b)
print(d_m)
print(d_b)
# OUTPUTS
# Sum(2*(b + m*x(i) - y(i))*x(i), (i, 0, n))
# Sum(2*b + 2*m*x(i) - 2*y(i), (i, 0, n))
你会看到分别为 m 和 b 的两个导数被打印出来。请注意 Sum()
函数将迭代并将项相加(在这种情况下是所有数据点),我们将 x 和 y 视为查找给定索引 i 处值的函数。
在数学符号中,其中 e(x) 代表平方和损失函数,这里是 m 和 b 的偏导数:
如果你想应用我们的数据集并使用梯度下降执行线性回归,你将需要执行一些额外的步骤,如示例 5-11 所示。我们需要替换n
、x(i)
和y(i)
的值,迭代所有数据点以计算d_m
和d_b
的导数函数。这样就只剩下m
和b
变量,我们将使用梯度下降寻找最优值。
示例 5-11. 使用 SymPy 解决线性回归
import pandas as pd
from sympy import *
# Import points from CSV
points = list(pd.read_csv("https://bit.ly/2KF29Bd").itertuples())
m, b, i, n = symbols('m b i n')
x, y = symbols('x y', cls=Function)
sum_of_squares = Sum((m*x(i) + b - y(i)) ** 2, (i, 0, n))
d_m = diff(sum_of_squares, m) \
.subs(n, len(points) - 1).doit() \
.replace(x, lambda i: points[i].x) \
.replace(y, lambda i: points[i].y)
d_b = diff(sum_of_squares, b) \
.subs(n, len(points) - 1).doit() \
.replace(x, lambda i: points[i].x) \
.replace(y, lambda i: points[i].y)
# compile using lambdify for faster computation
d_m = lambdify([m, b], d_m)
d_b = lambdify([m, b], d_b)
# Building the model
m = 0.0
b = 0.0
# The learning Rate
L = .001
# The number of iterations
iterations = 100_000
# Perform Gradient Descent
for i in range(iterations):
# update m and b
m -= d_m(m,b) * L
b -= d_b(m,b) * L
print("y = {0}x + {1}".format(m, b))
# y = 1.939393939393954x + 4.733333333333231
如示例 5-11 所示,对我们的偏导数函数都调用lambdify()
是个好主意,将它们从 SymPy 转换为优化的 Python 函数。这将使计算在执行梯度下降时更快。生成的 Python 函数由 NumPy、SciPy 或 SymPy 检测到的其他数值库支持。之后,我们可以执行梯度下降。
最后,如果你对这个简单线性回归的损失函数感兴趣,示例 5-12 展示了 SymPy 代码,将x
、y
和n
的值代入我们的损失函数,然后将m
和b
作为输入变量绘制出来。我们的梯度下降算法将我们带到了损失景观中的最低点,如图 5-9 所示。
示例 5-12. 绘制线性回归的损失函数
from sympy import *
from sympy.plotting import plot3d
import pandas as pd
points = list(pd.read_csv("https://bit.ly/2KF29Bd").itertuples())
m, b, i, n = symbols('m b i n')
x, y = symbols('x y', cls=Function)
sum_of_squares = Sum((m*x(i) + b - y(i)) ** 2, (i, 0, n)) \
.subs(n, len(points) - 1).doit() \
.replace(x, lambda i: points[i].x) \
.replace(y, lambda i: points[i].y)
plot3d(sum_of_squares)
![emds 0509]()
图 5-9. 简单线性回归的损失景观
过拟合和方差
你猜猜看:如果我们真的想最小化损失,即将平方和减少到 0,我们会怎么做?除了线性回归还有其他选择吗?你可能得出的一个结论就是简单地拟合一个触及所有点的曲线。嘿,为什么不只是连接点并用它来做预测,如图 5-10 所示?这样就得到了 0 的损失!
真糟糕,为什么我们要费力进行线性回归而不是做这个呢?嗯,记住我们的大局目标不是最小化平方和,而是在新数据上做出准确的预测。这种连接点模型严重过拟合,意味着它将回归形状调整得太精确到预测新数据时表现糟糕。这种简单的连接点模型对远离其他点的异常值敏感,意味着它在预测中具有很高的方差。虽然这个例子中的点相对接近一条直线,但在其他具有更广泛分布和异常值的数据集中,这个问题会更严重。因为过拟合增加了方差,预测结果将会到处都是!
![emds 0510]()
图 5-10. 通过简单连接点执行回归,导致损失为零
过拟合就是记忆
当有人说回归“记住”了数据而不是泛化它时,他们在谈论过拟合。
正如你所猜测的,我们希望在模型中找到有效的泛化,而不是记忆数据。否则,我们的回归模型简单地变成了一个数据库,我们只是查找数值。
在机器学习中,你会发现模型中添加了偏差,而线性回归被认为是一个高度偏置的模型。这与数据中的偏差不同,我们在第三章中有详细讨论。模型中的偏差意味着我们优先考虑一种方法(例如,保持一条直线),而不是弯曲和完全适应数据。一个有偏差的模型留有一些余地,希望在新数据上最小化损失以获得更好的预测,而不是在训练数据上最小化损失。我想你可以说,向模型添加偏差可以抵消过拟合,或者说对训练数据拟合较少。
你可以想象,这是一个平衡的过程,因为这是两个相互矛盾的目标。在机器学习中,我们基本上是在说,“我想要将回归拟合到我的数据,但我不想拟合得太多。我需要一些余地来预测新数据的不同之处。”
套索回归和岭回归
线性回归的两个比较流行的变体是套索回归和岭回归。岭回归在线性回归中添加了进一步的偏差,以一种惩罚的形式,因此导致它对数据拟合较少。套索回归将尝试边缘化嘈杂的变量,这在你想要自动删除可能不相关的变量时非常有用。
然而,我们不能仅仅将线性回归应用于一些数据,进行一些预测,并假设一切都没问题。即使是一条直线的线性回归也可能过拟合。因此,我们需要检查和缓解过拟合和欠拟合,以找到两者之间的平衡点。除非根本没有平衡点,否则你应该完全放弃该模型。
随机梯度下降
在机器学习的背景下,你不太可能像之前那样在实践中进行梯度下降,我们在所有训练数据上进行训练(称为批量梯度下降)。在实践中,你更有可能执行随机梯度下降,它将在每次迭代中仅对数据集的一个样本进行训练。在小批量梯度下降中,会使用数据集的多个样本(例如,10 或 100 个数据点)进行每次迭代。
为什么每次迭代只使用部分数据?机器学习从业者引用了一些好处。首先,它显著减少了计算量,因为每次迭代不必遍历整个训练数据集,而只需部分数据。第二个好处是减少过拟合。每次迭代只暴露训练算法于部分数据,使损失景观不断变化,因此不会稳定在损失最小值。毕竟,最小化损失是导致过拟合的原因,因此我们引入一些随机性来创建一点欠拟合(但希望不要太多)。
当然,我们的近似变得松散,所以我们必须小心。这就是为什么我们很快会谈论训练/测试拆分,以及其他评估我们线性回归可靠性的指标。
示例 5-13 展示了如何在 Python 中执行随机梯度下降。如果将样本大小改为大于 1,它将执行小批量梯度下降。
示例 5-13。执行线性回归的随机梯度下降
import pandas as pd
import numpy as np
# Input data
data = pd.read_csv('https://bit.ly/2KF29Bd', header=0)
X = data.iloc[:, 0].values
Y = data.iloc[:, 1].values
n = data.shape[0] # rows
# Building the model
m = 0.0
b = 0.0
sample_size = 1 # sample size
L = .0001 # The learning Rate
epochs = 1_000_000 # The number of iterations to perform gradient descent
# Performing Stochastic Gradient Descent
for i in range(epochs):
idx = np.random.choice(n, sample_size, replace=False)
x_sample = X[idx]
y_sample = Y[idx]
# The current predicted value of Y
Y_pred = m * x_sample + b
# d/dm derivative of loss function
D_m = (-2 / sample_size) * sum(x_sample * (y_sample - Y_pred))
# d/db derivative of loss function
D_b = (-2 / sample_size) * sum(y_sample - Y_pred)
m = m - L * D_m # Update m
b = b - L * D_b # Update b
# print progress
if i % 10000 == 0:
print(i, m, b)
print("y = {0}x + {1}".format(m, b))
当我运行这个时,我得到了一个线性回归 y = 1.9382830354181135x + 4.753408787648379。显然,你的结果会有所不同,由于随机梯度下降,我们实际上不会收敛到特定的最小值,而是会停留在一个更广泛的邻域。
随机性是坏事吗?
如果这种随机性让你感到不舒服,每次运行一段代码都会得到不同的答案,那么欢迎来到机器学习、优化和随机算法的世界!许多进行近似的算法都是基于随机性的,虽然有些非常有用,但有些可能效果不佳,正如你所预期的那样。
很多人把机器学习和人工智能看作是一种能够给出客观和精确答案的工具,但事实并非如此。机器学习产生的是带有一定不确定性的近似值,通常在生产中没有基本事实。如果不了解它的工作原理,机器学习可能会被滥用,不承认其非确定性和近似性质是不妥的。
虽然随机性可以创造一些强大的工具,但也可能被滥用。要小心不要使用种子值和随机性来 p-hack 一个“好”结果,并努力分析你的数据和模型。
相关系数
看看这个散点图 图 5-11 以及它的线性回归。为什么线性回归在这里效果不太好?
![emds 0511]()
图 5-11。具有高方差的数据的散点图
这里的问题是数据具有很高的方差。如果数据极为分散,它将使方差增加到使预测变得不太准确和有用的程度,导致大的残差。当然,我们可以引入更偏向的模型,如线性回归,以不那么容易弯曲和响应方差。然而,欠拟合也会削弱我们的预测,因为数据如此分散。我们需要数值化地衡量我们的预测有多“偏离”。
那么如何对这些残差进行整体测量呢?你又如何了解数据中方差的糟糕程度呢?让我向你介绍相关系数,也称为皮尔逊相关系数,它以-1 到 1 之间的值来衡量两个变量之间关系的强度。相关系数越接近 0,表示没有相关性。相关系数越接近 1,表示强正相关,意味着一个变量增加时,另一个变量成比例增加。如果接近-1,则表示强负相关,这意味着一个变量增加时,另一个成比例减少。
请注意,相关系数通常表示为r。在图 5-11 中高度分散的数据具有相关系数 0.1201。由于它比 1 更接近 0,我们可以推断数据之间关系很小。
这里是另外四个散点图,显示它们的相关系数。请注意,点越接近一条线,相关性越强。点更分散会导致相关性较弱。
![emds 0512]()
图 5-12。四个散点图的相关系数
可以想象,相关系数对于查看两个变量之间是否存在可能的关系是有用的。如果存在强正负关系,它将对我们的线性回归有所帮助。如果没有关系,它们可能只会添加噪音并损害模型的准确性。
我们如何使用 Python 计算相关系数?让我们使用之前使用的简单10 点数据集。分析所有变量对之间的相关性的快速简单方法是使用 Pandas 的corr()
函数。这使得轻松查看数据集中每对变量之间的相关系数,这种情况下只会是x
和y
。这被称为相关矩阵。在示例 5-14 中查看。
示例 5-14。使用 Pandas 查看每对变量之间的相关系数
import pandas as pd
# Read data into Pandas dataframe
df = pd.read_csv('https://bit.ly/2KF29Bd', delimiter=",")
# Print correlations between variables
correlations = df.corr(method='pearson')
print(correlations)
# OUTPUT:
# x y
# x 1.000000 0.957586
# y 0.957586 1.000000
正如您所看到的,x
和y
之间的相关系数0.957586
表明这两个变量之间存在强烈的正相关性。您可以忽略矩阵中x
或y
设置为自身且值为1.0
的部分。显然,当x
或y
设置为自身时,相关性将完美地为 1.0,因为值与自身完全匹配。当您有两个以上的变量时,相关性矩阵将显示更大的网格,因为有更多的变量进行配对和比较。
如果您更改代码以使用具有大量变化的不同数据集,其中数据分散,您将看到相关系数下降。这再次表明了较弱的相关性。
统计显著性
这里还有线性回归的另一个方面需要考虑:我的数据相关性是否巧合?在第三章中,我们研究了假设检验和 p 值,我们将在这里用线性回归扩展这些想法。
让我们从一个基本问题开始:我是否可能由于随机机会在我的数据中看到线性关系?我们如何能够确信这两个变量之间的相关性是显著的而不是巧合的 95%?如果这听起来像第三章中的假设检验,那是因为它就是!我们不仅需要表达相关系数,还需要量化我们对相关系数不是偶然发生的信心。
与我们在第三章中使用药物测试示例中所做的估计均值不同,我们正在基于样本估计总体相关系数。我们用希腊字母符号(Rho)表示总体相关系数,而我们的样本相关系数是r。就像我们在第三章中所做的那样,我们将有一个零假设和备择假设:
我们的零假设是两个变量之间没有关系,或更技术性地说,相关系数为 0。备择假设是存在关系,可以是正相关或负相关。这就是为什么备择假设被定义为,以支持正相关和负相关。
让我们回到我们的包含 10 个点的数据集,如图 5-13 所示。我们看到这些数据点是多大概率是偶然看到的?它们恰好产生了看起来是线性关系?
![emds 0513]()
图 5-13. 这些数据看起来具有线性相关性,我们有多大可能性是随机机会看到的?
我们已经在示例 5-14 中计算了这个数据集的相关系数为 0.957586。这是一个强有力的正相关。但是,我们需要评估这是否是由于随机运气。让我们以 95%的置信度进行双尾检验,探讨这两个变量之间是否存在关系。
我们在第三章中讨论了 T 分布,它有更厚的尾部以捕捉更多的方差和不确定性。我们使用 T 分布而不是正态分布进行线性回归的假设检验。首先,让我们绘制一个 T 分布,95%的临界值范围如图 5-14 所示。考虑到我们的样本中有 10 条记录,因此我们有 9 个自由度(10-1=9)。
![emds 0514]()
图 5-14. 9 个自由度的 T 分布,因为有 10 条记录,我们减去 1
临界值约为±2.262,我们可以在 Python 中计算如示例 5-16 所示。这捕捉了我们 T 分布中心区域的 95%。
示例 5-16. 从 T 分布计算临界值
from scipy.stats import t
n = 10
lower_cv = t(n-1).ppf(.025)
upper_cv = t(n-1).ppf(.975)
print(lower_cv, upper_cv)
# -2.262157162740992 2.2621571627409915
如果我们的检验值恰好落在(-2.262,2.262)的范围之外,那么我们可以拒绝我们的零假设。要计算检验值t,我们需要使用以下公式。再次,r是相关系数,n是样本大小:
让我们在 Python 中将整个测试放在一起,如示例 5-17 所示。如果我们的检验值落在 95%置信度的临界范围之外,我们接受我们的相关性不是偶然的。
示例 5-17. 测试看起来线性的数据的显著性
from scipy.stats import t
from math import sqrt
# sample size
n = 10
lower_cv = t(n-1).ppf(.025)
upper_cv = t(n-1).ppf(.975)
# correlation coefficient
# derived from data https://bit.ly/2KF29Bd
r = 0.957586
# Perform the test
test_value = r / sqrt((1-r**2) / (n-2))
print("TEST VALUE: {}".format(test_value))
print("CRITICAL RANGE: {}, {}".format(lower_cv, upper_cv))
if test_value < lower_cv or test_value > upper_cv:
print("CORRELATION PROVEN, REJECT H0")
else:
print("CORRELATION NOT PROVEN, FAILED TO REJECT H0 ")
# Calculate p-value
if test_value > 0:
p_value = 1.0 - t(n-1).cdf(test_value)
else:
p_value = t(n-1).cdf(test_value)
# Two-tailed, so multiply by 2
p_value = p_value * 2
print("P-VALUE: {}".format(p_value))
这里的检验值约为 9.39956,明显超出了(-2.262,2.262)的范围,因此我们可以拒绝零假设,并说我们的相关性是真实的。这是因为 p 值非常显著:0.000005976。这远低于我们的 0.05 阈值,因此这几乎不是巧合:存在相关性。p 值如此之小是有道理的,因为这些点强烈地类似于一条线。这些点随机地如此靠近一条线的可能性极小。
图 5-15 展示了一些其他数据集及其相关系数和 p 值。分析每一个。哪一个可能对预测最有用?其他数据集存在什么问题?
![emds 0515]()
图 5-15。不同数据集及其相关系数和 p 值
现在你有机会对来自图 5-15 的数据集进行分析后,让我们来看看结果。左侧图具有很高的正相关性,但只有三个数据点。数据不足显著提高了 p 值,达到 0.34913,并增加了数据发生偶然性的可能性。这是有道理的,因为只有三个数据点很可能会看到一个线性模式,但这并不比只有两个点好,这两个点只会连接一条直线。这提出了一个重要的规则:拥有更多数据将降低你的 p 值,特别是如果这些数据趋向于一条线。
第二幅图就是我们刚刚讨论的内容。它只有 10 个数据点,但形成了一个线性模式,我们不仅有很强的正相关性,而且 p 值极低。当 p 值如此之低时,你可以确定你正在测量一个经过精心设计和严格控制的过程,而不是某种社会学或自然现象。
图 5-15 中右侧的两幅图未能确定线性关系。它们的相关系数接近于 0,表明没有相关性,而 p 值不出所料地表明随机性起了作用。
规则是这样的:拥有更多数据且一致地类似于一条线,你的相关性的 p 值就会更显著。数据越分散或稀疏,p 值就会增加,从而表明你的相关性是由随机机会引起的。
决定系数
让我们学习一个在统计学和机器学习回归中经常见到的重要指标。决定系数,称为,衡量一个变量的变异有多少是由另一个变量的变异解释的。它也是相关系数的平方。当接近完美相关(-1 或 1)时,接近 1。基本上,显示了两个变量相互作用的程度。
让我们继续查看我们从图 5-13 中的数据。在示例 5-18 中,使用我们之前计算相关系数的数据框代码,然后简单地对其进行平方。这将使每个相关系数相互乘以自己。
示例 5-18。在 Pandas 中创建相关性矩阵
import pandas as pd
# Read data into Pandas dataframe
df = pd.read_csv('https://bit.ly/2KF29Bd', delimiter=",")
# Print correlations between variables
coeff_determination = df.corr(method='pearson') ** 2
print(coeff_determination)
# OUTPUT:
# x y
# x 1.000000 0.916971
# y 0.916971 1.000000
决定系数为 0.916971 被解释为x的变异的 91.6971%由y(反之亦然)解释,剩下的 8.3029%是由其他未捕获的变量引起的噪音;0.916971 是一个相当不错的决定系数,显示x和y解释彼此的方差。但可能有其他变量在起作用,占据了剩下的 0.083029。记住,相关性不等于因果关系,因此可能有其他变量导致我们看到的关系。
相关性不代表因果关系!
需要注意的是,虽然我们非常强调测量相关性并围绕其构建指标,请记住相关性不代表因果关系!你可能以前听过这句口头禅,但我想扩展一下统计学家为什么这么说。
仅仅因为我们看到x和y之间的相关性,并不意味着x导致y。实际上可能是y导致x!或者可能存在第三个未捕获的变量z导致x和y。也可能x和y根本不相互导致,相关性只是巧合,因此我们测量统计显著性非常重要。
现在我有一个更紧迫的问题要问你。计算机能区分相关性和因果关系吗?答案是“绝对不!”计算机有相关性的概念,但没有因果关系。假设我加载一个数据集到 scikit-learn,显示消耗的水量和我的水费。我的计算机,或包括 scikit-learn 在内的任何程序,都不知道更多的用水量是否导致更高的账单,或更高的账单是否导致更多的用水量。人工智能系统很容易得出后者的结论,尽管这是荒谬的。这就是为什么许多机器学习项目需要一个人来注入常识。
在计算机视觉中,这也会发生。计算机视觉通常会对数字像素进行回归以预测一个类别。如果我训练一个计算机视觉系统来识别牛,使用的是牛的图片,它可能会轻易地将场地与牛进行关联。因此,如果我展示一张空旷的场地的图片,它会将草地标记为牛!这同样是因为计算机没有因果关系的概念(牛的形状应该导致标签“牛”),而是陷入了我们不感兴趣的相关性中。
估计的标准误差
衡量线性回归整体误差的一种方法是SSE,或者平方误差和。我们之前学过这个概念,其中我们对每个残差进行平方并求和。如果(读作“y-hat”)是线上的每个预测值,而代表数据中的每个实际 y 值,这里是计算公式:
然而,所有这些平方值很难解释,所以我们可以使用一些平方根逻辑将事物重新缩放到它们的原始单位。我们还将所有这些值求平均,这就是估计的标准误差()的作用。如果n是数据点的数量,示例 5-19 展示了我们如何在 Python 中计算标准误差。
示例 5-19。计算估计的标准误差
Here is how we calculate it in Python:
import pandas as pd
from math import sqrt
# Load the data
points = list(pd.read_csv('https://bit.ly/2KF29Bd', delimiter=",").itertuples())
n = len(points)
# Regression line
m = 1.939
b = 4.733
# Calculate Standard Error of Estimate
S_e = sqrt((sum((p.y - (m*p.x +b))**2 for p in points))/(n-2))
print(S_e)
# 1.87406793500129
为什么是而不是像我们在第三章中的许多方差计算中所做的?不深入数学证明,这是因为线性回归有两个变量,而不只是一个,所以我们必须在自由度中再增加一个不确定性。
你会注意到估计的标准误差看起来与我们在第三章中学习的标准差非常相似。这并非偶然。这是因为它是线性回归的标准差。
预测区间
正如前面提到的,线性回归中的数据是从一个总体中取样得到的。因此,我们的回归结果只能和我们的样本一样好。我们的线性回归线也沿着正态分布运行。实际上,这使得每个预测的 y 值都像均值一样是一个样本统计量。事实上,“均值”沿着这条线移动。
还记得我们在第二章中讨论方差和标准差吗?这些概念在这里也适用。通过线性回归,我们希望数据以线性方式遵循正态分布。回归线充当我们钟形曲线的“均值”,数据围绕该线的分布反映了方差/标准差,如图 5-16 所示。
![emds 0516]()
图 5-16。线性回归假设正态分布遵循该线
当我们有一个正态分布遵循线性回归线时,我们不仅有一个变量,还有第二个变量引导着分布。每个y预测周围都有一个置信区间,这被称为预测区间。
让我们通过兽医示例重新带回一些背景,估计一只狗的年龄和兽医访问次数。我想知道一个 8.5 岁狗的兽医访问次数的 95%置信度的预测区间。这个预测区间的样子如图 5-17 所示。我们有 95%的信心,一个 8.5 岁的狗将有 16.462 到 25.966 次兽医访问。
![emds 0517]()
图 5-17。一个 8.5 岁狗的 95%置信度的预测区间
我们如何计算这个?我们需要得到误差边界,并在预测的 y 值周围加减。这是一个涉及 T 分布的临界值以及估计标准误差的庞大方程。让我们来看一下:
我们感兴趣的 x 值被指定为,在这种情况下是 8.5。这是我们如何在 Python 中解决这个问题的,如示例 5-20 所示。
示例 5-20。计算一只 8.5 岁狗的兽医访问预测区间
import pandas as pd
from scipy.stats import t
from math import sqrt
# Load the data
points = list(pd.read_csv('https://bit.ly/2KF29Bd', delimiter=",").itertuples())
n = len(points)
# Linear Regression Line
m = 1.939
b = 4.733
# Calculate Prediction Interval for x = 8.5
x_0 = 8.5
x_mean = sum(p.x for p in points) / len(points)
t_value = t(n - 2).ppf(.975)
standard_error = sqrt(sum((p.y - (m * p.x + b)) ** 2 for p in points) / (n - 2))
margin_of_error = t_value * standard_error * \
sqrt(1 + (1 / n) + (n * (x_0 - x_mean) ** 2) / \
(n * sum(p.x ** 2 for p in points) - \
sum(p.x for p in points) ** 2))
predicted_y = m*x_0 + b
# Calculate prediction interval
print(predicted_y - margin_of_error, predicted_y + margin_of_error)
# 16.462516875955465 25.966483124044537
哎呀!这是很多计算,不幸的是,SciPy 和其他主流数据科学库都不会做这个。但如果你倾向于统计分析,这是非常有用的信息。我们不仅基于线性回归创建预测(例如,一只 8.5 岁的狗将有 21.2145 次兽医访问),而且实际上能够说出一些远非绝对的东西:一个 8.5 岁的狗会在 16.46 到 25.96 次之间访问兽医的概率为 95%。很棒,对吧?这是一个更安全的说法,因为它涵盖了一个范围而不是一个单一值,因此考虑了不确定性。
训练/测试分割
我刚刚进行的这个分析,包括相关系数、统计显著性和决定系数,不幸的是并不总是由从业者完成。有时候他们处理的数据太多,没有时间或技术能力这样做。例如,一个 128×128 像素的图像至少有 16,384 个变量。你有时间对每个像素变量进行统计分析吗?可能没有!不幸的是,这导致许多数据科学家根本不学习这些统计指标。
在一个不知名的在线论坛上,我曾经看到一篇帖子说统计回归是手术刀,而机器学习是电锯。当处理大量数据和变量时,你无法用手术刀筛选所有这些。你必须求助于电锯,虽然你会失去可解释性和精度,但至少可以扩展到更多数据上进行更广泛的预测。话虽如此,抽样偏差和过拟合等统计问题并没有消失。但有一些实践方法可以用于快速验证。
为什么 scikit-learn 中没有置信区间和 P 值?
Scikit-learn 不支持置信区间和 P 值,因为这两种技术对于高维数据是一个悬而未决的问题。这只强调了统计学家和机器学习从业者之间的差距。正如 scikit-learn 的一位维护者 Gael Varoquaux 所说,“通常计算正确的 P 值需要对数据做出假设,而这些假设并不符合机器学习中使用的数据(没有多重共线性,与维度相比有足够的数据)....P 值是一种期望得到很好检查的东西(在医学研究中是一种保护)。实现它们会带来麻烦....我们只能在非常狭窄的情况下给出 P 值[有少量变量]。”
如果你想深入了解,GitHub 上有一些有趣的讨论:
如前所述,statsmodel 是一个为统计分析提供有用工具的库。只是要知道,由于前述原因,它可能不会适用于更大维度的模型。
机器学习从业者用于减少过拟合的基本技术之一是一种称为训练/测试分割的实践,通常将 1/3 的数据用于测试,另外的 2/3 用于训练(也可以使用其他比例)。训练数据集用于拟合线性回归,而测试数据集用于衡量线性回归在之前未见数据上的表现。这种技术通常用于所有监督学习,包括逻辑回归和神经网络。图 5-18 显示了我们如何将数据分割为 2/3 用于训练和 1/3 用于测试的可视化。
这是一个小数据集
正如我们将在后面学到的,有其他方法可以将训练/测试数据集分割为 2/3 和 1/3。如果你有一个这么小的数据集,你可能最好使用 9/10 和 1/10 与交叉验证配对,或者甚至只使用留一交叉验证。查看“训练/测试分割是否必须是三分之一?” 以了解更多。
![emds 0518]()
图 5-18. 将数据分割为训练/测试数据—使用最小二乘法将线拟合到训练数据(深蓝色),然后分析测试数据(浅红色)以查看预测在之前未见数据上的偏差
示例 5-21 展示了如何使用 scikit-learn 执行训练/测试分割,其中 1/3 的数据用于测试,另外的 2/3 用于训练。
训练即拟合回归
记住,“拟合”回归与“训练”是同义词。后者是机器学习从业者使用的词语。
示例 5-21. 在线性回归上进行训练/测试分割
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
# Load the data
df = pd.read_csv('https://bit.ly/3cIH97A', delimiter=",")
# Extract input variables (all rows, all columns but last column)
X = df.values[:, :-1]
# Extract output column (all rows, last column)
Y = df.values[:, -1]
# Separate training and testing data
# This leaves a third of the data out for testing
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=1/3)
model = LinearRegression()
model.fit(X_train, Y_train)
result = model.score(X_test, Y_test)
print("r²: %.3f" % result)
注意,train_test_split()
将获取我们的数据集(X 和 Y 列),对其进行洗牌,然后根据我们的测试数据集大小返回我们的训练和测试数据集。我们使用LinearRegression
的fit()
函数来拟合训练数据集X_train
和Y_train
。然后我们使用score()
函数在测试数据集X_test
和Y_test
上评估,从而让我们了解回归在之前未见过的数据上的表现。测试数据集的值越高,表示回归在之前未见过的数据上表现越好。具有更高数值的数字表示回归在之前未见过的数据上表现良好。
我们还可以在每个 1/3 折叠中交替使用测试数据集。这被称为交叉验证,通常被认为是验证技术的黄金标准。图 5-20 显示了数据的每个 1/3 轮流成为测试数据集。
![emds 0520]()
图 5-20. 三折交叉验证的可视化
示例 5-22 中的代码展示了跨三个折叠进行的交叉验证,然后评分指标(在本例中为均方和 [MSE])与其标准偏差一起平均,以展示每个测试的一致性表现。
示例 5-22. 使用三折交叉验证进行线性回归
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import KFold, cross_val_score
df = pd.read_csv('https://bit.ly/3cIH97A', delimiter=",")
# Extract input variables (all rows, all columns but last column)
X = df.values[:, :-1]
# Extract output column (all rows, last column)\
Y = df.values[:, -1]
# Perform a simple linear regression
kfold = KFold(n_splits=3, random_state=7, shuffle=True)
model = LinearRegression()
results = cross_val_score(model, X, Y, cv=kfold)
print(results)
print("MSE: mean=%.3f (stdev-%.3f)" % (results.mean(), results.std()))
当你开始关注模型中的方差时,你可以采用随机折叠验证,而不是简单的训练/测试拆分或交叉验证,重复地对数据进行洗牌和训练/测试拆分无限次,并汇总测试结果。在示例 5-23 中,有 10 次随机抽取数据的 1/3 进行测试,其余 2/3 进行训练。然后将这 10 个测试结果与它们的标准偏差平均,以查看测试数据集的表现一致性。
有什么问题?这在计算上非常昂贵,因为我们要多次训练回归。
示例 5-23. 使用随机折叠验证进行线性回归
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score, ShuffleSplit
df = pd.read_csv('https://bit.ly/38XwbeB', delimiter=",")
# Extract input variables (all rows, all columns but last column)
X = df.values[:, :-1]
# Extract output column (all rows, last column)\
Y = df.values[:, -1]
# Perform a simple linear regression
kfold = ShuffleSplit(n_splits=10, test_size=.33, random_state=7)
model = LinearRegression()
results = cross_val_score(model, X, Y, cv=kfold)
print(results)
print("mean=%.3f (stdev-%.3f)" % (results.mean(), results.std()))
因此,当你时间紧迫或数据量过大无法进行统计分析时,训练/测试拆分将提供一种衡量线性回归在未见过的数据上表现如何的方法。
训练/测试拆分并不保证结果
值得注意的是,仅仅因为你应用了机器学习的最佳实践,将训练和测试数据拆分,这并不意味着你的模型会表现良好。你很容易过度调整模型,并通过一些手段获得良好的测试结果,但最终发现在现实世界中并不奏效。这就是为什么有时候需要保留另一个数据集,称为验证集,特别是当你在比较不同模型或配置时。这样,你对训练数据的调整以获得更好的测试数据性能不会泄漏信息到训练中。你可以使用验证数据集作为最后一道防线,查看是否过度调整导致你对测试数据过拟合。
即使如此,你的整个数据集(包括训练、测试和验证)可能一开始就存在偏差,没有任何拆分可以减轻这种情况。Andrew Ng 在他与 DeepLearning.AI 和 Stanford HAI 的问答环节 中讨论了这个问题,他通过一个例子说明了为什么机器学习尚未取代放射科医生。
多元线性回归
在本章中,我们几乎完全专注于对一个输入变量和一个输出变量进行线性回归。然而,我们在这里学到的概念应该基本适用于多变量线性回归。像、标准误差和置信区间等指标可以使用,但随着变量的增加变得更加困难。示例 5-24 是一个使用 scikit-learn 进行的具有两个输入变量和一个输出变量的线性回归示例。
示例 5-24。具有两个输入变量的线性回归
import pandas as pd
from sklearn.linear_model import LinearRegression
# Load the data
df = pd.read_csv('https://bit.ly/2X1HWH7', delimiter=",")
# Extract input variables (all rows, all columns but last column)
X = df.values[:, :-1]
# Extract output column (all rows, last column)\
Y = df.values[:, -1]
# Training
fit = LinearRegression().fit(X, Y)
# Print coefficients
print("Coefficients = {0}".format(fit.coef_))
print("Intercept = {0}".format(fit.intercept_))
print("z = {0} + {1}x + {2}y".format(fit.intercept_, fit.coef_[0], fit.coef_[1]))
当模型变得充斥着变量以至于开始失去可解释性时,就会出现一定程度的不稳定性,这时机器学习实践开始将模型视为黑匣子。希望你相信统计问题并没有消失,随着添加的变量越来越多,数据变得越来越稀疏。但是,如果你退后一步,使用相关矩阵分析每对变量之间的关系,并寻求理解每对变量是如何相互作用的,这将有助于你努力创建一个高效的机器学习模型。
结论
在这一章中我们涵盖了很多内容。我们试图超越对线性回归的肤浅理解,并且不仅仅将训练/测试分割作为我们唯一的验证方式。我想向你展示割刀(统计学)和电锯(机器学习),这样你可以判断哪种对于你遇到的问题更好。仅在线性回归中就有许多指标和分析方法可用,我们涵盖了其中一些以了解线性回归是否可靠用于预测。你可能会发现自己处于一种情况,要么做出广泛的近似回归,要么使用统计工具仔细分析和整理数据。你使用哪种方法取决于情况,如果你想了解更多关于 Python 可用的统计工具,请查看statsmodel 库。
在第六章中涵盖逻辑回归时,我们将重新审视和统计显著性。希望这一章能让你相信有方法可以有意义地分析数据,并且这种投资可以在成功的项目中产生差异。
练习
提供了一个包含两个变量x和y的数据集这里。
-
执行简单的线性回归,找到最小化损失(平方和)的m和b值。
-
计算这些数据的相关系数和统计显著性(95%置信度)。相关性是否有用?
-
如果我预测x=50,那么y的预测值的 95%预测区间是多少?
-
重新开始你的回归,并进行训练/测试分割。随意尝试交叉验证和随机折叠验证。线性回归在测试数据上表现良好且一致吗?为什么?
答案在附录 B 中。
第六章:逻辑回归和分类
在本章中,我们将介绍逻辑回归,一种根据一个或多个自变量预测结果概率的回归类型。这反过来可以用于分类,即预测类别而不是像线性回归那样预测实数。
我们并不总是对将变量表示为连续感兴趣,其中它们可以表示无限数量的实数十进制值。有些情况下,我们更希望变量是离散的,或者代表整数、布尔值(1/0,真/假)。逻辑回归是在一个离散的输出变量上进行训练的(二进制 1 或 0)或一个分类数字(整数)。它输出一个概率的连续变量,但可以通过阈值转换为离散值。
逻辑回归易于实现,并且相对抗干扰和其他数据挑战。许多机器学习问题最好通过逻辑回归来解决,提供比其他类型的监督式机器学习更实用和更高性能的解决方案。
就像我们在第五章中讨论线性回归时所做的那样,我们将尝试在统计学和机器学习之间找到平衡,使用两个学科的工具和分析。逻辑回归将整合我们从本书中学到的许多概念,从概率到线性回归。
理解逻辑回归
想象一下,发生了一起小型工业事故,你正在尝试了解化学物质暴露的影响。有 11 名患者暴露于不同小时数的化学物质中(请注意这是虚构数据)。一些患者出现了症状(值为 1),而另一些没有出现症状(值为 0)。让我们在图 6-1 中绘制它们,其中 x 轴是暴露的小时数,y 轴是他们是否出现了症状(1 或 0)。
![emds 0601]()
图 6-1。绘制患者在* x *小时暴露后是否出现症状(1)或未出现症状(0)
患者在多长时间后开始出现症状?很容易看到,几乎在四小时后,我们立即从患者不出现症状(0)转变为出现症状(1)。在图 6-2 中,我们看到相同的数据带有一个预测曲线。
![emds 0602]()
图 6-2。四小时后,我们看到患者开始出现症状的明显跳跃
对这个样本进行粗略分析,我们可以说,暴露时间少于四小时的患者几乎不可能出现症状,但暴露时间超过四小时的患者出现症状的概率为 100%。在这两组之间,大约在四小时左右立即跳跃到出现症状。
当然,在现实世界中,没有什么是如此清晰明了的。假设你收集了更多数据,在范围的中间有一些患者表现出症状与不表现症状的混合,如图 6-3 所示。
![emds 0603]()
图 6-3. 中间存在一些表现出症状(1)和不表现症状(0)的患者混合
解释这一点的方式是,随着每小时的暴露,患者表现出症状的概率逐渐增加。让我们用一个逻辑函数或一个 S 形曲线来可视化这一点,其中输出变量被挤压在 0 和 1 之间,如图 6-4 所示。
![emds 0604]()
图 6-4. 将逻辑函数拟合到数据上
由于中间点的重叠,当患者表现出症状时没有明显的分界线,而是从 0%概率逐渐过渡到 100%概率(0 和 1)。这个例子展示了逻辑回归如何产生一个曲线,表示属于真实类别(患者表现出症状)的概率在一个独立变量(暴露小时数)上。
我们可以重新利用逻辑回归,不仅预测给定输入变量的概率,还可以添加一个阈值来预测它是否属于该类别。例如,如果我得到一个新患者,并发现他们暴露了六个小时,我预测他们有 71.1%的机会表现出症状,如图 6-5 所示。如果我的阈值至少为 50%的概率表现出症状,我将简单地分类为患者将表现出症状。
![emds 0605]()
图 6-5. 我们可以预期一个暴露了六个小时的患者有 71.1%的可能表现出症状,因为这大于 50%的阈值,我们预测他们将表现出症状
进行逻辑回归
那么我们如何进行逻辑回归呢?让我们首先看看逻辑函数,并探索其背后的数学。
逻辑函数
逻辑函数是一个 S 形曲线(也称为sigmoid 曲线),对于给定的一组输入变量,产生一个在 0 和 1 之间的输出变量。因为输出变量在 0 和 1 之间,它可以用来表示概率。
这是一个输出一个输入变量x的概率y的逻辑函数:
请注意,这个公式使用了欧拉数,我们在第一章中讨论过。x变量是独立/输入变量。和是我们需要解决的系数。
和 被打包在一个类似于线性函数的指数中,你可能会记得它看起来与 或 相同。这并非巧合;逻辑回归实际上与线性回归有着密切的关系,我们将在本章后面讨论。实际上, 是截距(在简单线性回归中我们称之为), 是x的斜率(在简单线性回归中我们称之为)。指数中的这个线性函数被称为对数几率函数,但现在只需知道整个逻辑函数产生了我们需要在 x 值上输出移动概率的 S 形曲线。
要在 Python 中声明逻辑函数,使用math
包中的exp()
函数声明e指数,如示例 6-1 所示。
示例 6-1. Python 中用于一个自变量的逻辑函数
import math
def predict_probability(x, b0, b1):
p = 1.0 / (1.0 + math.exp(-(b0 + b1 * x)))
return p
让我们绘制一下看看它是什么样子,并假设Β[0] = –2.823 和 Β[1] = 0.62。我们将在示例 6-2 中使用 SymPy,输出图形显示在图 6-6 中。
示例 6-2. 使用 SymPy 绘制逻辑函数
from sympy import *
b0, b1, x = symbols('b0 b1 x')
p = 1.0 / (1.0 + exp(-(b0 + b1 * x)))
p = p.subs(b0,-2.823)
p = p.subs(b1, 0.620)
print(p)
plot(p)
![emds 0606]()
图 6-6. 一个逻辑函数
在一些教科书中,你可能会看到逻辑函数被这样声明:
不要为此担心,因为这是相同的函数,只是代数上表达不同。注意,像线性回归一样,我们也可以将逻辑回归扩展到多于一个输入变量(),如此公式所示。我们只需添加更多的 系数:
拟合逻辑曲线
如何将逻辑曲线拟合到给定的训练数据集?首先,数据可以包含任意混合的十进制、整数和二进制变量,但输出变量必须是二进制(0 或 1)。当我们实际进行预测时,输出变量将在 0 和 1 之间,类似于概率。
数据提供了我们的输入和输出变量值,但我们需要解出和系数以拟合我们的逻辑函数。回想一下我们在 Chapter 5 中如何使用最小二乘法。然而,在这里不适用这种方法。相反,我们使用最大似然估计,顾名思义,最大化给定逻辑曲线输出观测数据的可能性。
要计算最大似然估计,实际上没有像线性回归那样的封闭形式方程。我们仍然可以使用梯度下降,或者让一个库来为我们做这件事。让我们从库 SciPy 开始涵盖这两种方法。
使用 SciPy
SciPy 的好处在于,模型通常具有一套标准化的函数和 API,这意味着在许多情况下,您可以复制/粘贴您的代码,然后在模型之间重复使用它。在 Example 6-3 中,您将看到我们的患者数据上执行的逻辑回归。如果您将其与我们在 Chapter 5 中的线性回归代码进行比较,您将看到在导入、分离和拟合数据方面几乎完全相同的代码。主要区别在于我使用LogisticRegression()
作为我的模型,而不是LinearRegression()
。
示例 6-3。在 SciPy 中使用普通逻辑回归
import pandas as pd
from sklearn.linear_model import LogisticRegression
# Load the data
df = pd.read_csv('https://bit.ly/33ebs2R', delimiter=",")
# Extract input variables (all rows, all columns but last column)
X = df.values[:, :-1]
# Extract output column (all rows, last column)
Y = df.values[:, -1]
# Perform logistic regression
# Turn off penalty
model = LogisticRegression(penalty='none')
model.fit(X, Y)
# print beta1
print(model.coef_.flatten()) # 0.69267212
# print beta0
print(model.intercept_.flatten()) # -3.17576395
在 SciPy 中运行模型后,我得到一个逻辑回归,其中β[0] = –3.17576395,β[1] = 0.69267212。当我绘制这个图时,应该看起来很好,就像在 Figure 6-7 中显示的那样。
![emds 0607]()
图 6-7。绘制逻辑回归
这里有几件事情需要注意。当我创建LogisticRegression()
模型时,我没有指定penalty
参数,这会选择像l1
或l2
这样的正则化技术。虽然这超出了本书的范围,但我在以下注释“学习关于 SciPy 参数”中包含了简要见解,以便您手边有有用的参考资料。
最后,我将flatten()
系数和截距,这些系数和截距出来时是多维矩阵,但只有一个元素。Flattening意味着将一组数字的矩阵折叠成较小的维度,特别是当元素少于维度时。例如,我在这里使用flatten()
来将嵌套在二维矩阵中的单个数字提取出来作为单个值。然后我有了我的和系数。
学习关于 SciPy 参数
SciPy 在其回归和分类模型中提供了许多选项。不幸的是,由于这不是一本专门关注机器学习的书籍,没有足够的带宽或页面来覆盖它们。
然而,SciPy 文档写得很好,逻辑回归页面在这里找到。
如果很多术语都很陌生,比如正则化和l1
和l2
惩罚,还有其他很棒的 O’Reilly 书籍探讨这些主题。我发现其中一本更有帮助的书是由 Aurélien Géron 撰写的Hands-On Machine Learning with Scikit-Learn, Keras, and TensorFlow。
使用最大似然和梯度下降
正如我在整本书中所做的,我旨在提供关于从头开始构建技术的见解,即使库可以为我们完成。有几种方法可以自己拟合逻辑回归,但所有方法通常都转向最大似然估计(MLE)。MLE 最大化了给定逻辑曲线输出观测数据的可能性。这与平方和不同,但我们仍然可以应用梯度下降或随机梯度下降来解决它。
我会尽量简化数学术语,并尽量减少线性代数的内容。基本上,这个想法是找到使我们的逻辑曲线尽可能接近这些点的和系数,表明它最有可能产生这些点。如果你还记得第二章中我们学习概率时,我们通过将多个事件的概率(或可能性)相乘来组合它们。在这个应用中,我们正在计算我们会看到所有这些点的可能性,对于给定的逻辑回归曲线。
应用联合概率的概念,每个患者都有一个基于拟合的逻辑函数的可能性,如图 6-8 所示。
![emds 0608]()
图 6-8. 每个输入值在逻辑曲线上都有相应的可能性
我们从逻辑回归曲线上方或下方获取每个点的可能性。如果点在逻辑回归曲线下方,我们需要从 1.0 中减去结果概率,因为我们也想最大化假阳性。
给定系数β[0] = –3.17576395 和β[1] = 0.69267212,示例 6-4 展示了我们如何在 Python 中计算这些数据的联合概率。
示例 6-4. 计算给定逻辑回归观察到所有点的联合概率
import math
import pandas as pd
patient_data = pd.read_csv('https://bit.ly/33ebs2R', delimiter=",").itertuples()
b0 = -3.17576395
b1 = 0.69267212
def logistic_function(x):
p = 1.0 / (1.0 + math.exp(-(b0 + b1 * x)))
return p
# Calculate the joint likelihood
joint_likelihood = 1.0
for p in patient_data:
if p.y == 1.0:
joint_likelihood *= logistic_function(p.x)
elif p.y == 0.0:
joint_likelihood *= (1.0 - logistic_function(p.x))
print(joint_likelihood) # 4.7911180221699105e-05
这里有一个数学技巧,我们可以用来压缩那个if
表达式。正如我们在第一章中讨论的,当你将任何数的幂设为 0 时,结果总是 1。看看这个公式,并注意指数中对真(1)和假(0)情况的处理:
要在 Python 中实现这一点,将for
循环中的所有内容压缩到示例 6-5 中。
示例 6-5. 在不使用if
表达式的情况下压缩联合概率计算
for p in patient_data:
joint_likelihood *= logistic_function(p.x) ** p.y * \
(1.0 - logistic_function(p.x)) ** (1.0 - p.y)
我到底做了什么?注意这个表达式有两个部分,一个是当时,另一个是当时。当任何数被提升到指数 0 时,结果将为 1。因此,无论y是 1 还是 0,它都会导致另一侧的条件评估为 1 并且在乘法中没有影响。我们可以用数学表达式完全表达我们的if
表达式。我们无法对使用if
的表达式进行导数,所以这将很有帮助。
请注意,计算机可能会因为乘以多个小小数而不堪重负,这被称为浮点下溢。这意味着随着小数变得越来越小,可能会在乘法中发生,计算机在跟踪那么多小数位数时会遇到限制。有一个巧妙的数学技巧可以解决这个问题。你可以对要相乘的每个小数取log()
,然后将它们相加。这要归功于我们在第一章中介绍的对数的加法性质。这样更加数值稳定,然后你可以调用exp()
函数将总和转换回来得到乘积。
让我们修改我们的代码,使用对数加法而不是乘法(参见示例 6-6)。请注意,log()
函数默认为基数e,虽然任何基数在技术上都可以工作,但这是首选,因为是其自身的导数,计算上更有效率。
示例 6-6. 使用对数加法
# Calculate the joint likelihood
joint_likelihood = 0.0
for p in patient_data:
joint_likelihood += math.log(logistic_function(p.x) ** p.y * \
(1.0 - logistic_function(p.x)) ** (1.0 - p.y))
joint_likelihood = math.exp(joint_likelihood)
要用数学符号表示前面的 Python 代码:
你想要计算前述表达式中的和的偏导数吗?我觉得不会。这太复杂了。天哪,在 SymPy 中表达那个函数本身就是一大口水!看看示例 6-7 中的内容。
示例 6-7. 在 SymPy 中表达逻辑回归的联合似然
joint_likelihood = Sum(log((1.0 / (1.0 + exp(-(b + m * x(i)))))**y(i) * \
(1.0 - (1.0 / (1.0 + exp(-(b + m * x(i))))))**(1-y(i))), (i, 0, n))
所以让我们让 SymPy 为我们做和的偏导数。然后我们将立即编译并使用它们进行梯度下降,如示例 6-8 所示。
示例 6-8. 在逻辑回归上使用梯度下降
from sympy import *
import pandas as pd
points = list(pd.read_csv("https://tinyurl.com/y2cocoo7").itertuples())
b1, b0, i, n = symbols('b1 b0 i n')
x, y = symbols('x y', cls=Function)
joint_likelihood = Sum(log((1.0 / (1.0 + exp(-(b0 + b1 * x(i))))) ** y(i) \
* (1.0 - (1.0 / (1.0 + exp(-(b0 + b1 * x(i)))))) ** (1 - y(i))), (i, 0, n))
# Partial derivative for m, with points substituted
d_b1 = diff(joint_likelihood, b1) \
.subs(n, len(points) - 1).doit() \
.replace(x, lambda i: points[i].x) \
.replace(y, lambda i: points[i].y)
# Partial derivative for m, with points substituted
d_b0 = diff(joint_likelihood, b0) \
.subs(n, len(points) - 1).doit() \
.replace(x, lambda i: points[i].x) \
.replace(y, lambda i: points[i].y)
# compile using lambdify for faster computation
d_b1 = lambdify([b1, b0], d_b1)
d_b0 = lambdify([b1, b0], d_b0)
# Perform Gradient Descent
b1 = 0.01
b0 = 0.01
L = .01
for j in range(10_000):
b1 += d_b1(b1, b0) * L
b0 += d_b0(b1, b0) * L
print(b1, b0)
# 0.6926693075370812 -3.175751550409821
在计算和的偏导数后,我们将 x 值和 y 值以及数据点数n代入。然后我们使用lambdify()
来编译导数函数以提高效率(它在幕后使用 NumPy)。之后,我们执行梯度下降,就像我们在第五章中所做的那样,但由于我们试图最大化而不是最小化,我们将每次调整添加到和中,而不是像最小二乘法中那样减去。
如您在示例 6-8 中所看到的,我们得到了β[0] = –3.17575 和β[1] = 0.692667。这与我们之前在 SciPy 中得到的系数值非常相似。
正如我们在第五章中学到的那样,我们也可以使用随机梯度下降,每次迭代只对一个或少数几个记录进行采样。这将延伸增加计算速度和性能的好处,同时防止过拟合。在这里再次覆盖将是多余的,所以我们将继续前进。
多变量逻辑回归
让我们尝试一个使用多个输入变量进行逻辑回归的示例。表 6-1 展示了一个虚构数据集中一些就业保留数据的样本(完整数据集在这里)。
表 6-1。就业保留数据样本
性别 |
年龄 |
晋升次数 |
工龄 |
是否离职 |
1 |
32 |
3 |
7 |
0 |
1 |
34 |
2 |
5 |
0 |
1 |
29 |
2 |
5 |
1 |
0 |
42 |
4 |
10 |
0 |
1 |
43 |
4 |
10 |
0 |
此数据集中有 54 条记录。假设我们想要用它来预测其他员工是否会离职,这里可以使用逻辑回归(尽管这不是一个好主意,稍后我会详细说明原因)。回想一下,我们可以支持多个输入变量,如下公式所示:
我将为每个变量sex
、age
、promotions
和years_employed
创建系数。输出变量did_quit
是二进制的,这将驱动我们正在预测的逻辑回归结果。因为我们处理多个维度,所以很难可视化我们的逻辑曲线所代表的曲线超平面。因此,我们将避免可视化。
让我们来点有趣的。我们将使用 scikit-learn,但创建一个交互式 shell,我们可以用来测试员工。示例 6-9 展示了代码,当我们运行它时,将执行逻辑回归,然后我们可以输入新员工以预测他们是否会离职。会出什么问题呢?我相信没有。我们只是根据人们的个人属性进行预测并做出相应决策。我相信一切都会好的。
(如果不清楚,我是在开玩笑)。
示例 6-9。对员工数据进行多变量逻辑回归
import pandas as pd
from sklearn.linear_model import LogisticRegression
employee_data = pd.read_csv("https://tinyurl.com/y6r7qjrp")
# grab independent variable columns
inputs = employee_data.iloc[:, :-1]
# grab dependent "did_quit" variable column
output = employee_data.iloc[:, -1]
# build logistic regression
fit = LogisticRegression(penalty='none').fit(inputs, output)
# Print coefficients:
print("COEFFICIENTS: {0}".format(fit.coef_.flatten()))
print("INTERCEPT: {0}".format(fit.intercept_.flatten()))
# Interact and test with new employee data
def predict_employee_will_stay(sex, age, promotions, years_employed):
prediction = fit.predict([[sex, age, promotions, years_employed]])
probabilities = fit.predict_proba([[sex, age, promotions, years_employed]])
if prediction == [[1]]:
return "WILL LEAVE: {0}".format(probabilities)
else:
return "WILL STAY: {0}".format(probabilities)
# Test a prediction
while True:
n = input("Predict employee will stay or leave {sex},
{age},{promotions},{years employed}: ")
(sex, age, promotions, years_employed) = n.split(",")
print(predict_employee_will_stay(int(sex), int(age), int(promotions),
int(years_employed)))
图 6-9 显示了员工是否被预测会离职的结果。员工的性别为“1”,年龄为 34 岁,晋升 1 次,公司工作了 5 年。果然,预测是“将离开”。
![emds 0609]()
图 6-9。预测 34 岁员工,1 次晋升和 5 年工作经验是否会离职
请注意,predict_proba()
函数将输出两个值,第一个是 0(假)的概率,第二个是 1(真)的概率。
您会注意到sex
、age
、promotions
和years_employed
的系数按照这个顺序显示。通过系数的权重,您可以看到sex
和age
在预测中起到很小的作用(它们的权重接近 0)。然而,promotions
和years_employed
具有显著的权重分别为-2.504
和0.97
。这个玩具数据集的一个秘密是,如果员工每两年没有晋升就会离职,我制造了这个模式,我的逻辑回归确实捕捉到了这个模式,您也可以尝试对其他员工进行测试。然而,如果您超出了它训练的数据范围,预测可能会开始失效(例如,如果放入一个 70 岁的员工,三年没有晋升,很难说这个模型会做出什么,因为它没有关于那个年龄的数据)。
当然,现实生活并不总是如此干净。一个在公司工作了八年,从未晋升过的员工很可能对自己的角色感到满意,不会很快离开。如果是这种情况,年龄等变量可能会发挥作用并被赋予权重。然后当然我们可能会担心其他未被捕捉到的相关变量。查看以下警告以了解更多信息。
谨慎对人进行分类!
一个快速而肯定地让自己陷入困境的方法是收集关于人们的数据,并随意地用它来做预测。不仅可能引发数据隐私问题,还可能出现法律和公关问题,如果模型被发现具有歧视性。像种族和性别这样的输入变量可能在机器学习训练中被赋予权重。之后,这些人口统计学数据可能会导致不良结果,比如不被录用或被拒绝贷款。更极端的应用包括被监视系统错误标记或被拒绝刑事假释。还要注意,看似无害的变量,比如通勤时间,已被发现与歧视性变量相关联。
在撰写本文时,已有多篇文章引用机器学习歧视作为一个问题:
随着数据隐私法律的不断发展,谨慎处理个人数据是明智的。考虑自动决策将如何传播以及如何造成伤害。有时最好的做法是让“问题”保持原样,继续手动处理。
最后,在这个员工留存的例子中,想想这些数据是从哪里来的。是的,我虚构了这个数据集,但在现实世界中,你总是要质疑数据是如何生成的。这个样本是从多长时间内得出的?我们要回溯多久来寻找已经离职的员工?什么构成了留下的员工?他们现在是当前员工吗?我们怎么知道他们不会马上离职,从而成为一个假阴性?数据科学家很容易陷入分析数据所说的内容,但不质疑数据来源和内置的假设。
获取这些问题的答案的最佳方法是了解预测用途。是用来决定何时提升人员以留住他们吗?这会产生一种循环偏见,促使具有一组属性的人员晋升吗?当这些晋升开始成为新的训练数据时,这种偏见会得到确认吗?
这些都是重要的问题,甚至可能是令人不快的问题,会导致不必要的范围渗入项目中。如果你的团队或领导不欢迎对项目进行这种审查,考虑让自己担任一个不同的角色,让好奇心成为一种优势。
理解对数几率
此时,是时候讨论逻辑回归及其数学构成了。这可能有点令人眩晕,所以在这里要花点时间。如果你感到不知所措,随时可以稍后回顾这一部分。
从 20 世纪开始,数学家一直对将线性函数的输出缩放到 0 和 1 之间感兴趣,因此对于预测概率是有用的。对数几率,也称为对数函数,适用于逻辑回归的这一目的。
记得之前我指出指数值 是一个线性函数吗?再看看我们的逻辑函数:
这个被提升到 e 的线性函数被称为 对数几率 函数,它取得了感兴趣事件的对数几率。你可能会说,“等等,我看不到 log()
或几率。我只看到一个线性函数!” 请耐心等待,我会展示隐藏的数学。
举个例子,让我们使用之前的逻辑回归,其中Β[0] = -3.17576395,Β[1] = 0.69267212。在六小时后出现症状的概率是多少,其中?我们已经知道如何做到这一点:将这些值代入我们的逻辑函数:
我们将这些值代入并输出概率为 0.72716。但让我们从赔率的角度来看这个问题。回想一下,在第二章中我们学习了如何从概率计算赔率:
因此,在六小时时,患者出现症状的可能性是不出现症状的 2.66517 倍。
当我们将赔率函数包装在一个自然对数(以e为底的对数)中时,我们称之为对数几率函数。这个公式的输出就是我们所说的对数几率,之所以这样命名...令人震惊...是因为我们取了赔率的对数:
我们在六小时时的对数几率为 0.9802687。这意味着什么,为什么我们要关心呢?当我们处于“对数几率领域”时,比较一组赔率相对容易。我们将大于 0 的任何值视为支持事件发生的赔率,而小于 0 的任何值则反对事件发生。对数几率为-1.05 与 0 的距离与 1.05 相同。然而,在普通赔率中,这些等价值分别为 0.3499 和 2.857,这并不容易解释。这就是对数几率的便利之处。
赔率和对数
对数和赔率有着有趣的关系。当赔率在 0.0 和 1.0 之间时,事件是不利的,但大于 1.0 的任何值都支持事件,并延伸到正无穷。这种缺乏对称性很尴尬。然而,对数重新调整了赔率,使其完全线性,其中对数几率为 0.0 表示公平的赔率。对数几率为-1.05 与 0 的距离与 1.05 相同,因此比较赔率更容易。
Josh Starmer 有一个很棒的视频讲述了赔率和对数之间的关系。
记得我说过我们逻辑回归公式中的线性函数是我们的对数几率函数。看看这个:
这与我们先前计算的值 0.98026877 相同,即在x = 6 时逻辑回归的赔率,然后取其log()
!那么这之间有什么联系?是什么将所有这些联系在一起?给定逻辑回归p的概率和输入变量x,就是这样:
让我们将对数几率线与逻辑回归一起绘制,如图 6-10 所示。
![emds 0610]()
图 6-10。对数几率线转换为输出概率的逻辑函数
每个逻辑回归实际上都由一个线性函数支持,而该线性函数是一个对数几率函数。请注意,在图 6-10 中,当对数几率在线上为 0.0 时,逻辑曲线的概率为 0.5。这是有道理的,因为当我们的赔率为 1.0 时,概率将为 0.50,如逻辑回归所示,而对数几率将为 0,如线所示。
从赔率的角度看逻辑回归的另一个好处是我们可以比较一个 x 值和另一个 x 值之间的效果。假设我想了解暴露于化学物质六小时和八小时之间我的赔率变化有多大。我可以取六小时和八小时的赔率,然后将两个赔率相对于彼此进行比较,得到一个赔率比。这不应与普通赔率混淆,是的,它是一个比率,但不是赔率比。
让我们首先找出分别为六小时和八小时的症状概率:
现在让我们将这些转换为几率,我们将其声明为:
最后,将两个几率相互对比作为一个几率比值,其中 8 小时的几率为分子,6 小时的几率为分母。我们得到一个约为 3.996 的数值,这意味着我们的症状出现的几率随着额外两小时的暴露增加了近四倍:
你会发现这个几率比值为 3.996 的数值在任何两小时范围内都成立,比如 2 小时到 4 小时,4 小时到 6 小时,8 小时到 10 小时等等。只要是两小时的间隔,你会发现几率比值保持一致。对于其他范围长度,它会有所不同。
R-平方
我们在第五章中已经涵盖了线性回归的许多统计指标,我们将尝试对 logistic 回归做同样的事情。我们仍然担心许多与线性回归相同的问题,包括过拟合和方差。事实上,我们可以借鉴并调整几个线性回归的指标,并将它们应用于 logistic 回归。让我们从开始。
就像线性回归一样,给定 logistic 回归也有一个。如果你回忆一下第五章,表示一个给定自变量解释因变量的程度。将这应用于我们的化学暴露问题,我们想要衡量化学暴露小时数解释症状出现的程度是有意义的。
实际上并没有关于如何计算 logistic 回归的 的最佳方法的共识,但一种被称为麦克法登伪 的流行技术紧密模仿了线性回归中使用的。我们将在以下示例中使用这种技术,以下是公式:
我们将学习如何计算“对数似然拟合”和“对数似然”,以便计算。
我们不能像线性回归那样在这里使用残差,但我们可以将结果投影回逻辑曲线,如图 6-11 所示,并查找它们在 0.0 和 1.0 之间的相应可能性。
![emds 0611]()
图 6-11. 将输出值投影回逻辑曲线
然后我们可以取这些可能性的log()
并将它们相加。这将是拟合的对数似然(示例 6-10)。就像我们计算最大似然一样,我们通过从 1.0 中减去“假”可能性来转换“假”可能性。
示例 6-10. 计算拟合的对数似然
from math import log, exp
import pandas as pd
patient_data = pd.read_csv('https://bit.ly/33ebs2R', delimiter=",").itertuples()
b0 = -3.17576395
b1 = 0.69267212
def logistic_function(x):
p = 1.0 / (1.0 + exp(-(b0 + b1 * x)))
return p
# Sum the log-likelihoods
log_likelihood_fit = 0.0
for p in patient_data:
if p.y == 1.0:
log_likelihood_fit += log(logistic_function(p.x))
elif p.y == 0.0:
log_likelihood_fit += log(1.0 - logistic_function(p.x))
print(log_likelihood_fit) # -9.946161673231583
使用一些巧妙的二进制乘法和 Python 推导,我们可以将那个for
循环和if
表达式整合成一行,返回log_likelihood_fit
。类似于我们在最大似然公式中所做的,我们可以使用一些真假病例之间的二进制减法来在数学上消除其中一个。在这种情况下,我们乘以 0,因此相应地将真或假情况应用于总和(示例 6-11)。
示例 6-11. 将我们的对数似然逻辑整合成一行
log_likelihood_fit = sum(log(logistic_function(p.x)) * p.y +
log(1.0 - logistic_function(p.x)) * (1.0 - p.y)
for p in patient_data)
如果我们要用数学符号来表达拟合的可能性,它会是这个样子。注意是给定输入变量的逻辑函数:
如在示例 6-10 和 6-11 中计算的,我们的拟合对数似然为-9.9461。我们需要另一个数据点来计算:即估计不使用任何输入变量,只使用真实病例数除以所有病例数(实际上只留下截距)。请注意,我们可以通过将所有 y 值相加来计算症状病例的数量,因为只有 1 会计入总和,而 0 不会。这是公式:
这是应用于示例 6-12 中的公式的 Python 等效展开。
示例 6-12. 患者的对数似然
import pandas as pd
from math import log, exp
patient_data = list(pd.read_csv('https://bit.ly/33ebs2R', delimiter=",") \
.itertuples())
likelihood = sum(p.y for p in patient_data) / len(patient_data)
log_likelihood = 0.0
for p in patient_data:
if p.y == 1.0:
log_likelihood += log(likelihood)
elif p.y == 0.0:
log_likelihood += log(1.0 - likelihood)
print(log_likelihood) # -14.341070198709906
为了整合这个逻辑并反映公式,我们可以将那个for
循环和if
表达式压缩成一行,使用一些二进制乘法逻辑来处理真假病例(示例 6-13)。
示例 6-13. 将对数似然整合成一行
log_likelihood = sum(log(likelihood)*p.y + log(1.0 - likelihood)*(1.0 - p.y) \
for p in patient_data)
最后,只需将这些值代入并获得你的:
下面是 Python 代码,显示在示例 6-14 中,完整计算。
示例 6-14. 计算逻辑回归的
import pandas as pd
from math import log, exp
patient_data = list(pd.read_csv('https://bit.ly/33ebs2R', delimiter=",") \
.itertuples())
# Declare fitted logistic regression
b0 = -3.17576395
b1 = 0.69267212
def logistic_function(x):
p = 1.0 / (1.0 + exp(-(b0 + b1 * x)))
return p
# calculate the log likelihood of the fit
log_likelihood_fit = sum(log(logistic_function(p.x)) * p.y +
log(1.0 - logistic_function(p.x)) * (1.0 - p.y)
for p in patient_data)
# calculate the log likelihood without fit
likelihood = sum(p.y for p in patient_data) / len(patient_data)
log_likelihood = sum(log(likelihood) * p.y + log(1.0 - likelihood) * (1.0 - p.y) \
for p in patient_data)
# calculate R-Square
r2 = (log_likelihood - log_likelihood_fit) / log_likelihood
print(r2) # 0.306456105756576
好的,所以我们得到一个 = 0.306456,那么化学暴露时间是否能解释某人是否出现症状?正如我们在第五章中学到的线性回归一样,拟合不好的情况下会接近 0.0,而拟合较好的情况下会接近 1.0。因此,我们可以得出结论,暴露时间对于预测症状是一般般的,因为为 0.30645。除了时间暴露之外,肯定还有其他变量更好地预测某人是否会出现症状。这是有道理的,因为我们观察到的大多数数据中,有很多患者出现症状和没有出现症状的混合,如图 6-12 所示。
![emds 0612]()
图 6-12. 我们的数据中间有很多方差,因此我们的数据有一个中等的为 0.30645
但是如果我们的数据中有一个明确的分界线,其中 1 和 0 的结果被清晰地分开,如图 6-13 所示,我们将有一个完美的为 1.0。
![emds 0613]()
图 6-13. 这个逻辑回归有一个完美的为 1.0,因为根据暴露时间预测的结果有一个清晰的分界线
P-值
就像线性回归一样,我们并不因为有一个就结束了。我们需要调查我们看到这些数据是因为偶然还是因为实际关系的可能性。这意味着我们需要一个 p 值。
为了做到这一点,我们需要学习一个称为卡方分布的新概率分布,标记为分布。它是连续的,在统计学的几个领域中使用,包括这个!
如果我们取标准正态分布中的每个值(均值为 0,标准差为 1)并平方,那将给我们一个自由度为 1 的分布。对于我们的目的,自由度将取决于我们逻辑回归中有多少参数,这将是。你可以在图 6-14 中看到不同自由度的例子。
![emds 0614]()
图 6-14。不同自由度的分布
由于我们有两个参数(暴露时间和是否出现症状),我们的自由度将为 1,因为。
我们需要在前一小节关于 R²中计算的对数似然拟合和对数似然。这是产生我们需要查找的值的公式:
然后我们取该值并从分布中查找概率。这将给我们我们的 p 值:
示例 6-15 展示了我们拟合的逻辑回归的 p 值。我们使用 SciPy 的chi2
模块使用卡方分布。
示例 6-15。计算给定逻辑回归的 p 值
import pandas as pd
from math import log, exp
from scipy.stats import chi2
patient_data = list(pd.read_csv('https://bit.ly/33ebs2R', delimiter=",").itertuples())
# Declare fitted logistic regression
b0 = -3.17576395
b1 = 0.69267212
def logistic_function(x):
p = 1.0 / (1.0 + exp(-(b0 + b1 * x)))
return p
# calculate the log likelihood of the fit
log_likelihood_fit = sum(log(logistic_function(p.x)) * p.y +
log(1.0 - logistic_function(p.x)) * (1.0 - p.y)
for p in patient_data)
# calculate the log likelihood without fit
likelihood = sum(p.y for p in patient_data) / len(patient_data)
log_likelihood = sum(log(likelihood) * p.y + log(1.0 - likelihood) * (1.0 - p.y) \
for p in patient_data)
# calculate p-value
chi2_input = 2 * (log_likelihood_fit - log_likelihood)
p_value = chi2.pdf(chi2_input, 1) # 1 degree of freedom (n - 1)
print(p_value) # 0.0016604875618753787
所以我们有一个 p 值为 0.00166,如果我们的显著性阈值为 0.05,我们说这些数据在统计上是显著的,不是随机事件。
训练/测试拆分
如第五章中所述,我们可以使用训练/测试拆分来验证机器学习算法。这是评估逻辑回归性能的更机器学习方法。虽然依赖传统统计指标如和 p 值是个好主意,但当你处理更多变量时,这变得不太实际。这时再次用到了训练/测试拆分。回顾一下,图 6-15 展示了一个三折交叉验证交替测试数据集。
![emds 0615]()
图 6-15。一个三折交叉验证,交替将数据集的每个第三部分作为测试数据集
在示例 6-16 中,我们对员工留存数据集执行逻辑回归,但将数据分成三部分。然后我们交替将每个部分作为测试数据。最后,我们用平均值和标准差总结三个准确性。
示例 6-16。执行带有三折交叉验证的逻辑回归
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import KFold, cross_val_score
# Load the data
df = pd.read_csv("https://tinyurl.com/y6r7qjrp", delimiter=",")
X = df.values[:, :-1]
Y = df.values[:, -1]
# "random_state" is the random seed, which we fix to 7
kfold = KFold(n_splits=3, random_state=7, shuffle=True)
model = LogisticRegression(penalty='none')
results = cross_val_score(model, X, Y, cv=kfold)
print("Accuracy Mean: %.3f (stdev=%.3f)" % (results.mean(), results.std()))
我们还可以使用随机折叠验证、留一交叉验证以及我们在第五章中执行的所有其他折叠变体。说完这些,让我们谈谈为什么准确率是分类的一个糟糕的度量。
混淆矩阵
假设一个模型观察到名为“迈克尔”的人离职。捕获名字作为输入变量的原因确实值得怀疑,因为一个人的名字是否会影响他们是否离职是可疑的。然而,为了简化例子,让我们继续。该模型随后预测任何名为“迈克尔”的人都会离职。
现在这就是准确率失效的地方。我有一百名员工,其中一个名叫“迈克尔”,另一个名叫“山姆”。迈克尔被错误地预测会离职,而实际上是山姆离职了。我的模型准确率是多少?是 98%,因为在一百名员工中只有两次错误预测,如图 6-16 所示。
![emds 0616]()
图 6-16. 预测名为“迈克尔”的员工会离职,但实际上是另一名员工离职,给我们带来了 98%的准确率
尤其对于数据不平衡的情况,其中感兴趣的事件(例如,员工离职)很少见,准确率指标对于分类问题是极其误导的。如果供应商、顾问或数据科学家试图通过准确性来销售分类系统,请要求一个混淆矩阵。
混淆矩阵是一个网格,将预测与实际结果进行对比,显示出真正的正例、真正的负例、假正例(I 型错误)和假负例(II 型错误)。这里是一个在图 6-17 中呈现的混淆矩阵。
![emds 0617]()
图 6-17. 一个简单的混淆矩阵
通常,我们希望对角线值(从左上到右下)更高,因为这些反映了正确的分类。我们想评估有多少被预测会离职的员工实际上确实离职(真正的正例)。相反,我们也想评估有多少被预测会留下的员工实际上确实留下(真正的负例)。
其他单元格反映了错误的预测,其中一个被预测会离职的员工最终留下(假正例),以及一个被预测会留下的员工最终离职(假负例)。
我们需要将准确率指标分解为针对混淆矩阵不同部分的更具体的准确率指标。让我们看看图 6-18,它添加了一些有用的指标。
从混淆矩阵中,我们可以得出除准确率之外的各种有用的指标。我们可以清楚地看到精确度(正面预测的准确性)和灵敏度(识别出的正面率)都为 0,这意味着这个机器学习模型在正面预测上完全失败。
![emds 0618]()
图 6-18. 在混淆矩阵中添加有用的指标
示例 6-17 展示了如何在 SciPy 中使用混淆矩阵 API 对具有训练/测试拆分的逻辑回归进行操作。请注意,混淆矩阵仅应用于测试数据集。
示例 6-17。在 SciPy 中为测试数据集创建混淆矩阵
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix
from sklearn.model_selection import train_test_split
# Load the data
df = pd.read_csv('https://bit.ly/3cManTi', delimiter=",")
# Extract input variables (all rows, all columns but last column)
X = df.values[:, :-1]
# Extract output column (all rows, last column)\
Y = df.values[:, -1]
model = LogisticRegression(solver='liblinear')
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=.33,
random_state=10)
model.fit(X_train, Y_train)
prediction = model.predict(X_test)
"""
The confusion matrix evaluates accuracy within each category.
[[truepositives falsenegatives]
[falsepositives truenegatives]]
The diagonal represents correct predictions,
so we want those to be higher
"""
matrix = confusion_matrix(y_true=Y_test, y_pred=prediction)
print(matrix)
贝叶斯定理和分类
你还记得第二章中的贝叶斯定理吗?你可以使用贝叶斯定理引入外部信息来进一步验证混淆矩阵的发现。图 6-19 展示了对一千名患者进行疾病测试的混淆矩阵。
![emds 0619]()
图 6-19。用于识别疾病的医学测试的混淆矩阵
我们被告知对于存在健康风险的患者,将成功识别 99%(敏感性)。使用混淆矩阵,我们可以数学上验证这一点:
但是如果我们颠倒条件呢?那么测试结果为阳性的人中有多少百分比存在健康风险(精确度)?虽然我们在颠倒条件概率,但在这里我们不必使用贝叶斯定理,因为混淆矩阵为我们提供了所有我们需要的数字:
好吧,79.8%并不算糟糕,这是测试为阳性的人中实际患有疾病的百分比。但请问自己...我们对数据做了什么假设?它是否代表人口?
一些快速研究发现 1%的人口实际上患有这种疾病。在这里有机会使用贝叶斯定理。我们可以考虑实际患有疾病的人口比例,并将其纳入我们的混淆矩阵结果中。然后我们发现了一些重要的东西。
当我们考虑到只有 1%的人口存在风险,并且我们测试的患者中有 20%存在风险时,接受阳性测试的人存在风险的概率为 3.39%!为什么从 99%下降到了这个数字?这只是展示了我们如何容易被高概率欺骗,这种概率只在特定样本中高,比如供应商的一千名测试患者。因此,如果这个测试只有 3.39%的概率成功识别真正的阳性,我们可能不应该使用它。
接收器操作特性/曲线下面积
当我们评估不同的机器学习配置时,可能会得到数十、数百或数千个混淆矩阵。这些可能很繁琐,因此我们可以用一个接收器操作特性(ROC)曲线来总结所有这些,如图 6-20 所示。这使我们能够看到每个测试实例(每个由一个黑点表示)并找到真正例和假正例之间的一个令人满意的平衡。
我们还可以通过为每个模型创建单独的 ROC 曲线来比较不同的机器学习模型。例如,如果在图 6-21 中,我们的顶部曲线代表逻辑回归,底部曲线代表决策树(这是本书中没有涵盖的一种机器学习技术),我们可以并排看到它们的性能。曲线下面积(AUC)是选择使用哪个模型的良好指标。由于顶部曲线(逻辑回归)的面积更大,这表明它是一个更优秀的模型。
![emds 0620]()
图 6-20. 接收器操作特性曲线
![emds 0621]()
图 6-21. 通过它们各自的 ROC 曲线比较两个模型的曲线下面积(AUC)
要将 AUC 作为评分指标使用,请在 scikit-learn API 中将scoring
参数更改为使用roc_auc
,如在示例 6-18 中所示进行交叉验证。
示例 6-18. 使用 AUC 作为 scikit-learn 参数
# put Scikit_learn model here
results = cross_val_score(model, X, Y, cv=kfold, scoring='roc_auc')
print("AUC: %.3f (%.3f)" % (results.mean(), results.std()))
# AUC: 0.791 (0.051)
类别不平衡
在我们结束本章之前,还有一件事情需要讨论。正如我们之前在讨论混淆矩阵时看到的那样,类别不平衡,即当数据在每个结果类别中没有得到平等表示时,是机器学习中的一个问题。不幸的是,许多感兴趣的问题都存在不平衡,比如疾病预测、安全漏洞、欺诈检测等等。类别不平衡仍然是一个没有很好解决方案的问题。然而,你可以尝试一些技术。
首先,你可以做一些明显的事情,比如收集更多数据或尝试不同的模型,以及使用混淆矩阵和 ROC/AUC 曲线。所有这些都将有助于跟踪糟糕的预测并主动捕捉错误。
另一种常见的技术是复制少数类别中的样本,直到它在数据集中得到平等的表示。你可以在 scikit-learn 中这样做,如示例 6-19 所示,在进行训练-测试拆分时传递包含类别值的列的stratify
选项,它将尝试平等地表示每个类别的数据。
示例 6-19. 在 scikit-learn 中使用stratify
选项来平衡数据中的类别
X, Y = ...
X_train, X_test, Y_train, Y_test = \
train_test_split(X, Y, test_size=.33, stratify=Y)
还有一类算法称为 SMOTE,它会生成少数类的合成样本。然而,最理想的方法是以一种利用异常检测模型的方式解决问题,这些模型专门设计用于寻找罕见事件。它们寻找异常值,但不一定是分类,因为它们是无监督算法。所有这些技术超出了本书的范围,但值得一提,因为它们可能为给定问题提供更好的解决方案。
结论
逻辑回归是预测数据概率和分类的主力模型。逻辑回归可以预测不止一个类别,而不仅仅是真/假。你只需构建单独的逻辑回归模型,模拟它是否属于该类别,产生最高概率的模型即为胜者。你可能会发现,大部分情况下,scikit-learn 会为你完成这一点,并在数据具有两个以上类别时进行检测。
在本章中,我们不仅介绍了如何使用梯度下降和 scikit-learn 拟合逻辑回归,还涉及了统计学和机器学习方法的验证。在统计学方面,我们涵盖了和 p 值,而在机器学习方面,我们探讨了训练/测试拆分、混淆矩阵和 ROC/AUC。
如果你想了解更多关于逻辑回归的知识,可能最好的资源是 Josh Starmer 的 StatQuest 播放列表关于逻辑回归的视频。我必须感谢 Josh 在协助本章某些部分方面的工作,特别是如何计算逻辑回归的和 p 值。如果没有其他,观看他的视频也是为了那些出色的开场曲!
一如既往,你将发现自己在统计学和机器学习两个世界之间徘徊。目前许多书籍和资源都从机器学习的角度讨论逻辑回归,但也尝试寻找统计学资源。两种思维方式各有优劣,只有适应两者才能取得胜利!
练习
提供了一个包含三个输入变量RED
、GREEN
和BLUE
以及一个输出变量LIGHT_OR_DARK_FONT_IND
的数据集,链接在这里。这将用于预测给定背景颜色(由 RGB 值指定)是否适合使用浅色/深色字体(分别为 0/1)。
-
对上述数据执行逻辑回归,使用三折交叉验证和准确率作为评估指标。
-
生成一个混淆矩阵,比较预测和实际数据。
-
选择几种不同的背景颜色(你可以使用像这个这样的 RGB 工具),看看逻辑回归是否明智地为每种颜色选择浅色(0)或深色(1)字体。
-
根据前面的练习,你认为逻辑回归对于预测给定背景颜色的浅色或深色字体有效吗?
答案在附录 B 中。
第七章:神经网络
在过去的 10 年里,一种经历复兴的回归和分类技术是神经网络。在最简单的定义中,神经网络是一个包含权重、偏差和非线性函数层的多层回归,位于输入变量和输出变量之间。深度学习是神经网络的一种流行变体,利用包含权重和偏差的多个“隐藏”(或中间)层节点。每个节点在传递给非线性函数(称为激活函数)之前类似于一个线性函数。就像我们在第五章中学到的线性回归一样,优化技术如随机梯度下降被用来找到最优的权重和偏差值以最小化残差。
神经网络为计算机以前难以解决的问题提供了令人兴奋的解决方案。从识别图像中的物体到处理音频中的单词,神经网络已经创造了影响我们日常生活的工具。这包括虚拟助手和搜索引擎,以及我们 iPhone 中的照片工具。
鉴于媒体的炒作和大胆宣称主导着关于神经网络的新闻头条,也许令人惊讶的是,它们自上世纪 50 年代以来就存在了。它们在 2010 年后突然变得流行的原因是由于数据和计算能力的不断增长。2011 年至 2015 年之间的 ImageNet 挑战赛可能是复兴的最大推动力,将对 140 万张图像进行一千个类别的分类准确率提高到了 96.4%。
然而,就像任何机器学习技术一样,它只适用于狭义定义的问题。即使是创建“自动驾驶”汽车的项目也不使用端到端的深度学习,主要使用手工编码的规则系统,卷积神经网络充当“标签制造机”来识别道路上的物体。我们将在本章后面讨论这一点,以了解神经网络实际上是如何使用的。但首先我们将在 NumPy 中构建一个简单的神经网络,然后使用 scikit-learn 作为库实现。
何时使用神经网络和深度学习
神经网络和深度学习可用于分类和回归,那么它们与线性回归、逻辑回归和其他类型的机器学习相比如何?你可能听说过“当你手中只有一把锤子时,所有事情看起来都像钉子”。每种算法都有其特定情况下的优势和劣势。线性回归、逻辑回归以及梯度提升树(本书未涵盖)在结构化数据上做出了相当出色的预测。将结构化数据视为可以轻松表示为表格的数据,具有行和列。但感知问题如图像分类则不太结构化,因为我们试图找到像素组之间的模糊相关性以识别形状和模式,而不是表格中的数据行。尝试预测正在输入的句子中的下四五个单词,或者解密音频剪辑中的单词,也是感知问题,是神经网络用于自然语言处理的例子。
在本章中,我们将主要关注只有一个隐藏层的简单神经网络。
使用神经网络是否有些大材小用?
对于即将介绍的例子来说,使用神经网络可能有些大材小用,因为逻辑回归可能更实用。甚至可以使用公式方法。然而,我一直是一个喜欢通过将复杂技术应用于简单问题来理解的人。你可以了解技术的优势和局限性,而不会被大型数据集所分散注意力。因此,请尽量不要在更实用的情况下使用神经网络。为了理解技术,我们将在本章中打破这个规则。
一个简单的神经网络
这里有一个简单的例子,让你对神经网络有所了解。我想要预测给定颜色背景下的字体应该是浅色(1)还是深色(0)。以下是不同背景颜色的几个示例,见图 7-1。顶部一行最适合浅色字体,底部一行最适合深色字体。
![emds 0701]()
图 7-1. 浅色背景颜色最适合深色字体,而深色背景颜色最适合浅色字体
在计算机科学中,表示颜色的一种方式是使用 RGB 值,即红色、绿色和蓝色值。每个值都介于 0 和 255 之间,表示这三种颜色如何混合以创建所需的颜色。例如,如果我们将 RGB 表示为(红色,绿色,蓝色),那么深橙色的 RGB 值为(255,140,0),粉色为(255,192,203)。黑色为(0,0,0),白色为(255,255,255)。
从机器学习和回归的角度来看,我们有三个数值输入变量red
、green
和blue
来捕捉给定背景颜色。我们需要对这些输入变量拟合一个函数,并输出是否应该为该背景颜色使用浅色(1)或深色(0)字体。
通过 RGB 表示颜色
在线有数百种颜色选择器调色板可供尝试 RGB 值。W3 Schools 有一个这里。
请注意,这个例子与神经网络识别图像的工作原理并不相去甚远,因为每个像素通常被建模为三个数值 RGB 值。在这种情况下,我们只关注一个“像素”作为背景颜色。
让我们从高层次开始,把所有的实现细节放在一边。我们将以洋葱的方式来处理这个主题,从更高的理解开始,然后慢慢剥离细节。目前,这就是为什么我们简单地将一个接受输入并产生输出的过程标记为“神秘数学”。我们有三个数值输入变量 R、G 和 B,这些变量被这个神秘的数学处理。然后它输出一个介于 0 和 1 之间的预测,如图 7-2 所示。
![emds 0702]()
图 7-2。我们有三个数值 RGB 值用于预测使用浅色或深色字体
这个预测输出表示一个概率。输出概率是使用神经网络进行分类的最常见模型。一旦我们用它们的数值替换 RGB,我们会发现小于 0.5 会建议使用深色字体,而大于 0.5 会建议使用浅色字体,如图 7-3 所示。
![emds 0703]()
图 7-3。如果我们输入一个粉色的背景色(255,192,203),那么神秘的数学会推荐使用浅色字体,因为输出概率 0.89 大于 0.5
那个神秘的数学黑匣子里到底发生了什么?让我们在图 7-4 中看一看。
我们还缺少神经网络的另一个部分,即激活函数,但我们很快会讨论到。让我们先了解这里发生了什么。左侧的第一层只是三个变量的输入,这些变量在这种情况下是红色、绿色和蓝色值。在隐藏(中间)层中,请注意我们在输入和输出之间产生了三个节点,或者说是权重和偏置的函数。每个节点本质上是一个线性函数,斜率为,截距为,与输入变量相乘并求和。每个输入节点和隐藏节点之间有一个权重,每个隐藏节点和输出节点之间有另一组权重。每个隐藏和输出节点都会额外添加一个偏置。
![emds 0704]()
图 7-4。神经网络的隐藏层对每个输入变量应用权重和偏置值,输出层对该输出应用另一组权重和偏置
注意,输出节点重复执行相同的操作,将隐藏层的加权和求和输出作为输入传递到最终层,其中另一组权重和偏置将被应用。
简而言之,这是一个回归问题,就像线性回归或逻辑回归一样,但需要解决更多参数。权重和偏置值类似于m和b,或者线性回归中的和参数。我们使用随机梯度下降和最小化损失,就像线性回归一样,但我们需要一种称为反向传播的额外工具来解开权重和偏置值,并使用链式法则计算它们的偏导数。我们将在本章后面详细讨论这一点,但现在让我们假设我们已经优化了权重和偏置值。我们需要先讨论激活函数。
激活函数
接下来让我们介绍激活函数。激活函数是一个非线性函数,它转换或压缩节点中的加权和值,帮助神经网络有效地分离数据,以便进行分类。让我们看一下图 7-5。如果没有激活函数,你的隐藏层将毫无生产力,表现不会比线性回归好。
![emds 0705]()
图 7-5. 应用激活函数
ReLU 激活函数将隐藏节点的任何负输出归零。如果权重、偏置和输入相乘并求和得到负数,它将被转换为 0。否则输出保持不变。这是使用 SymPy (示例 7-1) 绘制的 ReLU 图(图 7-6)。
示例 7-1. 绘制 ReLU 函数
from sympy import *
# plot relu
x = symbols('x')
relu = Max(0, x)
plot(relu)
![emds 0706]()
图 7-6. ReLU 函数图
ReLU 是“修正线性单元”的缩写,但这只是一种将负值转换为 0 的花哨方式。ReLU 在神经网络和深度学习中的隐藏层中变得流行,因为它速度快,并且缓解了梯度消失问题。梯度消失发生在偏导数斜率变得非常小,导致过早接近 0 并使训练停滞。
输出层有一个重要的任务:它接收来自神经网络隐藏层的大量数学,并将其转换为可解释的结果,例如呈现分类预测。对于这个特定的神经网络,输出层使用逻辑激活函数,这是一个简单的 S 形曲线。如果您阅读第六章,逻辑(或 S 形)函数应该很熟悉,它表明逻辑回归在我们的神经网络中充当一层。输出节点的权重、偏置和从隐藏层传入的每个值求和。之后,它通过逻辑函数传递结果值,以便输出介于 0 和 1 之间的数字。就像第六章中的逻辑回归一样,这代表了输入到神经网络的给定颜色建议使用浅色字体的概率。如果大于或等于 0.5,则神经网络建议使用浅色字体,否则建议使用深色字体。
这是使用 SymPy (示例 7-2) 绘制的逻辑函数图(图 7-7)。
示例 7-2. SymPy 中的逻辑激活函数
from sympy import *
# plot logistic
x = symbols('x')
logistic = 1 / (1 + exp(-x))
plot(logistic)
![emds 0707]()
图 7-7. 逻辑激活函数
当我们通过激活函数传递节点的加权、偏置和求和值时,我们现在称之为激活输出,意味着它已通过激活函数进行了过滤。当激活输出离开隐藏层时,信号准备好被馈送到下一层。激活函数可能会增强、减弱或保持信号不变。这就是神经网络中大脑和突触的比喻的来源。
鉴于复杂性的潜在可能性,您可能想知道是否还有其他激活函数。一些常见的激活函数显示在表 7-1 中。
表 7-1. 常见激活函数
名称 |
典型使用层 |
描述 |
注释 |
线性 |
输出 |
保持值不变 |
不常用 |
Logistic |
输出层 |
S 形 sigmoid 曲线 |
将值压缩在 0 和 1 之间,通常用于二元分类 |
双曲正切 |
隐藏层 |
tanh,在 -1 和 1 之间的 S 形 sigmoid 曲线 |
通过将均值接近 0 来“居中”数据 |
ReLU |
隐藏层 |
将负值转换为 0 |
比 sigmoid 和 tanh 更快的流行激活函数,缓解消失梯度问题,计算成本低廉 |
Leaky ReLU |
隐藏层 |
将负值乘以 0.01 |
ReLU 的有争议变体,边缘化而不是消除负值 |
Softmax |
输出层 |
确保所有输出节点加起来为 1.0 |
适用于多类别分类,重新缩放输出使其加起来为 1.0 |
这不是激活函数的全面列表,理论上神经网络中任何函数都可以是激活函数。
尽管这个神经网络表面上支持两类(浅色或深色字体),但实际上它被建模为一类:字体是否应该是浅色(1)或不是(0)。如果您想支持多个类别,可以为每个类别添加更多输出节点。例如,如果您试图识别手写数字 0-9,将有 10 个输出节点,代表给定图像是这些数字的概率。当有多个类别时,您可能还考虑在输出时使用 softmax 激活。图 7-8 展示了一个以像素化图像为输入的示例,其中像素被分解为单独的神经网络输入,然后通过两个中间层,最后一个输出层,有 10 个节点代表 10 个类别的概率(数字 0-9)。
![emds 0708]()
图 7-8. 一个神经网络,将每个像素作为输入,并预测图像包含的数字
在神经网络上使用 MNIST 数据集的示例可以在 附录 A 中找到。
我不知道要使用什么激活函数!
如果不确定要使用哪些激活函数,当前最佳实践倾向于在中间层使用 ReLU,在输出层使用 logistic(sigmoid)。如果输出中有多个分类,可以在输出层使用 softmax。
前向传播
让我们使用 NumPy 捕获到目前为止学到的知识。请注意,我尚未优化参数(我们的权重和偏置值)。我们将用随机值初始化它们。
示例 7-3 是创建一个简单的前馈神经网络的 Python 代码,尚未进行优化。前馈 意味着我们只是将一种颜色输入到神经网络中,看看它输出什么。权重和偏置是随机初始化的,并将在本章后面进行优化,因此暂时不要期望有用的输出。
示例 7-3. 一个具有随机权重和偏置值的简单前向传播网络
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
all_data = pd.read_csv("https://tinyurl.com/y2qmhfsr")
# Extract the input columns, scale down by 255
all_inputs = (all_data.iloc[:, 0:3].values / 255.0)
all_outputs = all_data.iloc[:, -1].values
# Split train and test data sets
X_train, X_test, Y_train, Y_test = train_test_split(all_inputs, all_outputs,
test_size=1/3)
n = X_train.shape[0] # number of training records
# Build neural network with weights and biases
# with random initialization
w_hidden = np.random.rand(3, 3)
w_output = np.random.rand(1, 3)
b_hidden = np.random.rand(3, 1)
b_output = np.random.rand(1, 1)
# Activation functions
relu = lambda x: np.maximum(x, 0)
logistic = lambda x: 1 / (1 + np.exp(-x))
# Runs inputs through the neural network to get predicted outputs
def forward_prop(X):
Z1 = w_hidden @ X + b_hidden
A1 = relu(Z1)
Z2 = w_output @ A1 + b_output
A2 = logistic(Z2)
return Z1, A1, Z2, A2
# Calculate accuracy
test_predictions = forward_prop(X_test.transpose())[3] # grab only output layer, A2
test_comparisons = np.equal((test_predictions >= .5).flatten().astype(int), Y_test)
accuracy = sum(test_comparisons.astype(int) / X_test.shape[0])
print("ACCURACY: ", accuracy)
这里有几点需要注意。 包含 RGB 输入值以及输出值(1 代表亮,0 代表暗)的数据集包含在此 CSV 文件中。 我将输入列 R、G 和 B 的值缩小了 1/255 的因子,使它们介于 0 和 1 之间。 这将有助于后续的训练,以便压缩数字空间。
请注意,我还使用 scikit-learn 将数据的 2/3 用于训练,1/3 用于测试,我们在第五章中学习了如何做到这一点。 n
简单地是训练数据记录的数量。
现在请注意示例 7-4 中显示的代码行。
示例 7-4。 NumPy 中的权重矩阵和偏置向量
# Build neural network with weights and biases
# with random initialization
w_hidden = np.random.rand(3, 3)
w_output = np.random.rand(1, 3)
b_hidden = np.random.rand(3, 1)
b_output = np.random.rand(1, 1)
这些声明了我们神经网络隐藏层和输出层的权重和偏置。 这可能还不明显,但矩阵乘法将使我们的代码变得强大简单,使用线性代数和 NumPy。
权重和偏置将被初始化为介于 0 和 1 之间的随机值。 让我们首先看一下权重矩阵。 当我运行代码时,我得到了这些矩阵:
请注意, 是隐藏层中的权重。第一行代表第一个节点的权重 , 和 。第二行是第二个节点,带有权重 , 和 。第三行是第三个节点,带有权重 , 和 。
输出层只有一个节点,意味着其矩阵只有一行,带有权重 , 和 。
看到规律了吗?每个节点在矩阵中表示为一行。如果有三个节点,就有三行。如果只有一个节点,就只有一行。每一列保存着该节点的权重值。
让我们也看看偏差。由于每个节点有一个偏差,隐藏层将有三行偏差,输出层将有一行偏差。每个节点只有一个偏差,所以只有一列:
现在让我们将这些矩阵值与我们在 图 7-9 中展示的神经网络进行比较。
![emds 0709]()
图 7-9。将我们的神经网络与权重和偏差矩阵值进行可视化
那么除了紧凑难懂之外,这种矩阵形式中的权重和偏差有什么好处呢?让我们把注意力转向 示例 7-5 中的这些代码行。
示例 7-5。我们神经网络的激活函数和前向传播函数
# Activation functions
relu = lambda x: np.maximum(x, 0)
logistic = lambda x: 1 / (1 + np.exp(-x))
# Runs inputs through the neural network to get predicted outputs
def forward_prop(X):
Z1 = w_hidden @ X + b_hidden
A1 = relu(Z1)
Z2 = w_output @ A1 + b_output
A2 = logistic(Z2)
return Z1, A1, Z2, A2
这段代码很重要,因为它使用矩阵乘法和矩阵-向量乘法简洁地执行我们整个神经网络。我们在第四章中学习了这些操作。它仅用几行代码将三个 RGB 输入的颜色通过权重、偏置和激活函数运行。
首先,我声明relu()
和logistic()
激活函数,它们分别接受给定的输入值并返回曲线的输出值。forward_prop()
函数为包含 R、G 和 B 值的给定颜色输入X
执行整个神经网络。它返回四个阶段的矩阵输出:Z1
、A1
、Z2
和A2
。数字“1”和“2”表示操作属于第 1 层和第 2 层。“Z”表示来自该层的未激活输出,“A”是来自该层的激活输出。
隐藏层由Z1
和A1
表示。Z1
是应用于X
的权重和偏置。然后A1
获取来自Z1
的输出,并通过激活 ReLU 函数。Z2
获取来自A1
的输出,并应用输出层的权重和偏置。该输出依次通过激活函数,即逻辑曲线,变为A2
。最终阶段,A2
,是输出层的预测概率,返回一个介于 0 和 1 之间的值。我们称之为A2
,因为它是来自第 2 层的“激活”输出。
让我们更详细地分解这个,从Z1
开始:
首先,我们在和输入颜色X
之间执行矩阵-向量乘法。我们将的每一行(每一行都是一个节点的权重集)与向量X
(RGB 颜色输入值)相乘。然后将偏置添加到该结果中,如图 7-10 所示。
![emds 0710]()
图 7-10。将隐藏层权重和偏置应用于输入X
,使用矩阵-向量乘法以及向量加法
那个向量是隐藏层的原始输出,但我们仍然需要通过激活函数将转换为。很简单。只需将该向量中的每个值通过 ReLU 函数传递,就会给我们。因为所有值都是正值,所以不应该有影响。
现在让我们将隐藏层输出传递到最终层,得到,然后。 成为输出层的输入。
最后,将这个单个值通过激活函数传递,得到。这将产生一个约为 0.95425 的预测:
执行我们整个神经网络,尽管我们尚未对其进行训练。但请花点时间欣赏,我们已经将所有这些输入值、权重、偏差和非线性函数转化为一个将提供预测的单个值。
再次强调,A2
是最终输出,用于预测背景颜色是否需要浅色(1)或深色(1)字体。尽管我们的权重和偏差尚未优化,让我们按照示例 7-6 计算准确率。取测试数据集X_test
,转置它,并通过forward_prop()
函数传递,但只获取A2
向量,其中包含每个测试颜色的预测。然后将预测与实际值进行比较,并计算正确预测的百分比。
示例 7-6. 计算准确率
# Calculate accuracy
test_predictions = forward_prop(X_test.transpose())[3] # grab only A2
test_comparisons = np.equal((test_predictions >= .5).flatten().astype(int), Y_test)
accuracy = sum(test_comparisons.astype(int) / X_test.shape[0])
print("ACCURACY: ", accuracy)
当我运行示例 7-3 中的整个代码时,我大致获得 55%到 67%的准确率。请记住,权重和偏差是随机生成的,因此答案会有所不同。虽然这个准确率似乎很高,考虑到参数是随机生成的,但请记住,输出预测是二元的:浅色或深色。因此,对于每个预测,随机抛硬币也可能产生这种结果,所以这个数字不应令人惊讶。
不要忘记检查是否存在数据不平衡!
正如在第六章中讨论的那样,不要忘记分析数据以检查是否存在不平衡的类别。整个背景颜色数据集有点不平衡:512 种颜色的输出为 0,833 种颜色的输出为 1。这可能会使准确率产生偏差,也可能是我们的随机权重和偏差倾向于高于 50%准确率的原因。如果数据极度不平衡(例如 99%的数据属于一类),请记住使用混淆矩阵来跟踪假阳性和假阴性。
到目前为止,结构上一切都合理吗?在继续之前,随时可以回顾一切。我们只剩下最后一个部分要讲解:优化权重和偏差。去冲杯浓缩咖啡或氮气咖啡吧,因为这是本书中我们将要做的最复杂的数学!
反向传播
在我们开始使用随机梯度下降来优化我们的神经网络之前,我们面临的挑战是如何相应地改变每个权重和偏差值,尽管它们都纠缠在一起以创建输出变量,然后用于计算残差。我们如何找到每个权重 和偏差 变量的导数?我们需要使用我们在第一章中讨论过的链式法则。
计算权重和偏差的导数
我们还没有准备好应用随机梯度下降来训练我们的神经网络。我们需要求出相对于权重 和偏差 的偏导数,而且我们有链式法则来帮助我们。
虽然过程基本相同,但在神经网络上使用随机梯度下降存在一个复杂性。一层中的节点将它们的权重和偏置传递到下一层,然后应用另一组权重和偏置。这创建了一个类似洋葱的嵌套结构,我们需要从输出层开始解开。
在梯度下降过程中,我们需要找出应该调整哪些权重和偏置,以及调整多少,以减少整体成本函数。单个预测的成本将是神经网络的平方输出减去实际值:
但让我们再往下一层。那个激活输出只是带有激活函数的:
转而,是应用于激活输出的输出权重和偏置,这些输出来自隐藏层:
是由通过 ReLU 激活函数传递的构建而成:
最后,是由隐藏层加权和偏置的输入 x 值:
我们需要找到包含在矩阵和向量,,和中的权重和偏置,以最小化我们的损失。通过微调它们的斜率,我们可以改变对最小化损失影响最大的权重和偏置。然而,对权重或偏置进行微小调整会一直传播到外层的损失函数。这就是链式法则可以帮助我们找出这种影响的地方。
让我们专注于找到输出层权重和成本函数之间的关系。权重的变化导致未激活的输出的变化。然后改变激活输出,进而改变成本函数。利用链式法则,我们可以定义关于的导数如下:
当我们将这三个梯度相乘在一起时,我们得到了改变将如何改变成本函数的度量。
现在我们将计算这三个导数。让我们使用 SymPy 计算在示例 7-7 中关于的成本函数的导数。
示例 7-7. 计算成本函数关于的导数
from sympy import *
A2, y = symbols('A2 Y')
C = (A2 - Y)**2
dC_dA2 = diff(C, A2)
print(dC_dA2) # 2*A2 - 2*Y
接下来,让我们找到关于的导数与的导数(示例 7-8)。记住是激活函数的输出,在这种情况下是逻辑函数。因此,我们实际上只是在求取 S 形曲线的导数。
示例 7-8. 找到关于的导数与
from sympy import *
Z2 = symbols('Z2')
logistic = lambda x: 1 / (1 + exp(-x))
A2 = logistic(Z2)
dA2_dZ2 = diff(A2, Z2)
print(dA2_dZ2) # exp(-Z2)/(1 + exp(-Z2))**2
关于的导数将会得到,因为它只是一个线性函数,将返回斜率(示例 7-9)。
示例 7-9. 关于的导数
from sympy import *
A1, W2, B2 = symbols('A1, W2, B2')
Z2 = A1*W2 + B2
dZ2_dW2 = diff(Z2, W2)
print(dZ2_dW2) # A1
将所有内容整合在一起,这里是找到改变中的权重会如何影响成本函数的导数:
当我们运行一个输入X
与三个输入 R、G 和 B 值时,我们将得到、、和的值。
不要在数学中迷失!
在这一点很容易在数学中迷失并忘记你最初想要实现的目标,即找到成本函数关于输出层中的权重()的导数。当你发现自己陷入困境并忘记你要做什么时,那就退一步,出去走走,喝杯咖啡,提醒自己你要达成的目标。如果你做不到,你应该从头开始,一步步地找回迷失的地方。
然而,这只是神经网络的一个组成部分,对于的导数。以下是我们在示例 7-10 中使用 SymPy 计算的其余部分偏导数,这些是我们在链式求导中需要的。
示例 7-10。计算我们神经网络所需的所有偏导数
from sympy import *
W1, W2, B1, B2, A1, A2, Z1, Z2, X, Y = \
symbols('W1 W2 B1 B2 A1 A2 Z1 Z2 X Y')
# Calculate derivative of cost function with respect to A2
C = (A2 - Y)**2
dC_dA2 = diff(C, A2)
print("dC_dA2 = ", dC_dA2) # 2*A2 - 2*Y
# Calculate derivative of A2 with respect to Z2
logistic = lambda x: 1 / (1 + exp(-x))
_A2 = logistic(Z2)
dA2_dZ2 = diff(_A2, Z2)
print("dA2_dZ2 = ", dA2_dZ2) # exp(-Z2)/(1 + exp(-Z2))**2
# Calculate derivative of Z2 with respect to A1
_Z2 = A1*W2 + B2
dZ2_dA1 = diff(_Z2, A1)
print("dZ2_dA1 = ", dZ2_dA1) # W2
# Calculate derivative of Z2 with respect to W2
dZ2_dW2 = diff(_Z2, W2)
print("dZ2_dW2 = ", dZ2_dW2) # A1
# Calculate derivative of Z2 with respect to B2
dZ2_dB2 = diff(_Z2, B2)
print("dZ2_dB2 = ", dZ2_dB2) # 1
# Calculate derivative of A1 with respect to Z1
relu = lambda x: Max(x, 0)
_A1 = relu(Z1)
d_relu = lambda x: x > 0 # Slope is 1 if positive, 0 if negative
dA1_dZ1 = d_relu(Z1)
print("dA1_dZ1 = ", dA1_dZ1) # Z1 > 0
# Calculate derivative of Z1 with respect to W1
_Z1 = X*W1 + B1
dZ1_dW1 = diff(_Z1, W1)
print("dZ1_dW1 = ", dZ1_dW1) # X
# Calculate derivative of Z1 with respect to B1
dZ1_dB1 = diff(_Z1, B1)
print("dZ1_dB1 = ", dZ1_dB1) # 1
注意,ReLU 是手动计算的,而不是使用 SymPy 的diff()
函数。这是因为导数适用于平滑曲线,而不是 ReLU 上存在的锯齿角。但通过简单地声明斜率为正数为 1,负数为 0,可以轻松解决这个问题。这是有道理的,因为负数具有斜率为 0 的平坦线。但正数保持不变,具有 1:1 的斜率。
这些偏导数可以链接在一起,以创建相对于权重和偏置的新偏导数。让我们为、、和相对于成本函数的四个偏导数。我们已经讨论了。让我们将其与其他三个链式求导一起展示: