从零开始的数据科学第二版（二）

原文：zh.annas-archive.org/md5/48ab308fc34189a6d7d26b91b72a6df9

译者：飞龙

协议：CC BY-NC-SA 4.0

第三章：数据可视化

我相信可视化是实现个人目标最强大的手段之一。

Harvey Mackay

数据科学家工具包的基本组成部分是数据可视化。虽然创建可视化非常容易，但要生成优质的可视化要困难得多。

数据可视化有两个主要用途：

• 探索数据

• 传达数据

在本章中，我们将集中于构建您开始探索自己数据所需的技能，并生成我们在本书其余部分中将使用的可视化内容。与大多数章节主题一样，数据可视化是一个值得拥有自己书籍的丰富领域。尽管如此，我将尝试给您一些关于好的可视化和不好的可视化的认识。

matplotlib

存在各种工具用于数据可视化。我们将使用广泛使用的matplotlib 库，尽管它有些显得老旧，但仍被广泛使用。如果您有兴趣生成精致的互动可视化内容，那么它可能不是正确的选择，但对于简单的柱状图、折线图和散点图，它运行得相当不错。

正如前面提到的，matplotlib 不是 Python 核心库的一部分。在激活虚拟环境后（要设置一个，请返回到“虚拟环境”并按照说明操作），使用以下命令安装它：

python -m pip install matplotlib

我们将使用matplotlib.pyplot模块。在最简单的用法中，pyplot保持一个内部状态，您可以逐步构建可视化内容。完成后，您可以使用savefig保存它，或者使用show显示它。

例如，制作简单的图表（如图 3-1）非常简单：

from matplotlib import pyplot as plt

years = [1950, 1960, 1970, 1980, 1990, 2000, 2010]
gdp = [300.2, 543.3, 1075.9, 2862.5, 5979.6, 10289.7, 14958.3]

# create a line chart, years on x-axis, gdp on y-axis
plt.plot(years, gdp, color='green', marker='o', linestyle='solid')

# add a title
plt.title("Nominal GDP")

# add a label to the y-axis
plt.ylabel("Billions of $")
plt.show()

图 3-1. 一个简单的折线图

制作看起来像出版物质量的图表更为复杂，超出了本章的范围。您可以通过诸如轴标签、线条样式和点标记等方式自定义图表的许多方面。与其尝试全面处理这些选项，我们将在示例中简单使用（并强调）其中一些。

注

尽管我们不会使用这种功能的大部分，但 matplotlib 能够在图形中生成复杂的细节图、高级格式化和交互式可视化内容。如果您希望比我们在本书中更深入地了解，请查阅其文档。

柱状图

当您想展示某些数量在一些离散项目中如何变化时，柱状图是一个不错的选择。例如，图 3-2 显示了各种电影赢得的学院奖数量：

movies = ["Annie Hall", "Ben-Hur", "Casablanca", "Gandhi", "West Side Story"]
num_oscars = [5, 11, 3, 8, 10]

# plot bars with left x-coordinates [0, 1, 2, 3, 4], heights [num_oscars]
plt.bar(range(len(movies)), num_oscars)

plt.title("My Favorite Movies")     # add a title
plt.ylabel("# of Academy Awards")   # label the y-axis

# label x-axis with movie names at bar centers
plt.xticks(range(len(movies)), movies)

plt.show()

图 3-2. 一个简单的柱状图

柱状图也可以是绘制分桶数值直方图的好选择，如图 3-3，以便直观地探索数值的分布：

from collections import Counter
grades = [83, 95, 91, 87, 70, 0, 85, 82, 100, 67, 73, 77, 0]

# Bucket grades by decile, but put 100 in with the 90s
histogram = Counter(min(grade // 10 * 10, 90) for grade in grades)

plt.bar([x + 5 for x in histogram.keys()],  # Shift bars right by 5
        histogram.values(),                 # Give each bar its correct height
        10,                                 # Give each bar a width of 10
        edgecolor=(0, 0, 0))                # Black edges for each bar

plt.axis([-5, 105, 0, 5])                  # x-axis from -5 to 105,
                                           # y-axis from 0 to 5

plt.xticks([10 * i for i in range(11)])    # x-axis labels at 0, 10, ..., 100
plt.xlabel("Decile")
plt.ylabel("# of Students")
plt.title("Distribution of Exam 1 Grades")
plt.show()

图 3-3. 使用条形图制作直方图

对 plt.bar 的第三个参数指定了条形宽度。在这里，我们选择宽度为 10，以填满整个十分位。我们还将条形向右移动了 5，因此，例如，“10” 条形（对应于十分位 10-20）其中心在 15 处，因此占据了正确的范围。我们还为每个条形添加了黑色边缘，使它们在视觉上更加突出。

对 plt.axis 的调用指示我们希望 x 轴范围从 -5 到 105（仅留出一点左右的空间），y 轴应该从 0 到 5。对 plt.xticks 的调用在 0、10、20、…、100 处放置 x 轴标签。

使用 plt.axis 时要慎重。在创建条形图时，如果你的 y 轴不从 0 开始，这被认为是一种特别糟糕的做法，因为这很容易误导人（图 3-4）：

mentions = [500, 505]
years = [2017, 2018]

plt.bar(years, mentions, 0.8)
plt.xticks(years)
plt.ylabel("# of times I heard someone say 'data science'")

# if you don't do this, matplotlib will label the x-axis 0, 1
# and then add a +2.013e3 off in the corner (bad matplotlib!)
plt.ticklabel_format(useOffset=False)

# misleading y-axis only shows the part above 500
plt.axis([2016.5, 2018.5, 499, 506])
plt.title("Look at the 'Huge' Increase!")
plt.show()

图 3-4. 具有误导性 y 轴的图表

在 图 3-5 中，我们使用更合理的轴，效果看起来不那么令人印象深刻：

plt.axis([2016.5, 2018.5, 0, 550])
plt.title("Not So Huge Anymore")
plt.show()

图 3-5. 具有非误导性 y 轴的相同图表

折线图

正如我们已经看到的，我们可以使用 plt.plot 制作折线图。这些图表适合展示趋势，如 图 3-6 所示：

variance     = [1, 2, 4, 8, 16, 32, 64, 128, 256]
bias_squared = [256, 128, 64, 32, 16, 8, 4, 2, 1]
total_error  = [x + y for x, y in zip(variance, bias_squared)]
xs = [i for i, _ in enumerate(variance)]

# We can make multiple calls to plt.plot
# to show multiple series on the same chart
plt.plot(xs, variance,     'g-',  label='variance')    # green solid line
plt.plot(xs, bias_squared, 'r-.', label='bias²')      # red dot-dashed line
plt.plot(xs, total_error,  'b:',  label='total error') # blue dotted line

# Because we've assigned labels to each series,
# we can get a legend for free (loc=9 means "top center")
plt.legend(loc=9)
plt.xlabel("model complexity")
plt.xticks([])
plt.title("The Bias-Variance Tradeoff")
plt.show()

图 3-6. 带有图例的多条线图

散点图

散点图是可视化两组配对数据之间关系的正确选择。例如，图 3-7 描述了用户拥有的朋友数量与他们每天在网站上花费的时间之间的关系：

friends = [ 70,  65,  72,  63,  71,  64,  60,  64,  67]
minutes = [175, 170, 205, 120, 220, 130, 105, 145, 190]
labels =  ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i']

plt.scatter(friends, minutes)

# label each point
for label, friend_count, minute_count in zip(labels, friends, minutes):
    plt.annotate(label,
        xy=(friend_count, minute_count), # Put the label with its point
        xytext=(5, -5),                  # but slightly offset
        textcoords='offset points')

plt.title("Daily Minutes vs. Number of Friends")
plt.xlabel("# of friends")
plt.ylabel("daily minutes spent on the site")
plt.show()

图 3-7. 朋友和网站上时间的散点图

如果你散布可比较的变量，如果让 matplotlib 自行选择比例，可能会得到一个误导性的图像，如 图 3-8 所示。

图 3-8. 具有不可比较轴的散点图

test_1_grades = [ 99, 90, 85, 97, 80]
test_2_grades = [100, 85, 60, 90, 70]

plt.scatter(test_1_grades, test_2_grades)
plt.title("Axes Aren't Comparable")
plt.xlabel("test 1 grade")
plt.ylabel("test 2 grade")
plt.show()

如果我们包含对 plt.axis("equal") 的调用，则绘图（图 3-9）更准确地显示了大部分变化发生在测试 2 上。

这已足够让你开始进行可视化了。我们会在本书中更多地学习可视化。

图 3-9. 具有相等轴的相同散点图

进一步探索

• matplotlib Gallery 将为你展示使用 matplotlib 可以做的各种事情（以及如何做）。

• seaborn 建立在 matplotlib 之上，使你能够轻松生成更漂亮（和更复杂）的可视化。

• Altair 是一个较新的 Python 库，用于创建声明式可视化。

• D3.js 是用于在 Web 上生成复杂交互式可视化的 JavaScript 库。虽然它不是 Python 的库，但它被广泛使用，你值得熟悉它。

• Bokeh 是一个将 D3 风格的可视化引入 Python 的库。

第四章：线性代数

还有比代数更无用或更少有用的东西吗？

比利·康诺利

线性代数是处理向量空间的数学分支。虽然我不能指望在简短的章节中教给您线性代数，但它是许多数据科学概念和技术的基础，这意味着我至少应该尝试一下。我们在本章学到的内容将在本书的其余部分中广泛使用。

向量

抽象地说，向量是可以相加以形成新向量的对象，并且可以乘以标量（即数字），也可以形成新向量。

具体来说（对于我们而言），向量是某个有限维空间中的点。尽管您可能不认为自己的数据是向量，但它们通常是表示数值数据的有用方式。

例如，如果您有许多人的身高、体重和年龄数据，您可以将您的数据视为三维向量[height, weight, age]。如果您正在教授一个有四次考试的班级，您可以将学生的成绩视为四维向量[exam1, exam2, exam3, exam4]。

从头开始的最简单方法是将向量表示为数字列表。三个数字的列表对应于三维空间中的一个向量，反之亦然。

我们将通过一个类型别名来实现，即 Vector 只是 float 的 list：

from typing import List

Vector = List[float]

height_weight_age = [70,  # inches,
                     170, # pounds,
                     40 ] # years

grades = [95,   # exam1
          80,   # exam2
          75,   # exam3
          62 ]  # exam4

我们还希望在向量上执行算术运算。因为 Python 的 list 不是向量（因此不提供向量算术的工具），我们需要自己构建这些算术工具。所以让我们从这里开始。

首先，我们经常需要添加两个向量。向量是分量相加的。这意味着如果两个向量v和w长度相同，则它们的和就是一个向量，其第一个元素是v[0] + w[0]，第二个元素是v[1] + w[1]，依此类推。（如果它们长度不同，则不允许将它们相加。）

例如，将向量 [1, 2] 和 [2, 1] 相加的结果是 [1 + 2, 2 + 1] 或 [3, 3]，如图 4-1 所示。

图 4-1. 添加两个向量

我们可以通过使用 zip 将向量一起并使用列表推导来添加相应元素来轻松实现这一点：

def add(v: Vector, w: Vector) -> Vector:
    """Adds corresponding elements"""
    assert len(v) == len(w), "vectors must be the same length"

    return [v_i + w_i for v_i, w_i in zip(v, w)]

assert add([1, 2, 3], [4, 5, 6]) == [5, 7, 9]

类似地，要减去两个向量，我们只需减去相应的元素：

def subtract(v: Vector, w: Vector) -> Vector:
    """Subtracts corresponding elements"""
    assert len(v) == len(w), "vectors must be the same length"

    return [v_i - w_i for v_i, w_i in zip(v, w)]

assert subtract([5, 7, 9], [4, 5, 6]) == [1, 2, 3]

我们有时还希望对向量列表进行分量求和，即创建一个新向量，其第一个元素是所有第一个元素的和，第二个元素是所有第二个元素的和，依此类推：

def vector_sum(vectors: List[Vector]) -> Vector:
    """Sums all corresponding elements"""
    # Check that vectors is not empty
    assert vectors, "no vectors provided!"

    # Check the vectors are all the same size
    num_elements = len(vectors[0])
    assert all(len(v) == num_elements for v in vectors), "different sizes!"

    # the i-th element of the result is the sum of every vector[i]
    return [sum(vector[i] for vector in vectors)
            for i in range(num_elements)]

assert vector_sum([[1, 2], [3, 4], [5, 6], [7, 8]]) == [16, 20]

我们还需要能够将向量乘以标量，这可以通过将向量的每个元素乘以该数来简单实现：

def scalar_multiply(c: float, v: Vector) -> Vector:
    """Multiplies every element by c"""
    return [c * v_i for v_i in v]

assert scalar_multiply(2, [1, 2, 3]) == [2, 4, 6]

这使我们能够计算（相同大小的）向量列表的分量均值：

def vector_mean(vectors: List[Vector]) -> Vector:
    """Computes the element-wise average"""
    n = len(vectors)
    return scalar_multiply(1/n, vector_sum(vectors))

assert vector_mean([[1, 2], [3, 4], [5, 6]]) == [3, 4]

一个不那么明显的工具是点积。两个向量的点积是它们各分量的乘积之和：

def dot(v: Vector, w: Vector) -> float:
    """Computes v_1 * w_1 + ... + v_n * w_n"""
    assert len(v) == len(w), "vectors must be same length"

    return sum(v_i * w_i for v_i, w_i in zip(v, w))

assert dot([1, 2, 3], [4, 5, 6]) == 32  # 1 * 4 + 2 * 5 + 3 * 6

如果w的大小为 1，则点积测量向量v在w方向延伸的距离。例如，如果w = [1, 0]，那么dot(v, w)就是v的第一个分量。另一种说法是，这是您在将 v 投影到 w上时获得的向量的长度（图 4-2）。

图 4-2。向量投影的点积

使用这个方法，计算一个向量的平方和很容易：

def sum_of_squares(v: Vector) -> float:
    """Returns v_1 * v_1 + ... + v_n * v_n"""
    return dot(v, v)

assert sum_of_squares([1, 2, 3]) == 14  # 1 * 1 + 2 * 2 + 3 * 3

我们可以用它来计算其大小（或长度）：

import math

def magnitude(v: Vector) -> float:
    """Returns the magnitude (or length) of v"""
    return math.sqrt(sum_of_squares(v))   # math.sqrt is square root function

assert magnitude([3, 4]) == 5

我们现在已经有了计算两个向量之间距离的所有要素，定义为：

(v 1 -w 1 ) 2 + ... + (v n -w n ) 2

在代码中：

def squared_distance(v: Vector, w: Vector) -> float:
    """Computes (v_1 - w_1) ** 2 + ... + (v_n - w_n) ** 2"""
    return sum_of_squares(subtract(v, w))

def distance(v: Vector, w: Vector) -> float:
    """Computes the distance between v and w"""
    return math.sqrt(squared_distance(v, w))

如果我们将其写为（等价的）形式，则可能更清楚：

def distance(v: Vector, w: Vector) -> float:
    return magnitude(subtract(v, w))

这应该足以让我们开始了。我们将在整本书中大量使用这些函数。

注意

将列表用作向量在阐述上很好，但在性能上很差。

在生产代码中，您会希望使用 NumPy 库，该库包括一个高性能的数组类，其中包括各种算术运算。

矩阵

矩阵是一种二维数字集合。我们将矩阵表示为列表的列表，每个内部列表具有相同的大小，并表示矩阵的一行。如果A是一个矩阵，那么A[i][j]是矩阵的第i行和第j列的元素。按照数学约定，我们经常使用大写字母表示矩阵。例如：

# Another type alias
Matrix = List[List[float]]

A = [[1, 2, 3],  # A has 2 rows and 3 columns
     [4, 5, 6]]

B = [[1, 2],     # B has 3 rows and 2 columns
     [3, 4],
     [5, 6]]

注意

在数学上，您通常会将矩阵的第一行命名为“第 1 行”，第一列命名为“第 1 列”。因为我们用 Python list表示矩阵，而 Python 中的列表是从零开始索引的，所以我们将矩阵的第一行称为“第 0 行”，第一列称为“第 0 列”。

给定这个列表的列表表示，矩阵A有len(A)行和len(A[0])列，我们将其称为其形状：

from typing import Tuple

def shape(A: Matrix) -> Tuple[int, int]:
    """Returns (# of rows of A, # of columns of A)"""
    num_rows = len(A)
    num_cols = len(A[0]) if A else 0   # number of elements in first row
    return num_rows, num_cols

assert shape([[1, 2, 3], [4, 5, 6]]) == (2, 3)  # 2 rows, 3 columns

如果一个矩阵有n行和k列，我们将其称为n × k 矩阵。我们可以（有时会）将n × k矩阵的每一行视为长度为k的向量，将每一列视为长度为n的向量：

def get_row(A: Matrix, i: int) -> Vector:
    """Returns the i-th row of A (as a Vector)"""
    return A[i]             # A[i] is already the ith row

def get_column(A: Matrix, j: int) -> Vector:
    """Returns the j-th column of A (as a Vector)"""
    return [A_i[j]          # jth element of row A_i
            for A_i in A]   # for each row A_i

我们还希望能够根据其形状和生成其元素的函数创建一个矩阵。我们可以使用嵌套列表推导来实现这一点：

from typing import Callable

def make_matrix(num_rows: int,
                num_cols: int,
                entry_fn: Callable[[int, int], float]) -> Matrix:
    """
 Returns a num_rows x num_cols matrix
 whose (i,j)-th entry is entry_fn(i, j)
 """
    return [[entry_fn(i, j)             # given i, create a list
             for j in range(num_cols)]  #   [entry_fn(i, 0), ... ]
            for i in range(num_rows)]   # create one list for each i

给定此函数，您可以制作一个 5 × 5 单位矩阵（对角线上为 1，其他位置为 0），如下所示：

def identity_matrix(n: int) -> Matrix:
    """Returns the n x n identity matrix"""
    return make_matrix(n, n, lambda i, j: 1 if i == j else 0)

assert identity_matrix(5) == [[1, 0, 0, 0, 0],
                              [0, 1, 0, 0, 0],
                              [0, 0, 1, 0, 0],
                              [0, 0, 0, 1, 0],
                              [0, 0, 0, 0, 1]]

矩阵对我们来说将是重要的，有几个原因。

首先，我们可以使用矩阵来表示由多个向量组成的数据集，简单地将每个向量视为矩阵的一行。例如，如果您有 1,000 人的身高、体重和年龄，您可以将它们放在一个 1,000 × 3 矩阵中：

data = [[70, 170, 40],
        [65, 120, 26],
        [77, 250, 19],
        # ....
       ]

其次，正如我们后面将看到的，我们可以使用n × k矩阵来表示将k维向量映射到n维向量的线性函数。我们的许多技术和概念都涉及这样的函数。

第三，矩阵可用于表示二元关系。在 第一章 中，我们将网络的边表示为一组对 (i, j)。另一种表示方法是创建一个矩阵 A，使得如果节点 i 和 j 连接，则 A[i][j] 为 1，否则为 0。

请记住，在此之前我们有：

friendships = [(0, 1), (0, 2), (1, 2), (1, 3), (2, 3), (3, 4),
               (4, 5), (5, 6), (5, 7), (6, 8), (7, 8), (8, 9)]

我们也可以将这表示为：

#            user 0  1  2  3  4  5  6  7  8  9
#
friend_matrix = [[0, 1, 1, 0, 0, 0, 0, 0, 0, 0],  # user 0
                 [1, 0, 1, 1, 0, 0, 0, 0, 0, 0],  # user 1
                 [1, 1, 0, 1, 0, 0, 0, 0, 0, 0],  # user 2
                 [0, 1, 1, 0, 1, 0, 0, 0, 0, 0],  # user 3
                 [0, 0, 0, 1, 0, 1, 0, 0, 0, 0],  # user 4
                 [0, 0, 0, 0, 1, 0, 1, 1, 0, 0],  # user 5
                 [0, 0, 0, 0, 0, 1, 0, 0, 1, 0],  # user 6
                 [0, 0, 0, 0, 0, 1, 0, 0, 1, 0],  # user 7
                 [0, 0, 0, 0, 0, 0, 1, 1, 0, 1],  # user 8
                 [0, 0, 0, 0, 0, 0, 0, 0, 1, 0]]  # user 9

如果连接很少，这种表示方法会更加低效，因为你最终需要存储很多个零。然而，使用矩阵表示法可以更快地检查两个节点是否连接——你只需要进行矩阵查找，而不是（可能）检查每条边：

assert friend_matrix[0][2] == 1, "0 and 2 are friends"
assert friend_matrix[0][8] == 0, "0 and 8 are not friends"

同样，要找到一个节点的连接，你只需要检查对应节点的列（或行）：

# only need to look at one row
friends_of_five = [i
                   for i, is_friend in enumerate(friend_matrix[5])
                   if is_friend]

对于一个小图，你可以简单地为每个节点对象添加一个连接列表以加速这个过程；但对于一个大型、不断演化的图，这可能会太昂贵和难以维护。

我们将在本书中多次涉及矩阵。

进一步探索

• 线性代数被数据科学家广泛使用（通常是隐式地使用，并且经常被不理解它的人使用）。阅读一本教材并不是一个坏主意。你可以在网上找到几本免费的教材：

  • 线性代数，作者 Jim Hefferon（圣迈克尔学院）

  • 线性代数，作者 David Cherney、Tom Denton、Rohit Thomas 和 Andrew Waldron（加州大学戴维斯分校）

  • 如果你感觉有冒险精神，线性代数错误方法，作者 Sergei Treil（布朗大学），是一个更高级的入门教程。

• 如果你使用 NumPy，那么我们在本章中构建的所有机制都是免费的。（你还会得到更多，包括更好的性能。）

第五章：统计学

事实是倔强的，但统计数据更加灵活。

马克·吐温

统计学是指我们理解数据的数学和技术。这是一个丰富而庞大的领域，更适合放在图书馆的书架（或房间）上，而不是书中的一章，因此我们的讨论必然不会很深入。相反，我会尽量教你足够的知识，让你有点危险，并激起你的兴趣，让你去学更多的东西。

描述单一数据集

通过口口相传和一点运气，DataSciencester 已经发展成了数十个成员，筹款副总裁要求你提供一些关于你的成员拥有多少朋友的描述，以便他可以在提出自己的想法时加入其中。

使用第一章中的技术，你很容易就能产生这些数据。但现在你面临着如何描述它的问题。

任何数据集的一个明显描述就是数据本身：

num_friends = [100, 49, 41, 40, 25,
               # ... and lots more
              ]

对于足够小的数据集，这甚至可能是最好的描述。但对于较大的数据集来说，这是笨重且可能晦涩的。（想象一下盯着一百万个数字的列表。）因此，我们使用统计数据来提炼和传达我们数据的相关特征。

作为第一步，你可以使用 Counter 和 plt.bar 将朋友数量制作成直方图（参见图 5-1）：

from collections import Counter
import matplotlib.pyplot as plt

friend_counts = Counter(num_friends)
xs = range(101)                         # largest value is 100
ys = [friend_counts[x] for x in xs]     # height is just # of friends
plt.bar(xs, ys)
plt.axis([0, 101, 0, 25])
plt.title("Histogram of Friend Counts")
plt.xlabel("# of friends")
plt.ylabel("# of people")
plt.show()

图 5-1 朋友数量直方图

不幸的是，这张图表仍然太难以插入对话中。所以你开始生成一些统计数据。可能最简单的统计数据是数据点的数量：

num_points = len(num_friends)               # 204

你可能也对最大值和最小值感兴趣：

largest_value = max(num_friends)            # 100
smallest_value = min(num_friends)           # 1

这些只是想要知道特定位置的值的特例：

sorted_values = sorted(num_friends)
smallest_value = sorted_values[0]           # 1
second_smallest_value = sorted_values[1]    # 1
second_largest_value = sorted_values[-2]    # 49

但我们只是刚刚开始。

中心趋势

通常，我们会想知道我们的数据位于何处。最常见的是我们会使用均值（或平均值），它只是数据之和除以其计数：

def mean(xs: List[float]) -> float:
    return sum(xs) / len(xs)

mean(num_friends)   # 7.333333

如果你有两个数据点，平均值就是它们之间的中点。随着你增加更多的点，平均值会移动，但它总是依赖于每个点的值。例如，如果你有 10 个数据点，而你增加其中任何一个的值 1，你就会增加平均值 0.1。

我们有时也会对中位数感兴趣，它是中间最值（如果数据点数量为奇数）或两个中间最值的平均值（如果数据点数量为偶数）。

例如，如果我们在一个排序好的向量 x 中有五个数据点，中位数是 x[5 // 2] 或 x[2]。如果我们有六个数据点，我们想要 x[2]（第三个点）和 x[3]（第四个点）的平均值。

注意——与平均值不同——中位数并不完全依赖于你的数据中的每个值。例如，如果你使最大值变大（或最小值变小），中间值保持不变，这意味着中位数也保持不变。

我们将为偶数和奇数情况编写不同的函数并将它们结合起来：

# The underscores indicate that these are "private" functions, as they're
# intended to be called by our median function but not by other people
# using our statistics library.
def _median_odd(xs: List[float]) -> float:
    """If len(xs) is odd, the median is the middle element"""
    return sorted(xs)[len(xs) // 2]

def _median_even(xs: List[float]) -> float:
    """If len(xs) is even, it's the average of the middle two elements"""
    sorted_xs = sorted(xs)
    hi_midpoint = len(xs) // 2  # e.g. length 4 => hi_midpoint 2
    return (sorted_xs[hi_midpoint - 1] + sorted_xs[hi_midpoint]) / 2

def median(v: List[float]) -> float:
    """Finds the 'middle-most' value of v"""
    return _median_even(v) if len(v) % 2 == 0 else _median_odd(v)

assert median([1, 10, 2, 9, 5]) == 5
assert median([1, 9, 2, 10]) == (2 + 9) / 2

现在我们可以计算朋友的中位数了：

print(median(num_friends))  # 6

显然，计算均值更简单，并且随着数据的变化而平稳变化。如果我们有n个数据点，其中一个增加了一小部分e，那么均值将增加e / n。（这使得均值适合各种微积分技巧。）然而，为了找到中位数，我们必须对数据进行排序。并且改变数据点中的一个小量e可能会使中位数增加e，少于e的某个数字，或者根本不增加（这取决于其他数据）。

注意

实际上，有一些不明显的技巧可以高效地计算中位数，而无需对数据进行排序。然而，这超出了本书的范围，所以我们必须对数据进行排序。

同时，均值对我们数据中的异常值非常敏感。如果我们最友好的用户有 200 个朋友（而不是 100 个），那么均值将上升到 7.82，而中位数将保持不变。如果异常值可能是不良数据（或者在其他方面不代表我们试图理解的现象），那么均值有时可能会给我们提供错误的图片。例如，通常会讲述这样一个故事：上世纪 80 年代中期，北卡罗来纳大学的起薪最高的专业是地理学，主要是因为 NBA 球星（及异常值）迈克尔·乔丹。

中位数的一种推广是分位数，它表示数据中位于某个百分位以下的值（中位数表示 50%的数据位于其以下）：

def quantile(xs: List[float], p: float) -> float:
    """Returns the pth-percentile value in x"""
    p_index = int(p * len(xs))
    return sorted(xs)[p_index]

assert quantile(num_friends, 0.10) == 1
assert quantile(num_friends, 0.25) == 3
assert quantile(num_friends, 0.75) == 9
assert quantile(num_friends, 0.90) == 13

更少见的是，您可能希望查看模式，或者最常见的值：

def mode(x: List[float]) -> List[float]:
    """Returns a list, since there might be more than one mode"""
    counts = Counter(x)
    max_count = max(counts.values())
    return [x_i for x_i, count in counts.items()
            if count == max_count]

assert set(mode(num_friends)) == {1, 6}

但是大多数情况下，我们会使用平均值。

离散度

分散是指我们的数据分散程度的度量。通常，这些统计量的值接近零表示完全不分散，而大值（无论其含义如何）则表示非常分散。例如，一个非常简单的度量是范围，它只是最大和最小元素之间的差异：

# "range" already means something in Python, so we'll use a different name
def data_range(xs: List[float]) -> float:
    return max(xs) - min(xs)

assert data_range(num_friends) == 99

当max和min相等时，范围恰好为零，这只可能发生在x的元素全部相同的情况下，这意味着数据尽可能不分散。相反，如果范围很大，则max远远大于min，数据更加分散。

像中位数一样，范围实际上并不依赖于整个数据集。一个数据集，其数据点全部为 0 或 100，其范围与一个数据集相同，其中数值为 0、100 和大量 50。但是第一个数据集看起来“应该”更加分散。

分散度的更复杂的测量是方差，其计算如下：

from scratch.linear_algebra import sum_of_squares

def de_mean(xs: List[float]) -> List[float]:
    """Translate xs by subtracting its mean (so the result has mean 0)"""
    x_bar = mean(xs)
    return [x - x_bar for x in xs]

def variance(xs: List[float]) -> float:
    """Almost the average squared deviation from the mean"""
    assert len(xs) >= 2, "variance requires at least two elements"

    n = len(xs)
    deviations = de_mean(xs)
    return sum_of_squares(deviations) / (n - 1)

assert 81.54 < variance(num_friends) < 81.55

注意

这看起来几乎是平均平方偏差，不过我们除以的是n - 1而不是n。事实上，当我们处理来自更大总体的样本时，x_bar只是实际均值的估计值，这意味着平均而言(x_i - x_bar) ** 2是x_i偏离均值的平方偏差的一个低估值，这就是为什么我们除以n - 1而不是n。参见Wikipedia。

现在，无论我们的数据以什么单位表示（例如，“朋友”），我们所有的集中趋势度量都是以相同的单位表示的。范围也将以相同的单位表示。另一方面，方差的单位是原始单位的平方（例如，“朋友的平方”）。由于这些单位可能难以理解，我们通常转而查看标准差：

import math

def standard_deviation(xs: List[float]) -> float:
    """The standard deviation is the square root of the variance"""
    return math.sqrt(variance(xs))

assert 9.02 < standard_deviation(num_friends) < 9.04

范围和标准差都有与我们之前对均值所看到的相同的异常值问题。以相同的例子，如果我们最友好的用户反而有 200 个朋友，标准差将为 14.89——高出 60%以上！

一种更加稳健的替代方法是计算第 75 百分位数值和第 25 百分位数值之间的差异：

def interquartile_range(xs: List[float]) -> float:
    """Returns the difference between the 75%-ile and the 25%-ile"""
    return quantile(xs, 0.75) - quantile(xs, 0.25)

assert interquartile_range(num_friends) == 6

这是一个显然不受少数异常值影响的数字。

相关性

DataSciencester 的增长副总裁有一个理论，即人们在网站上花费的时间与他们在网站上拥有的朋友数量有关（她不是无缘无故的副总裁），她请求您验证这一理论。

在分析了流量日志后，您得到了一个名为daily_minutes的列表，显示每个用户每天在 DataSciencester 上花费的分钟数，并且您已经按照前面的num_friends列表的顺序对其进行了排序。我们想要调查这两个指标之间的关系。

我们首先来看看协方差，它是方差的配对模拟。而方差衡量的是单个变量偏离其均值的程度，协方差衡量的是两个变量如何联动地偏离各自的均值：

from scratch.linear_algebra import dot

def covariance(xs: List[float], ys: List[float]) -> float:
    assert len(xs) == len(ys), "xs and ys must have same number of elements"

    return dot(de_mean(xs), de_mean(ys)) / (len(xs) - 1)

assert 22.42 < covariance(num_friends, daily_minutes) < 22.43
assert 22.42 / 60 < covariance(num_friends, daily_hours) < 22.43 / 60

请记住，dot求和了对应元素对的乘积。当x和y的对应元素都高于它们的均值或者都低于它们的均值时，一个正数进入总和。当一个高于其均值而另一个低于其均值时，一个负数进入总和。因此，“大”正协方差意味着x在y大时也大，在y小时也小。 “大”负协方差意味着相反的情况，即x在y大时较小，在y小时较大。接近零的协方差意味着没有这种关系存在。

尽管如此，由于一些原因，这个数字可能难以解释：

• 其单位是输入单位的乘积（例如，每天的朋友分钟），这可能难以理解。（“每天的朋友分钟”是什么意思？）

• 如果每个用户的朋友数量增加一倍（但花费的时间相同），协方差将增加一倍。但在某种意义上，变量之间的关联性将是一样的。换句话说，很难说什么才算是“大”协方差。

因此，更常见的是看相关性，它将两个变量的标准差除出：

def correlation(xs: List[float], ys: List[float]) -> float:
    """Measures how much xs and ys vary in tandem about their means"""
    stdev_x = standard_deviation(xs)
    stdev_y = standard_deviation(ys)
    if stdev_x > 0 and stdev_y > 0:
        return covariance(xs, ys) / stdev_x / stdev_y
    else:
        return 0    # if no variation, correlation is zero

assert 0.24 < correlation(num_friends, daily_minutes) < 0.25
assert 0.24 < correlation(num_friends, daily_hours) < 0.25

相关性是无单位的，并且总是介于–1（完全反相关）和 1（完全相关）之间。像 0.25 这样的数字代表相对较弱的正相关。

然而，我们忽略的一件事是检查我们的数据。查看图 5-2。

图 5-2. 含有离群值的相关性

拥有 100 个朋友（每天只在网站上花 1 分钟）的人是一个巨大的离群值，而相关性对离群值非常敏感。如果我们忽略他会发生什么？

outlier = num_friends.index(100)    # index of outlier

num_friends_good = [x
                    for i, x in enumerate(num_friends)
                    if i != outlier]

daily_minutes_good = [x
                      for i, x in enumerate(daily_minutes)
                      if i != outlier]

daily_hours_good = [dm / 60 for dm in daily_minutes_good]

assert 0.57 < correlation(num_friends_good, daily_minutes_good) < 0.58
assert 0.57 < correlation(num_friends_good, daily_hours_good) < 0.58

如果没有这个离群值，相关性会更加显著（图 5-3）。

图 5-3. 移除离群值后的相关性

你进一步调查并发现，这个离群值实际上是一个从未被删除的内部测试账户。因此，你认为排除它是合理的。

Simpson 悖论

在分析数据时，一个不常见的惊喜是Simpson 悖论，其中当忽略混淆变量时，相关性可能会误导。

例如，假设你可以将所有成员标识为东海岸数据科学家或西海岸数据科学家。你决定检查哪个海岸的数据科学家更友好：

海岸	成员数	平均朋友数
西海岸	101	8.2
东海岸	103	6.5

确实看起来西海岸的数据科学家比东海岸的数据科学家更友好。你的同事们对此提出了各种理论：也许是阳光、咖啡、有机农产品或者太平洋的悠闲氛围？

但在处理数据时，你发现了一些非常奇怪的事情。如果你只看有博士学位的人，东海岸的数据科学家平均有更多的朋友。如果你只看没有博士学位的人，东海岸的数据科学家平均也有更多的朋友！

海岸	学位	成员数	平均朋友数
西海岸	博士学位	35	3.1
东海岸	博士学位	70	3.2
西海岸	无博士学位	66	10.9
东海岸	无博士学位	33	13.4

一旦考虑到用户的学位，相关性的方向就会相反！将数据分桶为东海岸/西海岸掩盖了东海岸数据科学家更倾向于博士类型的事实。

这种现象在现实世界中经常发生。关键问题在于，相关性是在“其他一切相等”的情况下衡量你两个变量之间的关系。如果你的数据分类是随机分配的，如在一个设计良好的实验中可能发生的那样，“其他一切相等”可能不是一个糟糕的假设。但是，当类别分配有更深层次的模式时，“其他一切相等”可能是一个糟糕的假设。

避免这种情况的唯一真正方法是了解你的数据，并尽力确保你已经检查了可能的混淆因素。显然，这并不总是可能的。如果你没有这 200 名数据科学家的教育背景数据，你可能会简单地得出结论，西海岸有一些与社交性相关的因素。

其他一些相关性警告

相关系数为零表明两个变量之间没有线性关系。然而，可能存在其他类型的关系。例如，如果：

x = [-2, -1, 0, 1, 2]
y = [ 2,  1, 0, 1, 2]

那么x和y的相关性为零。但它们肯定有一种关系——y的每个元素等于x对应元素的绝对值。它们没有的是一种关系，即知道x_i与mean(x)的比较信息无法提供有关y_i与mean(y)的信息。这是相关性寻找的关系类型。

此外，相关性并不能告诉你关系的大小。以下变量：

x = [-2, -1, 0, 1, 2]
y = [99.98, 99.99, 100, 100.01, 100.02]

完全相关，但（根据你测量的内容）很可能这种关系并不那么有趣。

相关性与因果关系

你可能曾经听说过“相关不等于因果”，很可能是从某人那里听来的，他看到的数据挑战了他不愿质疑的世界观的某些部分。尽管如此，这是一个重要的观点——如果x和y强相关，这可能意味着x导致y，y导致x，彼此相互导致，第三个因素导致两者，或者根本不导致。

考虑num_friends和daily_minutes之间的关系。可能在这个网站上拥有更多的朋友导致DataSciencester 用户花更多时间在网站上。如果每个朋友每天发布一定数量的内容，这可能是情况，这意味着你拥有的朋友越多，就需要花费更多时间来了解他们的更新。

但是，在 DataSciencester 论坛上花费更多时间可能会导致用户遇到和结识志同道合的人。也就是说，花费更多时间在这个网站上导致用户拥有更多的朋友。

第三种可能性是，对数据科学最热衷的用户在网站上花费更多时间（因为他们觉得更有趣），并且更积极地收集数据科学朋友（因为他们不想与其他人交往）。

一种增强因果关系信心的方法是进行随机试验。如果你能随机将用户分为两组，这两组具有相似的人口统计学特征，并给其中一组提供稍有不同的体验，那么你通常可以相当确信不同的体验导致了不同的结果。

例如，如果你不介意被指责https://www.nytimes.com/2014/06/30/technology/facebook-tinkers-with-users-emotions-in-news-feed-experiment-stirring-outcry.html?r=0[对用户进行实验]，你可以随机选择你的一部分用户，并仅向他们展示部分朋友的内容。如果这部分用户随后在网站上花费的时间较少，那么这将使你相信拥有更多朋友 _ 会导致更多的时间花费在该网站上。

进一步探索

• SciPy，pandas，和StatsModels都提供各种各样的统计函数。

• 统计学是重要的。（或许统计是重要的？）如果你想成为更好的数据科学家，阅读统计学教材是个好主意。许多在线免费教材可供选择，包括：

  • Introductory Statistics，由道格拉斯·莎弗和张志义（Saylor Foundation）编写。

  • OnlineStatBook，由大卫·莱恩（莱斯大学）编写。

  • Introductory Statistics，由 OpenStax（OpenStax College）编写。

第六章：概率

概率定律，在一般情况下正确，在特定情况下谬误。

爱德华·吉本

没有对概率及其数学的某种理解，数据科学是很难做的。就像我们在第五章中处理统计学一样，我们将在很多地方简化和省略技术细节。

对于我们的目的，你应该把概率看作是对从某个事件宇宙中选择的事件的不确定性进行量化的一种方式。与其深究这些术语的技术含义，不如想象掷骰子。事件宇宙包含所有可能的结果。而这些结果的任何子集就是一个事件；例如，“骰子掷出 1”或“骰子掷出偶数”。

在符号上，我们写 P(E) 来表示“事件 E 发生的概率”。

我们将使用概率理论来建立模型。我们将使用概率理论来评估模型。我们将在很多地方使用概率理论。

如果愿意，可以深入探讨概率论的哲学含义。（最好在喝啤酒时进行。）我们不会做这件事。

依赖性与独立性

简而言之，如果知道 E 发生与否能够提供 F 发生与否的信息（反之亦然），我们称事件 E 和 F 是依赖的。否则，它们是独立的。

例如，如果我们抛一枚公平的硬币两次，知道第一次抛硬币是正面并不能提供关于第二次抛硬币是否正面的信息。这两个事件是独立的。另一方面，知道第一次抛硬币是正面肯定会影响到第二次抛硬币是否都是反面。（如果第一次抛硬币是正面，那么肯定不会是两次抛硬币都是反面。）这两个事件是依赖的。

数学上，我们说事件 E 和 F 是独立的，如果它们同时发生的概率等于它们各自发生的概率的乘积：

P ( E , F ) = P ( E ) P ( F )

在例子中，“第一次抛硬币为正面”的概率是 1/2，“两次抛硬币都为反面”的概率是 1/4，但“第一次抛硬币为正面且两次抛硬币都为反面”的概率是 0。

条件概率

当两个事件 E 和 F 是独立的时候，根据定义我们有：

P ( E , F ) = P ( E ) P ( F )

如果它们不一定是独立的（如果 F 的概率不为零），那么我们定义 E 在给定 F 的条件下的概率为：

P ( E | F ) = P ( E , F ) / P ( F )

你应该把这看作是在我们知道 F 发生的情况下，E 发生的概率。

我们经常将其重写为：

P ( E , F ) = P ( E | F ) P ( F )

当 E 和 F 是独立的时候，你可以检查这是否成立：

P ( E | F ) = P ( E )

这是数学上表达，即知道 F 发生了并不能给我们关于 E 是否发生额外的信息。

一个常见的棘手例子涉及一对（未知的）孩子的家庭。如果我们假设：

• 每个孩子都同等可能是男孩或女孩。

• 第二个孩子的性别与第一个孩子的性别是独立的。

那么事件“没有女孩”的概率为 1/4，事件“一个女孩一个男孩”的概率为 1/2，事件“两个女孩”的概率为 1/4。

现在我们可以问事件“两个孩子都是女孩”(B)在事件“老大是女孩”(G)条件下的概率是多少？使用条件概率的定义：

P ( B | G ) = P ( B , G ) / P ( G ) = P ( B ) / P ( G ) = 1 / 2

因为事件B和G（“两个孩子都是女孩且老大是女孩”）就是事件B。（一旦知道两个孩子都是女孩，老大是女孩就是必然的。）

大多数情况下，这个结果符合您的直觉。

我们还可以询问事件“两个孩子都是女孩”在事件“至少一个孩子是女孩”(L)条件下的概率。令人惊讶的是，答案与之前不同！

与之前一样，事件B和L（“两个孩子都是女孩且至少一个孩子是女孩”）就是事件B。这意味着我们有：

P ( B | L ) = P ( B , L ) / P ( L ) = P ( B ) / P ( L ) = 1 / 3

这是怎么回事？嗯，如果您所知道的是至少有一个孩子是女孩，那么这个家庭有一个男孩和一个女孩的可能性是有两倍多于两个女孩的可能性的。

我们可以通过“生成”许多家庭来验证这一点：

import enum, random

# An Enum is a typed set of enumerated values. We can use them
# to make our code more descriptive and readable.
class Kid(enum.Enum):
    BOY = 0
    GIRL = 1

def random_kid() -> Kid:
    return random.choice([Kid.BOY, Kid.GIRL])

both_girls = 0
older_girl = 0
either_girl = 0

random.seed(0)

for _ in range(10000):
    younger = random_kid()
    older = random_kid()
    if older == Kid.GIRL:
        older_girl += 1
    if older == Kid.GIRL and younger == Kid.GIRL:
        both_girls += 1
    if older == Kid.GIRL or younger == Kid.GIRL:
        either_girl += 1

print("P(both | older):", both_girls / older_girl)     # 0.514 ~ 1/2
print("P(both | either): ", both_girls / either_girl)  # 0.342 ~ 1/3

贝叶斯定理

数据科学家的最佳朋友之一是贝叶斯定理，这是一种“反转”条件概率的方法。假设我们需要知道某事件E在另一事件F发生条件下的概率。但我们只有关于E发生条件下F的概率信息。使用条件概率的定义两次告诉我们：

P ( E | F ) = P ( E , F ) / P ( F ) = P ( F | E ) P ( E ) / P ( F )

事件F可以分为两个互斥事件：“F和E”以及“F和非E”。如果我们用¬ E表示“非E”（即“E不发生”），那么：

P ( F ) = P ( F , E ) + P ( F , ¬ E )

所以：

P ( E | F ) = P ( F | E ) P ( E ) / [ P ( F | E ) P ( E ) + P ( F | ¬ E ) P ( ¬ E ) ]

这通常是贝叶斯定理的陈述方式。

这个定理经常被用来展示为什么数据科学家比医生更聪明。想象一种影响每 10,000 人中 1 人的特定疾病。再想象一种测试这种疾病的方法，它在 99%的情况下给出正确结果（如果您有疾病则为“患病”，如果您没有则为“未患病”）。

正测试意味着什么？让我们用T表示“您的测试呈阳性”事件，D表示“您患有疾病”事件。那么贝叶斯定理表明，在测试呈阳性的条件下，您患有疾病的概率是：

P ( D | T ) = P ( T | D ) P ( D ) / [ P ( T | D ) P ( D ) + P ( T | ¬ D ) P ( ¬ D ) ]

在这里我们知道，P ( T | D )，即染病者测试阳性的概率，为 0.99。P(D)，即任何给定人患病的概率，为 1/10,000 = 0.0001。 P ( T | ¬ D )，即没有患病者测试阳性的概率，为 0.01。而 P ( ¬ D )，即任何给定人没有患病的概率，为 0.9999。如果你把这些数字代入贝叶斯定理，你会发现：

P ( D | T ) = 0 . 98 %

也就是说，只有不到 1%的阳性测试者实际上患有这种疾病。

注

这假设人们更多或更少是随机参加测试的。如果只有具有某些症状的人参加测试，我们将需要在“阳性测试 和 症状”事件上进行条件判断，而阳性测试的人数可能会高得多。

更直观地看待这个问题的方法是想象一种 100 万人口的人群。你会预期其中有 100 人患有这种疾病，其中有 99 人测试结果呈阳性。另一方面，你会预期这中间有 999,900 人没有患有这种疾病，其中有 9,999 人测试结果呈阳性。这意味着你只会预期(99 + 9999)个阳性测试者中有 99 人实际上患有这种疾病。

随机变量

随机变量是具有相关概率分布的可能值的变量。一个非常简单的随机变量，如果抛硬币正面朝上则值为 1，如果反面朝上则值为 0。一个更复杂的随机变量可能会测量你在抛硬币 10 次时观察到的头像数量，或者从range(10)中选取的值，其中每个数字都是同样可能的。

相关的分布给出了变量实现其可能值的概率。抛硬币变量等于 0 的概率为 0.5，等于 1 的概率为 0.5。range(10)变量具有分布，将 0 到 9 中的每个数字分配概率 0.1。

有时我们会谈论一个随机变量的期望值，这是其值按其概率加权的平均值。抛硬币变量的期望值为 1/2（= 0 * 1/2 + 1 * 1/2），而range(10)变量的期望值为 4.5。

随机变量可以像其他事件一样被条件化。回到“条件概率”中的两个孩子示例，如果X是表示女孩数量的随机变量，那么X等于 0 的概率为 1/4，等于 1 的概率为 1/2，等于 2 的概率为 1/4。

我们可以定义一个新的随机变量Y，条件是至少有一个孩子是女孩。然后Y以 2/3 的概率等于 1，以 1/3 的概率等于 2。还有一个变量Z，条件是较大的孩子是女孩，则以 1/2 的概率等于 1，以 1/2 的概率等于 2。

大多数情况下，在我们进行的操作中，我们将隐含地使用随机变量，而不特别关注它们。但是如果你深入研究，你会发现它们。

连续分布

抛硬币对应于离散分布——它将离散结果与正概率关联起来。通常我们希望模拟跨越连续结果的分布。（对我们而言，这些结果将始终是实数，尽管在现实生活中并非总是如此。）例如，均匀分布将在 0 到 1 之间所有数字上赋予相等的权重。

因为在 0 和 1 之间有无穷多个数，这意味着它分配给单个点的权重必然为零。因此，我们用概率密度函数（PDF）表示连续分布，使得在某个区间内看到一个值的概率等于该区间上密度函数的积分。

注意

如果你的积分微积分有点生疏，一个更简单的理解方式是，如果一个分布有密度函数f，那么在x和x + h之间看到一个值的概率大约是h * f(x)，如果h很小的话。

均匀分布的密度函数就是：

def uniform_pdf(x: float) -> float:
    return 1 if 0 <= x < 1 else 0

按照该分布，随机变量落在 0.2 到 0.3 之间的概率是 1/10，正如你所期望的那样。Python 的random.random是一个均匀密度的（伪）随机变量。

我们通常更感兴趣的是累积分布函数（CDF），它给出随机变量小于或等于某个值的概率。对于均匀分布，创建 CDF 并不困难（见图 6-1）：

def uniform_cdf(x: float) -> float:
    """Returns the probability that a uniform random variable is <= x"""
    if x < 0:   return 0    # uniform random is never less than 0
    elif x < 1: return x    # e.g. P(X <= 0.4) = 0.4
    else:       return 1    # uniform random is always less than 1

图 6-1. 均匀分布的累积分布函数

正态分布

正态分布是经典的钟形曲线分布，完全由两个参数确定：其均值μ（mu）和标准差σ（sigma）。均值表示钟形曲线的中心位置，标准差表示其“宽度”。

它的 PDF 是：

f ( x | μ , σ ) = 1 2πσ exp ( - (x-μ) 2 2σ 2 )

我们可以这样实现它：

import math
SQRT_TWO_PI = math.sqrt(2 * math.pi)

def normal_pdf(x: float, mu: float = 0, sigma: float = 1) -> float:
    return (math.exp(-(x-mu) ** 2 / 2 / sigma ** 2) / (SQRT_TWO_PI * sigma))

在图 6-2 中，我们绘制了一些这些 PDF，看看它们的样子：

import matplotlib.pyplot as plt
xs = [x / 10.0 for x in range(-50, 50)]
plt.plot(xs,[normal_pdf(x,sigma=1) for x in xs],'-',label='mu=0,sigma=1')
plt.plot(xs,[normal_pdf(x,sigma=2) for x in xs],'--',label='mu=0,sigma=2')
plt.plot(xs,[normal_pdf(x,sigma=0.5) for x in xs],':',label='mu=0,sigma=0.5')
plt.plot(xs,[normal_pdf(x,mu=-1)   for x in xs],'-.',label='mu=-1,sigma=1')
plt.legend()
plt.title("Various Normal pdfs")
plt.show()

图 6-2. 各种正态分布的 PDF

当μ = 0 且σ = 1 时，称为标准正态分布。如果Z是一个标准正态随机变量，则有：

X = σ Z + μ

也服从正态分布，但其均值为μ，标准差为σ 。反之，如果X是均值为μ，标准差为σ的正态随机变量，

Z = ( X - μ ) / σ

是一个标准正态变量。

正态分布的累积分布函数不能用“基本”的方式写出，但我们可以使用 Python 的math.erf 误差函数来表示它：

def normal_cdf(x: float, mu: float = 0, sigma: float = 1) -> float:
    return (1 + math.erf((x - mu) / math.sqrt(2) / sigma)) / 2

再次，在图 6-3 中，我们绘制了几个累积分布函数（CDF）：

xs = [x / 10.0 for x in range(-50, 50)]
plt.plot(xs,[normal_cdf(x,sigma=1) for x in xs],'-',label='mu=0,sigma=1')
plt.plot(xs,[normal_cdf(x,sigma=2) for x in xs],'--',label='mu=0,sigma=2')
plt.plot(xs,[normal_cdf(x,sigma=0.5) for x in xs],':',label='mu=0,sigma=0.5')
plt.plot(xs,[normal_cdf(x,mu=-1) for x in xs],'-.',label='mu=-1,sigma=1')
plt.legend(loc=4) # bottom right
plt.title("Various Normal cdfs")
plt.show()

图 6-3. 各种正态分布的累积分布函数

有时我们需要求逆normal_cdf以找到对应于指定概率的值。虽然没有简单的方法来计算其逆，但normal_cdf是连续且严格递增的，因此我们可以使用二分查找：

def inverse_normal_cdf(p: float,
                       mu: float = 0,
                       sigma: float = 1,
                       tolerance: float = 0.00001) -> float:
    """Find approximate inverse using binary search"""

    # if not standard, compute standard and rescale
    if mu != 0 or sigma != 1:
        return mu + sigma * inverse_normal_cdf(p, tolerance=tolerance)

    low_z = -10.0                      # normal_cdf(-10) is (very close to) 0
    hi_z  =  10.0                      # normal_cdf(10)  is (very close to) 1
    while hi_z - low_z > tolerance:
        mid_z = (low_z + hi_z) / 2     # Consider the midpoint
        mid_p = normal_cdf(mid_z)      # and the CDF's value there
        if mid_p < p:
            low_z = mid_z              # Midpoint too low, search above it
        else:
            hi_z = mid_z               # Midpoint too high, search below it

    return mid_z

该函数重复地将区间二分，直到找到接近所需概率的Z。

中心极限定理

正态分布如此有用的一个原因是中心极限定理，它（本质上）表明大量独立同分布的随机变量的平均值本身近似服从正态分布。

特别是，如果x 1 , ... , x n是均值为μ，标准差为σ的随机变量，并且n很大，则

1 n ( x 1 + ... + x n )

大致上符合正态分布，均值为μ，标准差为σ / n 。同样（但通常更有用），

(x 1 +...+x n )-μn σn

大致上符合均值为 0，标准差为 1 的正态分布。

一个易于说明这一点的方法是看二项式随机变量，它具有两个参数n和p。一个二项式(n,p)随机变量简单地是n个独立的伯努利(p)随机变量的和，其中每个随机变量以概率p等于 1，以概率 1 - p等于 0：

def bernoulli_trial(p: float) -> int:
    """Returns 1 with probability p and 0 with probability 1-p"""
    return 1 if random.random() < p else 0

def binomial(n: int, p: float) -> int:
    """Returns the sum of n bernoulli(p) trials"""
    return sum(bernoulli_trial(p) for _ in range(n))

伯努利变量（Bernoulli(p)）的均值是 p，标准差是 p ( 1 - p ) 。 中心极限定理表明，当 n 变大时，二项分布（Binomial(n,p)）变量近似于均值 μ = n p 和标准差 σ = n p ( 1 - p ) 的正态随机变量。 如果我们同时绘制它们，你可以很容易地看到它们的相似之处：

from collections import Counter

def binomial_histogram(p: float, n: int, num_points: int) -> None:
    """Picks points from a Binomial(n, p) and plots their histogram"""
    data = [binomial(n, p) for _ in range(num_points)]

    # use a bar chart to show the actual binomial samples
    histogram = Counter(data)
    plt.bar([x - 0.4 for x in histogram.keys()],
            [v / num_points for v in histogram.values()],
            0.8,
            color='0.75')

    mu = p * n
    sigma = math.sqrt(n * p * (1 - p))

    # use a line chart to show the normal approximation
    xs = range(min(data), max(data) + 1)
    ys = [normal_cdf(i + 0.5, mu, sigma) - normal_cdf(i - 0.5, mu, sigma)
          for i in xs]
    plt.plot(xs,ys)
    plt.title("Binomial Distribution vs. Normal Approximation")
    plt.show()

例如，当你调用 make_hist(0.75, 100, 10000) 时，你会得到图 6-4 中的图表。

图 6-4. 调用 binomial_histogram 的输出

这个近似的道理是，如果你想知道（比如说）在 100 次投掷中，一个公平硬币出现超过 60 次正面的概率，你可以估计为一个正态分布（Normal(50,5)）大于 60 的概率，这比计算二项分布（Binomial(100,0.5)）的累积分布函数要简单。 （尽管在大多数应用中，你可能会使用愿意计算任何你想要的概率的统计软件。）

进一步探索

• scipy.stats 包含大多数流行概率分布的概率密度函数（PDF）和累积分布函数（CDF）。

• 还记得在第 5 章的结尾处我说过，学习一本统计学教材是个好主意吗？学习一本概率论教材也是个好主意。我知道的最好的在线教材是 Introduction to Probability，由 Charles M. Grinstead 和 J. Laurie Snell（美国数学学会）编写。

第七章：假设与推断

被统计数据感动是真正聪明的人的标志。

乔治·伯纳德·肖

我们会用所有这些统计数据和概率理论做什么呢？数据科学的科学部分经常涉及形成和测试关于数据及其生成过程的假设。

统计假设检验

作为数据科学家，我们经常想测试某个假设是否可能成立。对于我们来说，假设是一些断言，比如“这枚硬币是公平的”或“数据科学家更喜欢 Python 而不是 R”或“人们更有可能在看不见关闭按钮的烦人插页广告出现后，不阅读内容直接离开页面”。这些可以转化为关于数据的统计信息。在各种假设下，这些统计信息可以被视为来自已知分布的随机变量的观察结果，这使我们能够对这些假设的可能性进行陈述。

在经典设置中，我们有一个零假设，H 0 ，代表某些默认位置，以及我们想要与之比较的一些备选假设，H 1 。我们使用统计数据来决定我们是否可以拒绝H 0 作为虚假的。这可能通过一个例子更容易理解。

示例：抛硬币

想象我们有一枚硬币，我们想测试它是否公平。我们假设硬币有一定概率 p 出现正面，因此我们的零假设是硬币是公平的——也就是说，p = 0.5。我们将这个假设与备选假设 p ≠ 0.5 进行测试。

特别是，我们的测试将涉及抛硬币 n 次，并计算正面出现的次数 X。每次抛硬币都是伯努利试验，这意味着 X 是一个二项分布变量，我们可以用正态分布（如我们在 第六章 中看到的）来近似它：

from typing import Tuple
import math

def normal_approximation_to_binomial(n: int, p: float) -> Tuple[float, float]:
    """Returns mu and sigma corresponding to a Binomial(n, p)"""
    mu = p * n
    sigma = math.sqrt(p * (1 - p) * n)
    return mu, sigma

每当一个随机变量遊服从正态分布时，我们可以使用 normal_cdf 来确定其实现值在特定区间内或外的概率：

from scratch.probability import normal_cdf

# The normal cdf _is_ the probability the variable is below a threshold
normal_probability_below = normal_cdf

# It's above the threshold if it's not below the threshold
def normal_probability_above(lo: float,
                             mu: float = 0,
                             sigma: float = 1) -> float:
    """The probability that an N(mu, sigma) is greater than lo."""
    return 1 - normal_cdf(lo, mu, sigma)

# It's between if it's less than hi, but not less than lo
def normal_probability_between(lo: float,
                               hi: float,
                               mu: float = 0,
                               sigma: float = 1) -> float:
    """The probability that an N(mu, sigma) is between lo and hi."""
    return normal_cdf(hi, mu, sigma) - normal_cdf(lo, mu, sigma)

# It's outside if it's not between
def normal_probability_outside(lo: float,
                               hi: float,
                               mu: float = 0,
                               sigma: float = 1) -> float:
    """The probability that an N(mu, sigma) is not between lo and hi."""
    return 1 - normal_probability_between(lo, hi, mu, sigma)

我们也可以反过来——找到非尾部区域或（对称的）包含一定概率水平的平均值的区间。例如，如果我们想找到一个以均值为中心且包含 60% 概率的区间，那么我们找到上下尾部各包含 20% 概率的截止点（留下 60%）：

from scratch.probability import inverse_normal_cdf

def normal_upper_bound(probability: float,
                       mu: float = 0,
                       sigma: float = 1) -> float:
    """Returns the z for which P(Z <= z) = probability"""
    return inverse_normal_cdf(probability, mu, sigma)

def normal_lower_bound(probability: float,
                       mu: float = 0,
                       sigma: float = 1) -> float:
    """Returns the z for which P(Z >= z) = probability"""
    return inverse_normal_cdf(1 - probability, mu, sigma)

def normal_two_sided_bounds(probability: float,
                            mu: float = 0,
                            sigma: float = 1) -> Tuple[float, float]:
    """
 Returns the symmetric (about the mean) bounds
 that contain the specified probability
 """
    tail_probability = (1 - probability) / 2

    # upper bound should have tail_probability above it
    upper_bound = normal_lower_bound(tail_probability, mu, sigma)

    # lower bound should have tail_probability below it
    lower_bound = normal_upper_bound(tail_probability, mu, sigma)

    return lower_bound, upper_bound

特别是，假设我们选择抛硬币 n = 1,000 次。如果我们的公平假设成立，X 应该近似正态分布，均值为 500，标准差为 15.8：

mu_0, sigma_0 = normal_approximation_to_binomial(1000, 0.5)

我们需要就显著性做出决策——我们愿意做第一类错误的程度，即拒绝H 0即使它是真的。由于历史记载已经遗失，这种愿意通常设定为 5%或 1%。让我们选择 5%。

考虑到如果X落在以下界限之外则拒绝H 0的测试：

# (469, 531)
lower_bound, upper_bound = normal_two_sided_bounds(0.95, mu_0, sigma_0)

假设p真的等于 0.5（即H 0为真），我们只有 5%的机会观察到一个落在这个区间之外的X，这正是我们想要的显著性水平。换句话说，如果H 0为真，那么大约有 20 次中有 19 次这个测试会给出正确的结果。

我们还经常关注测试的功效，即不犯第二类错误（“假阴性”）的概率，即我们未能拒绝H 0，尽管它是错误的。为了衡量这一点，我们必须明确H 0假设为假的含义。 （仅知道p不等于 0.5 并不能给我们关于X分布的大量信息。）特别是，让我们检查p真的是 0.55 的情况，这样硬币稍微偏向正面。

在这种情况下，我们可以计算测试的功效：

# 95% bounds based on assumption p is 0.5
lo, hi = normal_two_sided_bounds(0.95, mu_0, sigma_0)

# actual mu and sigma based on p = 0.55
mu_1, sigma_1 = normal_approximation_to_binomial(1000, 0.55)

# a type 2 error means we fail to reject the null hypothesis,
# which will happen when X is still in our original interval
type_2_probability = normal_probability_between(lo, hi, mu_1, sigma_1)
power = 1 - type_2_probability      # 0.887

假设相反，我们的零假设是硬币不偏向正面，或者p ≤ 0 . 5。在这种情况下，我们需要一个单边测试，当X远大于 500 时拒绝零假设，但当X小于 500 时不拒绝。因此，一个 5%显著性测试涉及使用normal_probability_below找到概率分布中 95%以下的截止值：

hi = normal_upper_bound(0.95, mu_0, sigma_0)
# is 526 (< 531, since we need more probability in the upper tail)

type_2_probability = normal_probability_below(hi, mu_1, sigma_1)
power = 1 - type_2_probability      # 0.936

这是一个更强大的测试，因为当X低于 469（如果H 1为真的话，这几乎不太可能发生）时，它不再拒绝H 0，而是在X在 526 到 531 之间时拒绝H 0（如果H 1为真的话，这有一定可能发生）。

p-值

对前述测试的另一种思考方式涉及p 值。我们不再基于某个概率截断选择边界，而是计算概率——假设H 0为真——我们会看到一个至少与我们实际观察到的值一样极端的值的概率。

对于我们的双边测试，检验硬币是否公平，我们计算：

def two_sided_p_value(x: float, mu: float = 0, sigma: float = 1) -> float:
    """
 How likely are we to see a value at least as extreme as x (in either
 direction) if our values are from an N(mu, sigma)?
 """
    if x >= mu:
        # x is greater than the mean, so the tail is everything greater than x
        return 2 * normal_probability_above(x, mu, sigma)
    else:
        # x is less than the mean, so the tail is everything less than x
        return 2 * normal_probability_below(x, mu, sigma)

如果我们看到 530 次正面朝上，我们会计算：

two_sided_p_value(529.5, mu_0, sigma_0)   # 0.062

注意

为什么我们使用529.5而不是530？这就是所谓的连续性校正。这反映了normal_probability_between(529.5, 530.5, mu_0, sigma_0)比normal_probability_between(530, 531, mu_0, sigma_0)更好地估计看到 530 枚硬币正面的概率。

相应地，normal_probability_above(529.5, mu_0, sigma_0)更好地估计了至少看到 530 枚硬币正面的概率。你可能已经注意到，我们在生成图 6-4 的代码中也使用了这个。

说服自己这是一个明智的估计的一种方法是通过模拟：

import random

extreme_value_count = 0
for _ in range(1000):
    num_heads = sum(1 if random.random() < 0.5 else 0    # Count # of heads
                    for _ in range(1000))                # in 1000 flips,
    if num_heads >= 530 or num_heads <= 470:             # and count how often
        extreme_value_count += 1                         # the # is 'extreme'

# p-value was 0.062 => ~62 extreme values out of 1000
assert 59 < extreme_value_count < 65, f"{extreme_value_count}"

由于p-值大于我们的 5%显著性水平，我们不拒绝零假设。如果我们看到 532 枚硬币正面，p-值将是：

two_sided_p_value(531.5, mu_0, sigma_0)   # 0.0463

这比 5%的显著性水平还要小，这意味着我们将拒绝零假设。这和之前完全一样的检验，只是统计学上的不同方法而已。

同样，我们会得到：

upper_p_value = normal_probability_above
lower_p_value = normal_probability_below

对于我们的单边检验，如果我们看到 525 枚硬币正面，我们会计算：

upper_p_value(524.5, mu_0, sigma_0) # 0.061

这意味着我们不会拒绝零假设。如果我们看到 527 枚硬币正面，计算将是：

upper_p_value(526.5, mu_0, sigma_0) # 0.047

我们将拒绝零假设。

警告

在使用normal_probability_above计算p-值之前，请确保你的数据大致服从正态分布。糟糕的数据科学充满了人们声称某些观察到的事件在随机发生的机会是百万分之一的例子，当他们真正的意思是“机会，假设数据是正态分布的”，如果数据不是的话，这是相当毫无意义的。

有各种检验正态性的统计方法，但即使是绘制数据也是一个很好的开始。

置信区间

我们一直在测试关于硬币正面概率p值的假设，这是未知“正面”分布的一个参数。在这种情况下，第三种方法是围绕参数的观察值构建一个置信区间。

例如，我们可以通过观察与每次翻转相对应的伯努利变量的平均值来估计不公平硬币的概率—如果是正面则为 1，如果是反面则为 0。如果我们在 1,000 次翻转中观察到 525 枚硬币正面，则我们估计p等于 0.525。

我们对这个估计有多么有信心？嗯，如果我们知道p的确切值，那么中心极限定理（回顾“中心极限定理”）告诉我们，这些伯努利变量的平均值应该近似正态分布，均值为p，标准差为：

math.sqrt(p * (1 - p) / 1000)

在这里我们不知道p的值，所以我们使用我们的估计值：

p_hat = 525 / 1000
mu = p_hat
sigma = math.sqrt(p_hat * (1 - p_hat) / 1000)   # 0.0158

这并不完全合理，但人们似乎仍然这样做。使用正态近似，我们得出“95%的置信度”下以下区间包含真实参数p的结论：

normal_two_sided_bounds(0.95, mu, sigma)        # [0.4940, 0.5560]

注意

这是关于区间而不是关于p的说法。你应该理解它是这样的断言：如果你重复进行实验很多次，那么 95%的时间，“真实”的参数（每次都相同）会落在观察到的置信区间内（每次可能不同）。

特别是，我们不会得出硬币不公平的结论，因为 0.5 位于我们的置信区间内。

如果我们看到了 540 次正面，那么我们将会：

p_hat = 540 / 1000
mu = p_hat
sigma = math.sqrt(p_hat * (1 - p_hat) / 1000) # 0.0158
normal_two_sided_bounds(0.95, mu, sigma) # [0.5091, 0.5709]

在这里，“公平硬币”的置信区间中不存在。（如果真的是公平硬币假设，它不会通过一个测试。）

p-操纵

一个过程，只有 5%的时间错误地拒绝原假设，按定义来说：

from typing import List

def run_experiment() -> List[bool]:
    """Flips a fair coin 1000 times, True = heads, False = tails"""
    return [random.random() < 0.5 for _ in range(1000)]

def reject_fairness(experiment: List[bool]) -> bool:
    """Using the 5% significance levels"""
    num_heads = len([flip for flip in experiment if flip])
    return num_heads < 469 or num_heads > 531

random.seed(0)
experiments = [run_experiment() for _ in range(1000)]
num_rejections = len([experiment
                      for experiment in experiments
                      if reject_fairness(experiment)])

assert num_rejections == 46

这意味着，如果你试图找到“显著”结果，通常你可以找到。对数据集测试足够多的假设，几乎肯定会出现一个显著结果。移除正确的异常值，你可能会将p-值降低到 0.05 以下。（我们在“相关性”中做了类似的事情；你注意到了吗？）

这有时被称为p-操纵，在某种程度上是“从p-值框架推断”的结果。一篇批评这种方法的好文章是“地球是圆的”，作者是雅各布·科恩。

如果你想进行良好的科学研究，你应该在查看数据之前确定你的假设，清理数据时不要考虑假设，并且要记住p-值并不能替代常识。（另一种方法在“贝叶斯推断”中讨论。）

例如：运行 A/B 测试

你在 DataSciencester 的主要职责之一是体验优化，这是一个试图让人们点击广告的委婉说法。你的一个广告客户开发了一种新的面向数据科学家的能量饮料，广告部副总裁希望你帮助选择广告 A（“味道棒！”）和广告 B（“更少偏见！”）之间的区别。

作为科学家，你决定通过随机向站点访问者展示两个广告，并跟踪点击每个广告的人数来运行实验。

如果在 1,000 个 A 观众中有 990 个点击他们的广告，而在 1,000 个 B 观众中只有 10 个点击他们的广告，你可以相当有信心认为 A 是更好的广告。但如果差异不那么明显呢？这就是你会使用统计推断的地方。

假设N A人看到广告 A，并且其中n A人点击了它。我们可以将每次广告浏览视为伯努利试验，其中p A是某人点击广告 A 的概率。那么（如果N A很大，这里是这样），我们知道n A / N A大致上是一个均值为p A，标准差为σ A = p A ( 1 - p A ) / N A的正态随机变量。

同样，n B / N B大致上是一个均值为p B，标准差为σ B = p B ( 1 - p B ) / N B。我们可以用代码表示这个过程：

def estimated_parameters(N: int, n: int) -> Tuple[float, float]:
    p = n / N
    sigma = math.sqrt(p * (1 - p) / N)
    return p, sigma

如果我们假设这两个正态分布是独立的（这似乎是合理的，因为单个伯努利试验应该是独立的），那么它们的差异也应该是均值为p B - p A，标准差为σ A 2 + σ B 2。

注意

这有点作弊。只有当您知道标准偏差时，数学才能完全工作。在这里，我们是从数据中估计它们，这意味着我们确实应该使用t分布。但是对于足够大的数据集，它接近标准正态分布，所以差别不大。

这意味着我们可以测试零假设，即p A和p B相同（即p A - p B为 0），使用的统计量是：

def a_b_test_statistic(N_A: int, n_A: int, N_B: int, n_B: int) -> float:
    p_A, sigma_A = estimated_parameters(N_A, n_A)
    p_B, sigma_B = estimated_parameters(N_B, n_B)
    return (p_B - p_A) / math.sqrt(sigma_A ** 2 + sigma_B ** 2)

应该大约是一个标准正态分布。

例如，如果“味道好”的广告在 1,000 次浏览中获得了 200 次点击，“更少偏见”的广告在 1,000 次浏览中获得了 180 次点击，则统计量等于：

z = a_b_test_statistic(1000, 200, 1000, 180)    # -1.14

如果实际上平均值相等，观察到这么大差异的概率将是：

two_sided_p_value(z)                            # 0.254

这个概率足够大，我们无法得出有太大差异的结论。另一方面，如果“更少偏见”的点击仅为 150 次，则有：

z = a_b_test_statistic(1000, 200, 1000, 150)    # -2.94
two_sided_p_value(z)                            # 0.003

这意味着如果广告效果相同，我们看到这么大的差异的概率只有 0.003。

贝叶斯推断

我们所看到的程序涉及对我们的测试做出概率陈述：例如，“如果我们的零假设成立，你观察到如此极端的统计量的概率只有 3%”。

推理的另一种方法涉及将未知参数本身视为随机变量。分析员（也就是你）从参数的先验分布开始，然后使用观察到的数据和贝叶斯定理来获得参数的更新后验分布。与对测试进行概率判断不同，你对参数进行概率判断。

例如，当未知参数是概率时（如我们抛硬币的示例），我们通常使用Beta 分布中的先验，该分布将其所有概率都放在 0 和 1 之间：

def B(alpha: float, beta: float) -> float:
    """A normalizing constant so that the total probability is 1"""
    return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta)

def beta_pdf(x: float, alpha: float, beta: float) -> float:
    if x <= 0 or x >= 1:          # no weight outside of [0, 1]
        return 0
    return x ** (alpha - 1) * (1 - x) ** (beta - 1) / B(alpha, beta)

一般来说，这个分布将其权重集中在：

alpha / (alpha + beta)

而且alpha和beta越大，分布就越“紧密”。

例如，如果alpha和beta都为 1，那就是均匀分布（以 0.5 为中心，非常分散）。如果alpha远大于beta，大部分权重都集中在 1 附近。如果alpha远小于beta，大部分权重都集中在 0 附近。图 7-1 展示了几种不同的 Beta 分布。

图 7-1。示例 Beta 分布

假设我们对p有一个先验分布。也许我们不想对硬币是否公平发表立场，我们选择alpha和beta都等于 1。或者我们非常相信硬币 55%的时间会正面朝上，我们选择alpha等于 55，beta等于 45。

然后我们多次抛硬币，看到h次正面和t次反面。贝叶斯定理（以及一些在这里过于繁琐的数学）告诉我们p的后验分布再次是 Beta 分布，但参数为alpha + h和beta + t。

注意

后验分布再次是 Beta 分布并非巧合。头数由二项式分布给出，而 Beta 是与二项式分布共轭先验。这意味着每当您使用相应的二项式观察更新 Beta 先验时，您将得到一个 Beta 后验。

假设你抛了 10 次硬币，只看到 3 次正面。如果你从均匀先验开始（在某种意义上拒绝对硬币的公平性发表立场），你的后验分布将是一个 Beta(4, 8)，中心在 0.33 左右。由于你认为所有概率都是同等可能的，你的最佳猜测接近观察到的概率。

如果你最初使用一个 Beta(20, 20)（表达了一种认为硬币大致公平的信念），你的后验分布将是一个 Beta(23, 27)，中心在 0.46 左右，表明修正后的信念可能硬币稍微偏向反面。

如果你最初假设一个 Beta(30, 10)（表达了一种认为硬币有 75%概率翻转为正面的信念），你的后验分布将是一个 Beta(33, 17)，中心在 0.66 左右。在这种情况下，你仍然相信硬币有正面的倾向，但不像最初那样强烈。这三种不同的后验分布在图 7-2 中绘制出来。

图 7-2。来自不同先验的后验分布

如果你不断地翻转硬币，先验的影响将越来越小，最终你将拥有（几乎）相同的后验分布，无论你最初选择了哪个先验。

例如，无论你最初认为硬币有多么偏向正面，看到 2000 次翻转中有 1000 次正面后，很难维持那种信念（除非你选择像 Beta(1000000,1)这样的先验，这是疯子的行为）。

有趣的是，这使我们能够对假设进行概率性陈述：“基于先验和观察数据，硬币的正面概率在 49%到 51%之间的可能性仅为 5%。”这在哲学上与“如果硬币是公平的，我们预期观察到极端数据的概率也仅为 5%”这样的陈述有很大的不同。

使用贝叶斯推断来测试假设在某种程度上被认为是有争议的——部分原因是因为数学可以变得相当复杂，部分原因是因为选择先验的主观性质。我们在本书中不会进一步使用它，但了解这一点是很好的。

进一步探索

• 我们只是浅尝辄止地探讨了统计推断中你应该了解的内容。第五章末推荐的书籍详细讨论了这些内容。

• Coursera 提供了一门数据分析与统计推断课程，涵盖了许多这些主题。

第八章：梯度下降

那些夸口自己的衰退的人，吹嘘他们欠别人的东西。

塞内加

在做数据科学时，我们经常会试图找到某种情况下的最佳模型。通常，“最佳”意味着“最小化预测误差”或“最大化数据的可能性”。换句话说，它将代表某种优化问题的解决方案。

这意味着我们需要解决许多优化问题。特别是，我们需要从头开始解决它们。我们的方法将是一种称为梯度下降的技术，它非常适合从头开始的处理。你可能不会觉得它本身特别令人兴奋，但它将使我们能够在整本书中做一些令人兴奋的事情，所以请忍耐一下。

梯度下降背后的思想

假设我们有一个函数f，它以一个实数向量作为输入并输出一个实数。一个简单的这样的函数是：

from scratch.linear_algebra import Vector, dot

def sum_of_squares(v: Vector) -> float:
    """Computes the sum of squared elements in v"""
    return dot(v, v)

我们经常需要最大化或最小化这样的函数。也就是说，我们需要找到产生最大（或最小）可能值的输入v。

对于像我们这样的函数，梯度（如果你记得你的微积分，这是偏导数的向量）给出了函数增加最快的输入方向。（如果你不记得你的微积分，相信我或者去互联网上查一查。）

相应地，最大化一个函数的一种方法是选择一个随机起点，计算梯度，沿着梯度的方向迈出一小步（即导致函数增加最多的方向），然后用新的起点重复这个过程。类似地，你可以尝试通过向相反方向迈小步来最小化一个函数，如图 8-1 所示。

图 8-1. 使用梯度下降找到最小值

注意

如果一个函数有一个唯一的全局最小值，这个过程很可能会找到它。如果一个函数有多个（局部）最小值，这个过程可能会“找到”它们中的错误之一，这种情况下你可能需要从不同的起点重新运行该过程。如果一个函数没有最小值，那么可能这个过程会永远进行下去。

估计梯度

如果f是一个关于一个变量的函数，在点x处的导数测量了当我们对x做一个非常小的变化时f(x)如何改变。导数被定义为差商的极限：

from typing import Callable

def difference_quotient(f: Callable[[float], float],
                        x: float,
                        h: float) -> float:
    return (f(x + h) - f(x)) / h

当h趋近于零时。

（许多想要学习微积分的人都被极限的数学定义所困扰，这是美丽的但可能看起来有些令人生畏的。在这里，我们将作弊并简单地说，“极限”意味着你认为的意思。）

导数是切线在 ( x , f ( x ) ) 处的斜率，而差商是穿过 ( x + h , f ( x + h ) ) 的不完全切线的斜率。随着 h 变得越来越小，这个不完全切线越来越接近切线 (图 8-2)。

图 8-2. 使用差商近似导数

对于许多函数来说，精确计算导数是很容易的。例如，square 函数：

def square(x: float) -> float:
    return x * x

具有导数：

def derivative(x: float) -> float:
    return 2 * x

对于我们来说很容易检查的是，通过明确计算差商并取极限。(这只需要高中代数就可以做到。)

如果你不能（或不想）找到梯度怎么办？尽管我们不能在 Python 中取极限，但我们可以通过评估一个非常小的 e 的差商来估计导数。图 8-3 显示了这样一个估计的结果：

xs = range(-10, 11)
actuals = [derivative(x) for x in xs]
estimates = [difference_quotient(square, x, h=0.001) for x in xs]

# plot to show they're basically the same
import matplotlib.pyplot as plt
plt.title("Actual Derivatives vs. Estimates")
plt.plot(xs, actuals, 'rx', label='Actual')       # red  x
plt.plot(xs, estimates, 'b+', label='Estimate')   # blue +
plt.legend(loc=9)
plt.show()

图 8-3. 差商近似的好处

当 f 是多变量函数时，它具有多个 偏导数，每个偏导数指示当我们只对其中一个输入变量进行微小变化时 f 如何变化。

我们通过将其视为仅其第 i 个变量的函数，并固定其他变量来计算其 ith 偏导数：

def partial_difference_quotient(f: Callable[[Vector], float],
                                v: Vector,
                                i: int,
                                h: float) -> float:
    """Returns the i-th partial difference quotient of f at v"""
    w = [v_j + (h if j == i else 0)    # add h to just the ith element of v
         for j, v_j in enumerate(v)]

    return (f(w) - f(v)) / h

然后我们可以以同样的方式估计梯度：

def estimate_gradient(f: Callable[[Vector], float],
                      v: Vector,
                      h: float = 0.0001):
    return [partial_difference_quotient(f, v, i, h)
            for i in range(len(v))]

注意

这种“使用差商估计”的方法的一个主要缺点是计算成本高昂。如果 v 的长度为 n，estimate_gradient 必须在 2n 个不同的输入上评估 f。如果你反复估计梯度，那么你需要做大量额外的工作。在我们所做的一切中，我们将使用数学来显式计算我们的梯度函数。

使用梯度

很容易看出，当其输入 v 是一个零向量时，sum_of_squares 函数最小。但想象一下，如果我们不知道这一点。让我们使用梯度来找到所有三维向量中的最小值。我们只需选择一个随机起点，然后沿着梯度的相反方向迈出微小的步骤，直到我们到达梯度非常小的点：

import random
from scratch.linear_algebra import distance, add, scalar_multiply

def gradient_step(v: Vector, gradient: Vector, step_size: float) -> Vector:
    """Moves `step_size` in the `gradient` direction from `v`"""
    assert len(v) == len(gradient)
    step = scalar_multiply(step_size, gradient)
    return add(v, step)

def sum_of_squares_gradient(v: Vector) -> Vector:
    return [2 * v_i for v_i in v]

# pick a random starting point
v = [random.uniform(-10, 10) for i in range(3)]

for epoch in range(1000):
    grad = sum_of_squares_gradient(v)    # compute the gradient at v
    v = gradient_step(v, grad, -0.01)    # take a negative gradient step
    print(epoch, v)

assert distance(v, [0, 0, 0]) < 0.001    # v should be close to 0

如果你运行这个程序，你会发现它总是最终得到一个非常接近 [0,0,0] 的 v。你运行的周期越多，它就会越接近。

选择合适的步长

尽管反对梯度的理由很明确，但移动多远并不清楚。事实上，选择正确的步长更多地是一门艺术而不是一门科学。流行的选择包括：

• 使用固定步长

• 随着时间逐渐缩小步长

• 在每一步中，选择使目标函数值最小化的步长。

最后一种方法听起来很好，但在实践中是一种昂贵的计算。为了保持简单，我们将主要使用固定步长。“有效”的步长取决于问题——太小，你的梯度下降将永远不会结束；太大，你将采取巨大的步骤，可能会使你关心的函数变得更大甚至未定义。因此，我们需要进行实验。

使用梯度下降拟合模型

在这本书中，我们将使用梯度下降来拟合数据的参数化模型。通常情况下，我们会有一些数据集和一些（假设的）依赖于一个或多个参数的数据模型（以可微分的方式）。我们还会有一个损失函数，用于衡量模型拟合数据的好坏程度（越小越好）。

如果我们将数据视为固定的，那么我们的损失函数告诉我们任何特定模型参数的好坏程度。这意味着我们可以使用梯度下降找到使损失尽可能小的模型参数。让我们看一个简单的例子：

# x ranges from -50 to 49, y is always 20 * x + 5
inputs = [(x, 20 * x + 5) for x in range(-50, 50)]

在这种情况下，我们知道线性关系的参数x和y，但想象一下我们希望从数据中学习这些参数。我们将使用梯度下降来找到最小化平均平方误差的斜率和截距。

我们将以一个基于单个数据点误差的梯度确定梯度的函数开始：

def linear_gradient(x: float, y: float, theta: Vector) -> Vector:
    slope, intercept = theta
    predicted = slope * x + intercept    # The prediction of the model.
    error = (predicted - y)              # error is (predicted - actual).
    squared_error = error ** 2           # We'll minimize squared error
    grad = [2 * error * x, 2 * error]    # using its gradient.
    return grad

让我们考虑一下梯度的含义。假设对于某个x，我们的预测值过大。在这种情况下，error是正的。第二个梯度项，2 * error，是正的，这反映了截距的小增加会使（已经过大的）预测值变得更大，进而导致该x的平方误差变得更大。

第一个梯度项2 * error * x与x的符号相同。确实，如果x为正，那么斜率的小增加会再次使预测（因此误差）变大。但是如果x为负，斜率的小增加会使预测（因此误差）变小。

现在，这个计算是针对单个数据点的。对于整个数据集，我们将看均方误差。均方误差的梯度就是单个梯度的均值。

所以，我们要做的是：

1. 从一个随机值theta开始。

2. 计算梯度的均值。

3. 在那个方向上调整theta。

4. 重复。

经过许多epochs（我们称每次通过数据集的迭代），我们应该学到一些正确的参数：

from scratch.linear_algebra import vector_mean

# Start with random values for slope and intercept
theta = [random.uniform(-1, 1), random.uniform(-1, 1)]

learning_rate = 0.001

for epoch in range(5000):
    # Compute the mean of the gradients
    grad = vector_mean([linear_gradient(x, y, theta) for x, y in inputs])
    # Take a step in that direction
    theta = gradient_step(theta, grad, -learning_rate)
    print(epoch, theta)

slope, intercept = theta
assert 19.9 < slope < 20.1,   "slope should be about 20"
assert 4.9 < intercept < 5.1, "intercept should be about 5"

小批量和随机梯度下降

前述方法的一个缺点是，我们必须在可以采取梯度步骤并更新参数之前对整个数据集进行梯度评估。在这种情况下，这是可以接受的，因为我们的数据集只有 100 对，梯度计算是廉价的。

然而，您的模型通常会有大型数据集和昂贵的梯度计算。在这种情况下，您会希望更频繁地进行梯度步骤。

我们可以使用一种称为小批量梯度下降的技术来实现这一点，通过从更大的数据集中抽取“小批量”来计算梯度（并执行梯度步骤）：

from typing import TypeVar, List, Iterator

T = TypeVar('T')  # this allows us to type "generic" functions

def minibatches(dataset: List[T],
                batch_size: int,
                shuffle: bool = True) -> Iterator[List[T]]:
    """Generates `batch_size`-sized minibatches from the dataset"""
    # start indexes 0, batch_size, 2 * batch_size, ...
    batch_starts = [start for start in range(0, len(dataset), batch_size)]

    if shuffle: random.shuffle(batch_starts)  # shuffle the batches

    for start in batch_starts:
        end = start + batch_size
        yield dataset[start:end]

注意

TypeVar(T)允许我们创建一个“泛型”函数。它表示我们的dataset可以是任何一种类型的列表——str、int、list等等——但无论是哪种类型，输出都将是它们的批次。

现在我们可以再次使用小批量来解决我们的问题：

theta = [random.uniform(-1, 1), random.uniform(-1, 1)]

for epoch in range(1000):
    for batch in minibatches(inputs, batch_size=20):
        grad = vector_mean([linear_gradient(x, y, theta) for x, y in batch])
        theta = gradient_step(theta, grad, -learning_rate)
    print(epoch, theta)

slope, intercept = theta
assert 19.9 < slope < 20.1,   "slope should be about 20"
assert 4.9 < intercept < 5.1, "intercept should be about 5"

另一种变体是随机梯度下降，在这种方法中，您基于一个训练样本一次进行梯度步骤：

theta = [random.uniform(-1, 1), random.uniform(-1, 1)]

for epoch in range(100):
    for x, y in inputs:
        grad = linear_gradient(x, y, theta)
        theta = gradient_step(theta, grad, -learning_rate)
    print(epoch, theta)

slope, intercept = theta
assert 19.9 < slope < 20.1,   "slope should be about 20"
assert 4.9 < intercept < 5.1, "intercept should be about 5"

在这个问题上，随机梯度下降在更少的迭代次数内找到了最优参数。但总是存在权衡。基于小型小批量（或单个数据点）进行梯度步骤可以让您执行更多次梯度步骤，但单个点的梯度可能与整个数据集的梯度方向差异很大。

此外，如果我们不是从头开始进行线性代数运算，通过“向量化”我们在批次中的计算而不是逐点计算梯度，会有性能提升。

在整本书中，我们将尝试找到最佳的批量大小和步长。

注意

关于各种梯度下降方法的术语并不统一。通常称为批量梯度下降的方法是“对整个数据集计算梯度”，有些人在提到小批量版本时称为随机梯度下降（其中逐点版本是一种特例）。

进一步探索

• 继续阅读！我们将在本书的其余部分使用梯度下降来解决问题。

• 此时，您无疑对我建议您阅读教科书感到厌倦。如果能稍微安慰一下的话，Active Calculus 1.0，由 Matthew Boelkins、David Austin 和 Steven Schlicker（Grand Valley State University Libraries）编写的书，似乎比我学习的微积分教科书更好。

• Sebastian Ruder 在他的 epic 博客文章 中比较了梯度下降及其许多变体。