Python-统计和微积分研讨会（三）

Python 统计和微积分研讨会（三）

原文：zh.annas-archive.org/md5/6cbaed7d834977b8ea96cc7aa6d8a083

译者：飞龙

协议：CC BY-NC-SA 4.0

第五章：用 Python 进行更多数学

概述

在本章结束时，你将能够掌握序列和级数的基本概念，并编写实现这些概念的 Python 函数。你将了解基本三角函数及其应用之间的关系，比如著名的毕达哥拉斯定理。你将练习向量微积分，并通过在 Python 中进行向量代数来了解它的适用性。最后，你会感到高兴，因为复数并不比其他类型的数字差；它们与三角学密切相关，并且对现实世界的应用很有用。

介绍

在上一章中，我们用 Python 介绍了函数和代数，从基本函数开始，然后进行变换和解方程。在本章中，我们将介绍序列和级数，在现实世界中有许多应用，比如金融，也是理解微积分的基础。此外，我们还将探讨三角学、向量和复数，以便更好地理解数学世界。

任何出色的 Python 程序员的核心技能包括对背景数学的深刻理解以及有效地应用它们。把向量和复数看作是我们数学工具箱的有价值扩展，它们将有助于有效地描述、量化和解决来自金融、科学、商业和社会领域的现实问题。

序列和级数等在利润、损失、股利或其他支付定期发生的情况中出现。三角学和三角函数对解决地理空间问题至关重要，而向量在物理学和工程学、机器学习等领域得到广泛应用，其中多个不同的值被分组在一起，方向的概念至关重要。复数是一些最基本的概念，在电磁学、光学、量子力学和计算机科学等领域有广泛的应用。

序列和级数

如果你参加一个电视节目，其中 1 万美元的问题是“给定数字 2、4、8、16 和 32，序列中接下来是什么？”，你最好的猜测是什么？如果你的回答是 64，那么恭喜你——你刚刚更接近理解数学抽象中的一个关键概念：序列。序列就像普通意义上的顺序一样，是一种事物相互跟随的特定顺序。在这里，事物（在大多数情况下）是相关的整数或实数。元素的顺序很重要。元素也被称为序列的成员或术语。

例如，在你参与的电视节目的前述序列中，每个术语都是前一个数字乘以 2 得到的；这个序列没有结束，因为可以无限产生整数数字。在其他情况下，序列中的元素可以出现多次。想想一年中每个月的天数，或者随机事件的结果序列，比如抛硬币的结果。自古印度以来就已知的一个著名序列是斐波那契数列——1、1、2、3、5、8、13……这个序列中，每个新术语都是前两个术语的和。

也就是说，我们需要知道至少两个术语才能推导出其他任何术语。换句话说，我们需要读取前两个数字（在前述序列中是 1 和 1，但通常是任意两个数字）才能推导和预测第三个数字。我们知道一些序列，比如斐波那契数列，其中包含一些内在的逻辑；我们可以遵循一个基本规则，并推导出序列的任何术语。

在本章中，我们将重点关注基本序列，也称为级数，这些序列在应用数学和编程的许多领域中反复出现，属于三种基本类别之一：等差、等比和递归。这些不是唯一的可能性；然而，它们是最受欢迎的序列家族，并且说明了它们所蕴含的逻辑。

一系列数字{αn} = {α1, α2, α3, ..., αΝ, ...}是一个有序的术语（元素或成员）的集合，对于这个集合有一个规则，将每个自然数 n = 1, 2, 3, ..., N 与序列中的一个术语关联起来。序列的长度（即其术语的数量）可以是有限的或无限的，因此相应地称为有限或无限序列。

一个级数是一个按照以下方式求和的数学序列：

图 5.1：级数方程

这也可以使用求和符号进行求和，如下所示：

图 5.2：无限级数的方程

在前面的情况下，我们的级数是无限的（即，它是无限序列的所有项的和）。然而，一个级数，就像一个序列一样，也可以是有限的。为什么一个总和会有无限个项？因为事实证明，在许多情况下，通过应用已知的公式进行计算更有效。此外，即使序列是无限的，求和也可以收敛到一个数字（不是无穷大）或某个函数。由于这个原因，级数可以被认为是已知函数的构建块，它们的术语可以用来表示越来越复杂的函数，从而使得对它们的性质的研究变得直观。级数和序列在数学、物理、工程、金融等领域中无处不在，并且自古以来就为人所知。它们出现并且作为无限和的定义在导数和其他函数中特别有用。

等差数列

像大多数数学概念一样，序列无处不在我们的日常生活中。你可能以前没有想过，但每次你乘坐出租车时，后台都会运行一个序列来计算你的乘车总费用。有一个初始费用，每公里（或英里）增加一个固定金额。因此，在任何给定时刻，都有一个实际的对应数字（到目前为止的乘车费用）。所有这些小计的有序集合形成一个序列。同样，随着你的成长，你的身高是一个随时间（天或月）变化的实数序列（以厘米或英寸表示的身高）。这两个例子都构成了随时间非递减的序列——换句话说，每个术语要么大于等于任何先前的术语，但从不小于。然而，这两个例子之间存在微妙的差异：随着我们的成长，我们增加身高的速度不同（儿童增长速度快，青少年增长速度慢，成年人增长速度为零），而出租车费用增加的速度是恒定的。这导致我们需要引入一种特殊类别的序列——等差数列，其定义如下。

任意两个连续项之间的差是恒定的序列称为等差。因此，等差数列的公式如下：αn+1- αn = d

在这里，d是常数，必须对所有n成立。当然，很明显，如果你知道参数d和一些（任何）项αn，那么项αn+1 可以通过前面的关系直接找到。通过重复，所有项αn+2，αn+3 ...，以及项αn-1，αn-2 都可以找到。换句话说，如果你知道参数d和序列的第一项α1，那么所有项（即唯一确定的）的序列都是已知的。给出我们序列的第n*项的一般公式如下：

αn = α1 + (n – 1)d

在这里，d被称为公差。

相反，要测试一个通用序列是否是等差数列，只需检查其项的所有成对差异αn+1 – αn，并查看这些是否是相同的常数。在相应的等差数列中，前面序列的和变成了下面的形式：

Σnj αj = Σnj [ α1 + (j – 1)d ] = n(α1 + αn)/2

这意味着通过知道长度n，序列的第一个和最后一个项，我们可以确定从α1 到αn 的所有项的和。请注意，和（α1 + αn）给出整个序列的算术平均数的两倍，因此该序列不过是算术平均数的n倍。

现在，我们知道算术序列的主要逻辑和组成部分。现在，让我们看一些具体的例子。目前，我们不需要在 Python 中导入任何特定的库，因为我们将创建自己的函数。让我们提醒自己，这些函数总是需要以def开头，后面跟着一个空格，函数名（任何我们喜欢的东西），以及函数在括号内接受的参数列表，后面跟着一个分号。接下来的行是缩进的（向右缩进四个空格），并且是函数的逻辑，也就是函数的算法或方法。例如，考虑以下例子：

def my_function(arg1, arg2):
    '''Write a function that adds two numbers 
       and returns their sum'''
    result = arg1 + arg2
    return result

在最终语句result之后，返回的是函数的返回值。因此，例如，如果我们正在编写前面的my_function定义，该定义接收两个输入数字arg1和arg2，那么我们可以将它传递给一个新变量，比如以下变量：

summed = my_function(2,9)
print(summed)

输出将如下所示：

11

在这里，summed是一个新变量，它正是由my_function返回（生成）的。请注意，如果在函数定义中缺少return语句，则语法仍然正确，函数仍然可以被调用。但是，summed变量将等于None。

现在，如果我们想创建一个（任何）数字序列，我们应该在我们的函数内包含一个迭代。在 Python 中，可以通过for或while循环来实现这一点。让我们看一个例子，一个函数输出一个n个和的序列：

def my_sequence(arg1, arg2, n):
    '''Write a function that adds two numbers n times and 
       prints their sum'''
    result = 0
    for i in range(n):
        result = result + arg1 + arg2
        print(result)

在这里，我们初始化变量 result（为零），然后迭代地将arg1 + arg2加到它上面。这个迭代发生了n次，其中n也是我们新函数my_sequence的一个参数。每次循环（跟在for语句后面的）执行时，result增加了arg1 + arg2，然后打印在屏幕上。为简单起见，我们在这里省略了return语句。在这里，我们使用了 Python 内置的range()方法，它生成一个整数序列，从 0 开始，以给定的停止整数之前的一个数结束（我们提供的数字）。让我们调用我们的函数：

my_sequence(2,9,4)

我们将获得以下输出：

11
22
33
44

如果我们使用了while循环，我们将得到相同的结果：

def my_sequence(arg1, arg2, n):
    '''Write a function that adds two numbers n times 
       and prints their sum'''
    i = 0
    result = 0
    while i < n:
        result = result + arg1 + arg2
        i += 1
        print(result)

如果我们调用my_sequence函数，我们将得到与先前相同输入的相同输出。

生成器

在 Python 中进行顺序操作的另一个有趣选项是使用生成器。生成器是类似于函数的对象，返回一组可迭代的项目，一次一个值。简单地说，如果一个函数包含至少一个yield语句，它就成为一个生成器函数。使用生成器而不是函数的好处是，我们可以根据需要（这里是无限次数）调用生成器，而不会使系统的内存过载。在某些情况下，它们可以是无价的工具。要获取一系列项的一个项，我们使用next()方法。首先，让我们定义我们的函数：

def my_generator(arg1, arg2, n):
    '''Write a generator function that adds 
       two numbers n times and prints their sum'''
    i = 0
    result = 0
    while i < n:
        result = result + arg1 + arg2
        i += 1
        yield result

现在，让我们多次调用next()方法：

my_gen = my_generator(2,9,4)
next(my_gen)

以下是输出：

11

第二次调用该方法：

next(my_gen)

以下是输出：

22

第三次调用它：

next(my_gen)

以下是输出：

33

第四次调用该方法：

next(my_gen)

以下是输出：

44

因此，我们得到了与上一个示例相同的结果，但是一次一个。如果我们重复调用next()方法，我们将收到错误消息，因为我们已经耗尽了我们的生成器：

next(my_gen)
Traceback (most recent call last):
    File "<stdin>", line 1, in <module>
StopIteration

现在，我们准备在 Python 代码中实现我们学到的序列关系。

练习 5.01：确定算术序列和算术级数的第 n 项

在这个练习中，我们将使用一个简单的 Python 函数创建有限和无限的算术序列。作为输入，我们希望提供序列的第一项a1，公差d和序列的长度n。我们的目标是获得以下内容：

• 序列的一个项（第n项）。

• 完整的数字序列。

• 算术序列的n项的和，以便将其与先前给出的算术级数的结果进行比较。

为了计算前面的目标，我们需要提供序列的第一项a1，公差d和序列的长度n作为输入。让我们实现这个练习：

1. 首先，我们想要编写一个函数，根据通用公式αn = α1 + (n – 1)d 返回第n项：

def a_n(a1, d, n):
    '''Return the n-th term of the arithmetic sequence.
    :a1: first term of the sequence. Integer or real.
    :n: the n-th term in sequence
    returns: n-th term. Integer or real.'''
    an = a1 + (n - 1)*d
    return an

通过这样做，我们获得了序列的第n项，而无需知道任何其他前面的项。例如，让我们用参数(4, 3, 10)调用我们的函数：

a_n(4, 3, 10)

我们将得到以下输出：

31

1. 现在，让我们编写一个函数，它将初始项a1递增d，n次，并将所有项存储在列表中：

def a_seq(a1, d, n):
    '''Obtain the whole arithmetic sequence up to n.
    :a1: first term of the sequence. Integer or real.
    :d: common difference of the sequence. Integer or real.
    :n: length of sequence
    returns: sequence as a list.'''
    sequence = []
    for _ in range(n):
        sequence.append(a1)
        a1 = a1 + d
    return sequence

1. 要检查结果列表，请添加以下代码：

a_seq(4, 3, 10)

输出如下：

[4, 7, 10, 13, 16, 19, 22, 25, 28, 31]

在这里，我们得到了一个长度为 10 的算术序列，从4开始，增加 3。

1. 现在，让我们生成无限序列。我们可以使用之前介绍的 Python 生成器来实现这一点：

def infinite_a_sequence(a1, d):
    while True:
        yield a1
        a1 = a1 + d
for i in infinite_a_sequence(4,3):
    print(i, end=" ")

如果运行上述代码，您会注意到我们必须手动中止执行；否则，for循环将永远打印出序列的元素。使用 Python 生成器的另一种方法是，如前所述，直接在生成器对象（这里是infinite_a_sequence()）上调用next()方法。

1. 让我们通过调用sum() Python 方法来计算我们序列的项的和：

sum(a_seq(4, 3, 10))

输出如下：

175

1. 最后，实现αn = α1 + (n – 1)d 公式，这给我们算术级数，以便我们可以将其与我们的求和结果进行比较：

def a_series(a1, d, n):
    result = n * (a1 + a_n(a1, d, n)) / 2
    return result

1. 运行函数，如下：

a_series(4, 3, 10)

输出如下：

175.0

注意

要访问此特定部分的源代码，请参阅packt.live/2D2S52c.

您还可以在packt.live/31DjRfO.上在线运行此示例

有了这个，我们通过使用序列或级数得到了相同的结果，得到了算术序列元素的总和。能够用两种独立的数学方法交叉验证给定结果的能力对于所有级别的程序员来说都是非常有用的，并且是科学验证的核心。此外，了解不同的方法（在这里，我们用来得到级数结果的两种方法），以及每种方法的优势（以及劣势）对于在高级水平上编写代码至关重要。

我们将学习另一种不同但同样基本的序列类别：几何序列。

几何序列

传染病从一个人传播到另一个或更多人，取决于给定社区的人口密度。在像大流行这样的情况下，对于一种中等传染性的疾病，平均每个患病者每天感染两个人是现实的。因此，如果在第 1 天只有一个感染者，第 2 天将有两个新感染者，第 3 天，对于每个先前感染的两个人，另外两个人将感染疾病，使新感染者的数量增加到四个。同样，在第 4 天，会出现八个新病例，依此类推。我们可以看到疾病扩展的速度并不是恒定的，因为新病例的数量取决于给定时刻现有病例的数量，这解释了大流行是如何呈指数增长和传播的。

前面的数字（1, 2, 4, 8...）形成一个序列。请注意，现在，算术序列的要求没有得到满足：两个连续项之间的差异不是恒定的。然而，比率是恒定的。这说明了前面的序列是一种特殊类型的序列，称为几何，并且被定义为一个序列或一组有序数字，其中任意两个连续项的比率是恒定的。

在数学的简洁语言中，我们可以将前面的行为写成αn+1 = r αn。

在这里，αn 是第n天的病例数量，αn+1 是第n+1天的新病例数量，r>0是定义增长速度（或减速速度）的系数。这被称为公比。前面的公式是通用的，这意味着它适用于所有成员n。因此，如果它对n成立，那么它对n-1，n-2等也成立。通过递归地处理前面的关系，我们可以轻松地得到αn = rn-1α方程。

在这里，我们给出了几何序列的第n项，一旦给出了第一项α=α1 和公比r。术语α被称为比例因子。

请注意，r可以有任何非零值。如果r>1，每一代，αn+1 都比前一代大，因此序列是不断增加的，而如果r<1，则相反：随着n的增加，αn+1 趋于零。因此，在传染病的初始例子中，r>1意味着传播是增加的，而r<1则导致传播减少。

让我们编写一个 Python 函数，根据αn = rn-1α的公式计算几何函数的第n项：

def n_geom_seq(r, a, n):
    an = r**(n-1) * a 
    return an 

该函数中的输入是r，公比，a，比例因子，和n，我们想要找到的第n项。让我们用一些参数调用这个函数，(2, 3, 10)：

n_geom_seq(2, 3, 10) 

输出如下：

1536

同样，对于算术序列的情况，我们将几何级数定义为序列长度n的项的总和：

图 5.3：几何序列

或者，我们可以这样表达：

图 5.4：几何序列的另一种表达

为了更好地理解几何级数，让我们看看它在 Python 中是如何工作并进行可视化的。我们需要定义一个函数，接受r、a和n（就像之前一样）作为输入，并计算第二个公式，即到第n项的级数：

def sum_n(r, a, n):
    sum_n = a*(1 - r**n) / (1 - r) 
    return sum_n

现在，像之前一样，为参数(2, 3, 10)调用函数：

sum_n(2, 3, 10)

输出如下：

3069.0

看一下几何序列的以下示例图，其中值对r>1递增：

图 5.5：几何序列 r>1 递增

看一下几何序列的以下示例图，其中值对r<1递减：

图 5.6：几何序列 r<1 递减

在本节中，我们已经看到了几何序列的进展以及如何在 Python 中轻松找到它的项，以及几何级数。我们现在准备在练习中实现我们所学到的内容，以便更好地理解序列及其应用。

练习 5.02：编写一个函数来找到序列的下一个项

培养皿中细菌数量以几何序列增加。给定每天的细菌数量，跨越一定天数n，编写一个函数，计算第n+1天的细菌数量。按照以下步骤完成这个练习：

1. 编写一个函数，接受可变数量的参数（*args）并计算任何元素与其前一个元素之间的比率（从第二个元素开始）。然后，检查找到的所有比率是否相同，并返回它们的唯一值。否则，函数返回-1（序列没有唯一的公共比率）：

def find_ratio(*args):
    arg0=args[0]
    ratios = []
    for arg in args[1:]:
        ratio = round(arg/arg0,8)
        arg0=arg
        ratios.append(ratio)
    if len(set(ratios)) == 1:
        return ratio
    else:
        return -1

1. 现在，检查find_ratio函数的两种不同情况。首先，让我们使用以下序列：

find_ratio(1,2,4,8,16,32,64,128,256,512)

输出如下：

2.0

1. 现在，让我们使用以下序列：

find_ratio(1,2,3)

输出如下：

-1

如前面的输出所示，find_ratio函数打印出比率，如果存在的话，或者如果序列不是几何序列，则打印-1。

1. 现在，编写第二个函数，读取一个序列并打印出下一个项将是什么。为此，读取一个（逗号分隔的）数字列表，找到它们的比率，然后预测下一个项：

def find_next(*args):
    if find_ratio(*args) == -1:
        raise ValueError('The sequence you entered' \
                         'is not a geometric sequence. '\
                         'Please check input.')
    else:
        return args[-1]*find_ratio(*args)

请注意，我们要通过调用我们之前编写的find_ratio()函数来检查序列是否具有公共比率。如果没有，就引发一个错误；如果有，就找到下一个项并返回它。

1. 通过使用以下序列来检查它是否有效：

find_next(1,2,4)

以下是前面代码的输出：

8.0

1. 现在，尝试使用不同的序列：

find_next(1.36,0.85680,0.539784,0.34006392)

输出如下：

0.2142402696

它有效。在第一种情况下，明显的结果8.0被打印出来。在第二种情况下，找到并打印出了递减的几何序列的不太明显的结果。总之，我们能够编写一个函数，检测几何序列，找到它的比率，并用它来预测下一个序列项。这在现实生活中非常有用，比如需要验证复利利率的情况。

注意

要访问此特定部分的源代码，请参阅packt.live/2NUyT8N。

您也可以在packt.live/3dRMwQV上在线运行此示例。

在前面的章节中，我们看到序列，无论是等差还是等比，都可以用两种等价的方式来定义。我们看到序列的第n项是通过知道序列的给定项（通常是第一项，但不一定）和公差或公比来确定的。更有趣的是，我们看到序列的第n项可以通过知道第(n-1)项来找到，而第(n-1)项又可以通过知道第(n-2)项来找到，依此类推。因此，这里有一个有趣的模式，决定了我们研究的两种序列类型，实际上超出了它们。事实证明，我们可以概括这种行为，并以纯粹递归的方式定义序列，而不一定是等差或等比的。现在，让我们继续下一节，我们将了解递归序列。

递归序列

递归序列是一个元素序列，υn，通过递归关系产生，也就是说，每个元素都是唯一地从前面的元素产生出来的。

υn 可以依赖于一个或多个在它之前的元素。例如，我们在本章前面看到的斐波那契数列是一个递归序列，其中知道第n项需要知道第(n-1)项和第(n-2)项。另一方面，阶乘只需要前面的元素。具体来说，它由递归关系定义，n! = n(n-1)! , n > 0，和初始条件，0! = 1。

让我们将前面的公式转换成 Python 代码：

def factorial(n):
    if n == 0 or n ==1:
        return 1
    elif n == 2:
        return 2
    else:
        return n*factorial(n - 1)

前面的代码是阶乘函数的递归实现：为了计算n的结果，我们调用n-1的函数，然后n-1调用n-2的函数，依此类推，直到n=2。

如果我们对n=11的情况执行前面的函数，我们得到以下结果：

factorial(11)

输出如下：

39916800

请注意，到目前为止我们看到的前两类序列（等差和等比）是互斥的，但是递归序列不是，也就是说，序列既可以是递归的，也可以是等差的或等比的。按照惯例，我们用术语递归来表示这些类型的序列，与等比和等差不同，它们不能以非递归的方式表示。

现在我们已经探讨了递归序列的基本概念，我们可以在 Python 中实现这一点，并编写计算递归定义的任何序列的任意数量的元素的代码。

练习 5.03：创建自定义递归序列

在这个练习中，我们将使用我们在前一节中解释的概念创建一个自定义的递归序列。给定序列的前三个元素，Pn，也就是P1=1，P2=7，和P3=2，找出递归定义的序列的接下来七个项，关系为Pn+3= (3Pn+1 - Pn+2)/(P*n – 1)。按照以下步骤完成这个练习：

1. 定义一个 Python 函数，它是递归的，并实现了前面给出的序列的第n个元素的关系：

def p_n(n):
    if n < 1:
        return -1
    elif n == 1:
        return 1
    elif n == 2:
        return 7
    elif n == 3:
        return 2
    else:
        pn = (3*p_n(n-2) - p_n(n-1) )/ (p_n(n-3) + 1)
        return pn

在这里，我们首先定义了基本情况，也就是在简介中给出的已知结果：如果n=1，那么P=1，如果n=2，那么P=7，如果n=3，那么P=2。我们还包括了n<1的情况。这是无效的输入，惯例上，我们的函数返回值为-1。这使得我们的函数有界，并且受到保护，不会进入无限循环和无效输入。一旦处理了这些情况，我们就定义了递归关系。

1. 现在，让我们测试我们的函数，并打印出序列的前 10 个值（三个对应于基本情况，另外七个是我们要找的）：

for i in range(1,11):
    print(p_n(i))

输出如下：

1
7
2
9.5
-0.4375
9.645833333333334
-1.0436507936507937
53.29982363315697
-5.30073825572847
-3784.586609737289

从前面的输出中可以看出，我们的函数有效，并返回了序列的已知值（P1 = 1，P2 = 7，和P3 = 2）以及我们正在寻找的下一个项（P_1到P_10）。

1. 作为奖励，让我们使用matplotlib模块绘制我们的发现。我们将创建一个包含序列前九个值的列表，然后用pyplot绘制它：

from matplotlib import pyplot as plt
plist = []
for i in range(1,40):
    plist.append(p_n(i))

plt.plot(plist, linestyle='--', marker='o', color='b')
plt.show()

输出如下：

图 5.7：使用 matplotlib 库创建的图表

注意

要访问本节的源代码，请参阅packt.live/2D3vlPF。

您也可以在packt.live/3eY05Q4上在线运行此示例。

我们可以看到，一个简单而明确定义的递归关系可以导致明显随机或混沌的结果。实际上，如果您继续绘制前述序列的项，很快就会注意到在项的模式中没有明显的规律，它们在 0 周围广泛且不对称地振荡。这促使我们得出结论，即使定义了一个递归序列并预测了它的第 n 项是直接的，相反的情况并不总是成立。正如我们所看到的，鉴于一个序列（一系列数字），检查它是否形成等差数列、等比数列或两者都不是，是非常简单的。然而，要回答一个给定的序列是否由递归关系导出——更不用说这个递归是什么——是一个非平凡的任务，在大多数情况下都无法回答。

在本节中，我们介绍了序列是什么，为什么它们很重要，以及它们与数学中另一个重要概念的联系：级数。我们研究了三种一般类型的序列，即等差、等比和递归，并看到它们如何可以在 Python 中通过几个简单的步骤实现。在下一节中，我们将深入研究三角学，并学习如何使用 Python 轻松解决三角学问题。

三角学

三角学是研究三角形，特别是它们的角与边的关系。三角形的三条边（边）中的两条的比值提供了关于特定角的信息，并且对于这样一对边，我们给它一个特定的名称，并称之为函数。三角学和数学的美妙之处在于这些函数，它们诞生于三角形内部，在三角形不存在的任何其他情况下都有（抽象的）意义，并且作为独立的数学对象运行。因此，诸如正切、余弦和正弦之类的函数在大多数数学、物理和工程领域都可以找到，而无需参考三角形。

让我们看看最基本的三角函数及其用法。

基本三角函数

我们将从定义一个直角三角形（或简称直角三角形）开始，三角形 ABC。它的一个角（下图中的角 BCA）是一个直角，即 90 度角。直角的对边称为斜边（下图中的边h），而其他两边（a和b）称为腿。它们也被称为相对于各自角的对边和邻边。例如，边b是相邻于下图中的右下角的（角 CAB 或θ），而当我们提到顶角（角 CBA）时，它是对边：

图 5.8：直角三角形

最常见的三角函数是通过前面的图表定义的，并且定义如下：

图 5.9：三角函数

对于正切函数，也成立 tanθ = sinθ/cosθ。

此外，对于任何角度θ，以下恒等式始终成立：sinθ2 + cosθ2 = 1。

根据构造，三角函数是周期性的。这意味着，无论三角形的边的大小如何，前述函数都会在每 2π重复一次。这将在下一个练习中变得明显，我们将在其中绘制它们。正弦和余弦函数的范围是区间[-1,1]。这意味着它们可以获得的最小值是-1，最大值是 1，无论输入θ是什么。

最后但并非最不重要的，直角三角形的边是通过著名的毕达哥拉斯定理连接的：h2 = a2 + b2

在 Python 代码中，毕达哥拉斯定理的一个简单实现是编写一个函数，利用math模块的平方根（sqrt）方法，计算h，给定a和b；例如：

from math import sqrt
def hypotenuse(a,b):
    h = sqrt(a**2 + b**2)
    return h

对a=3和b=4调用此函数会给我们以下结果：

hypotenuse(a = 3, b = 4)

输出如下：

5.0

现在，让我们看一些具体的例子，以便我们能够掌握这些想法。

练习 5.04：绘制直角三角形

在这个练习中，我们将编写 Python 函数，用于绘制给定点p1 和p2*的直角三角形。直角三角形将对应于三角形腿的端点。我们还将计算非直角的三角函数。让我们绘制基本的三角函数：

1. 导入numpy和pyplot库：

import numpy as np
from matplotlib import pyplot as plt

现在，编写一个函数，当给定两边p1 和p2 时，使用毕达哥拉斯定理返回斜边：

def find_hypotenuse(p1, p2):
    p3 = round( (p1**2 + p2**2)**0.5, 8)
    return p3

1. 现在，让我们编写另一个函数，实现sin、cos和tan函数的关系。输入是给定角的邻边、对边和斜边的长度，结果是三角函数值的元组：

def find_trig(adjacent, opposite, hypotenuse):
    '''Returns the tuple (sin, cos, tan)'''
    return opposite/hypotenuse, adjacent/hypotenuse, \
           opposite/adjacent

1. 现在，编写绘制三角形的函数。为简单起见，将直角放在坐标轴的原点（0,0），第一个输入点沿x轴放在（p1,0），第二个输入点沿y轴放在（0,p2）：

def plot_triangle(p1, p2, lw=5):
    x = [0, p1, 0]
    y = [0, 0, p2]
    n = ['0', 'p1', 'p2']
    fig, ax = plt.subplots(figsize=(p1,p2))
    # plot points
    ax.scatter(x, y,  s=400, c="#8C4799", alpha=0.4)
    ax.annotate(find_hypotenuse(p1,p2),(p1/2,p2/2))

    # plot edges
    ax.plot([0, p1], [0, 0], lw=lw, color='r')
    ax.plot([0, 0], [0, p2], lw=lw, color='b')
    ax.plot([0, p1], [p2, 0], lw=lw, color='y')
    for i, txt in enumerate(n):
        ax.annotate(txt, (x[i], y[i]), va='center')

在这里，我们创建了包含点的列表x和y，还有一个标签列表n。然后，我们创建了一个pyplot对象，首先绘制点，然后绘制边。最后两行用于注释我们的绘图；即，在我们的点旁边添加标签（从列表n中）。

1. 我们需要选择两个点来定义一个三角形。然后，我们需要调用我们的函数来显示绘图：

p01 = 4
p02 = 4
print(find_trig(p01,p02,find_hypotenuse(p01,p02)))
plot_triangle(p01,p02)

第一行打印三角函数sin、cos和tan的值。然后，我们绘制我们的三角形，在这种情况下是等腰三角形，因为它有两条相等长度的边。

输出将如下所示：

图 5.10：绘制等腰三角形

结果是预期的和正确的——在四舍五入误差后——因为这种特定形状的几何形状很简单（一个等腰直角三角形，两个角相等于π/4）。然后，我们检查了结果（请注意，在 NumPy 中，π的值可以直接调用np.pi）。

1. 最后，为了对sin和cos三角函数有一个总体概述，让我们绘制它们：

x = np.linspace(0, 10, 200)
sin = np.sin(x)
cos = np.cos(x)
plt.xticks([0, np.pi/2, np.pi, 3*np.pi/2, 2*np.pi, \
            5*np.pi/2, 3*np.pi], \
           ['0','','\u03C0','','2\u03C0','','3\u03C0'])
plt.plot(x, sin, marker='o', label='sin')
plt.plot(x, cos, marker='x', label='cos')
plt.legend(loc="upper left")
plt.ylim(-1.1, 1.6)
plt.show()

输出将如下所示：

图 5.11：正弦和余弦三角函数的绘图

在这个练习中，我们启动了对三角学领域的探索，并看到如何在 Python 中得到有用的可视化。

注意

要访问此特定部分的源代码，请参阅packt.live/2Zz0TnU

您也可以在packt.live/2AoxS63上在线运行此示例

有了这个，我们已经建立了主要的三角函数，并看到它们如何在角度和相关的三角函数值之间提供操作，由 sin、cos 或 tan 函数给出。此外，我们看到这三个函数是周期性的，即每 2π重复一次，而前两个是有界的，即它们可以取的值永远不会超出区间[-1,1]。这些值可以直接在 Python 或科学口袋计算器中找到。然而，在许多情况下，需要进行反向过程：如果我给出 sin、cos 或 tan 的值，我能找到角度吗？这样的函数存在吗？我们将在下一节中回答这些问题。

反三角函数

反三角函数是三角函数的反函数，与它们的对应函数一样有用。反函数是一个反转原始函数操作或结果的函数。请记住，三角函数接受角度作为输入值，并输出纯数（比值）。反三角函数则相反：它们接受纯数作为输入，并给出角度作为输出。因此，例如，如果一个点π被映射到点-1（如 cos 函数所做的那样），那么它的反函数需要完全相反。这种映射需要对反函数定义的每个点都成立。

sin(x)函数的反函数称为arcsin(x)：如果y=sin(x)，那么x=arcsin(y)。请记住sin是一个周期函数，因此许多不同的x被映射到相同的y。因此，反函数会将一个点映射到几个不同的点。这是不允许的，因为它与函数的定义相冲突。为了避免这个缺点，我们需要限制我们的arcsin（以及类似地arccos）的定义域为区间[-1,1]，而映射y=arcsin(x)和y=arccos(x)则限制在范围[-π/2,π/2]和[0, π]。

我们可以定义三个基本的反三角函数如下：

• arcsin(x) = y，使得 arcsin(sin(x)) = x

• arccos(x) = y，使得 arccos(cos(x)) = x

• arctan(x) = y，使得 arctan(tan(x)) = x

在 Python 中，这些函数可以从math模块或numpy库中调用。由于大多数 Python 实现的三角反函数返回弧度，我们可能希望将结果转换为度数。我们可以通过将弧度乘以 180 然后除以π来实现这一点。

让我们看看这如何在代码中编写。请注意，输入x表示为-1 和 1 之间的纯数，而输出表示为弧度。让我们导入所需的库并声明x的值：

from math import acos, asin, atan, cos
x = 0.5

现在，要打印余弦的反函数，添加以下代码：

print(acos(x))

输出如下：

1.0471975511965979

要打印正弦的反函数，添加以下代码：

print(asin(x))

输出如下：

0.5235987755982989

要打印正切的反函数，添加以下代码：

print(atan(x))

输出如下：

0.4636476090008061

让我们尝试为acos函数添加一个超出范围[-1,1]的输入：

x = -1.2
print(acos(x))

我们将会得到一个错误，如下所示：

Traceback (most recent call last):
    File "<stdin>", line 1, in <module>

asin也会发生类似的情况。这是可以预料到的，因为不存在任何角度φ可以返回-1.2作为 cos（或 sin）。然而，这个输入在atan函数中是允许的：

x = -1.2
print(atan(x))

输出如下：

-0.8760580505981934

最后，让我们来看看反函数arccos(cos(x))的反函数给我们带来了什么：

print(acos(cos(0.2)))

输出如下：

0.2

正如预期的那样，我们检索到cos函数的输入值。

反三角函数在数学、物理和工程学中有各种应用。例如，可以使用反三角函数来计算积分。不定积分如下：

图 5.12：反三角函数

在这里，a是一个参数，C是一个常数，积分立即可以通过反三角函数得到解决。

练习 5.05：使用反三角函数找到通往宝藏的最短路径

在这个练习中，您将获得一张指向B的秘密地图，几个世纪以来一些宝贵的宝藏一直在那里。您在点A，指令很明确：您必须向南导航 20 公里，然后向西导航 33 公里，以便到达宝藏。然而，直线段AB是最短的。您需要找到地图上的角度θ，以便您的导航正确定位：

图 5.13：点 A、B 和 C 的图形表示

我们需要找到角θ，即线段AB和AC之间的角度。按照以下步骤进行：

1. 导入atan（arctan 或反正切）函数：

from math import atan, pi

1. 使用BC和AC找到θ的正切：

AC = 33
BC = 20
tan_th = BC/AC
print(tan_th)

输出如下：

0.6060606060606061

1. 然后，通过取反正切函数来找到角度。其参数是θ的正切：

theta = atan(tan_th)

1. 将其转换为度并打印出值：

theta_degrees = theta*180/pi
print(theta_degrees)

输出如下：

31.218402764346372

因此，答案是我们需要转动 31.22 度才能正确导航。

1. 作为奖励分，计算我们将沿着路径AB行进的距离。这只是由勾股定理给出的：

AB2 = AC2 + BC2

在 Python 中，使用以下代码：

AB = (AC**2 + BC**2)**0.5
print(AB)

输出如下：

38.58756276314948

课程将是 38.59 公里。

在 Python 中通过调用find_hypotenuse()函数很容易计算。正如预期的那样，这比路径AC + BC = 53公里要短得多。

注意

要访问此特定部分的源代码，请参考packt.live/31CF4qr。

您还可以在packt.live/38jfVlI上在线运行此示例。

练习 5.06：找到与对象的最佳距离

您正在参观当地的竞技场观看您最喜欢的节目，您站在竞技场的中间。除了主舞台外，还有一个观看屏幕，这样人们就可以观看并不会错过节目的细节。屏幕底部距离您的眼睛高 3 米，屏幕本身高 7 米。视野角是通过观看屏幕的底部和顶部形成的。找到最佳距离x，使视野角最大化：

图 5.14：眼睛和屏幕之间形成的视野角

这是一个稍微复杂的问题，需要一些代数，但我们将把它分解成简单的步骤并解释逻辑。首先注意问题的情节如何引导我们并帮助我们找到解决方案。这个明显复杂的现实世界问题转化为一个更抽象和简单的几何图像。按照以下步骤完成此练习：

1. 计算x。这是三角形的底边，也是角θ1（也是θ=θ1+θ2）的邻边。答案x将由观看角θ2 或等效地，tan(θ2)最大化的条件给出。从屏幕的前述图中，我们可以立即得出三个角度的以下关系：θ1（内角）、θ2（外角）和θ=θ1+θ2：

tan(θ1) = 对边/邻边 = 3/x

tan(θ) = tan(θ1+θ2) = 对边/邻边 = (7+3)/x .

现在，使用代数来处理这两个关系，并得到θ2 的条件。

1. 两个角的正切和的已知身份如下：

图 5.15：两个角的正切的公式

通过在后一个关系中代入我们找到的tan(θ)和tan(θ1)，并经过代数运算，我们得到以下结果：

tan(θ2) = 7x/(30+x2) 或

θ2 = arctan(7x/(30+x2)).

换句话说，我们已经结合了问题的要素，并发现角度θ1 应该随着距离x的变化而变化，这是前一行给出的函数x的函数。

1. 让我们绘制这个函数，看看它是如何变化的。首先，加载必要的库：

from matplotlib import pyplot as plt
import numpy as np

1. 然后，通过使用numpy的arctan方法，通过定义域x和值y来绘制函数。这些可以使用pyplot的plot()方法轻松绘制，如下所示：

x = np.linspace(0.1, 50, 2000)
y = np.arctan(7*x / (30+x**2) )
plt.plot(x,y)
plt.show()

输出将如下所示：

图 5.16：使用 arctan 方法绘制函数的图形

从前面的图表中，我们可以看到函数获得了最大值。

1. 确定函数的最大值y和位置x，以及发生这种情况的位置：

ymax = max(y)
xmax = x[list(y).index(ymax)]
print(round(xmax,2), round(ymax,2))

输出如下：

5.47 0.57

1. 最后，将找到的角度转换为度数：

ymax_degrees = round(ymax * 180 / np.pi, 2)
print(ymax_degrees)

输出如下：

32.58

因此，观察角度θ2 在 32.58 度时达到最大值，并且当我们站在距屏幕 5.47 米的地方时发生。我们使用三角函数和反三角函数，在 Python 中实现它们，并找到了一个来自现实生活中几何设置的问题的答案。这更加清楚地说明了几何和三角学概念如何被有用地和轻松地编码以提供预期的结果。

注意

要访问此特定部分的源代码，请参阅packt.live/2VB3Oez

您也可以在packt.live/2VG9x2T上在线运行此示例。

现在，我们将继续研究数学中的另一个核心概念，它在代数、物理学、计算机科学和应用数据科学中有着广泛的应用：向量。

向量

向量是具有大小（大小）和方向（方向）的抽象数学对象。向量由一个箭头表示，它有一个基（尾部）和一个头。箭头的头部显示向量的方向，而箭头的长度显示它的大小。

标量与向量相反，是一个单独的数字。它是一个非向量，即一个纯整数、实数或复数（我们稍后会看到），它没有元素，因此没有方向。

向量通常用粗体字母A、带箭头的字母或普通字母表示，如果在讨论中关于符号的表示没有歧义。向量A的大小被写成|A|或简单地写成A。现在，让我们来看看各种向量运算。

向量运算

简而言之，向量是由两个、三个或更多个数字组成的集合（可以想象成列表或数组），形成一个数学对象。这个对象存在于一个特定的几何空间中，称为向量空间，具有一些属性，如度量属性和维度。向量空间可以是二维的（想象一下你的书页上的平面），三维的（我们周围的普通欧几里德空间），或者在数学和物理学中的许多抽象情况下更高维度的。用于识别向量的元素或数字等于空间的维度。现在我们已经定义了一个向量空间——向量的游乐场——我们可以用一组坐标轴（通常的x、y和z轴）来装备它，标记原点并测量空间。在这样一个明确定义的空间中，我们需要确定一组数字（两个、三个或更多）来唯一定义一个向量，因为向量被假定从坐标轴的原点开始。向量的元素可以是整数、有理数、实数或（很少）复数。在 Python 中，它们通常由列表或 NumPy 数组表示。

与实数类似，向量上定义了一组线性运算。在两个向量 A = (a1, a2, a3)和 B = (b1, b2, b3)之间，我们可以定义以下内容：

图 5.17：执行向量操作时点 A、B 和 C 及其关系

现在让我们看看可以对这些向量执行的各种操作：

• 加法作为导致向量C = A + B = (a1 + b1, a2 + b2, a3 + b3)的操作。

• 减法作为导致向量C = A - B = (a1 - b1, a2 - b2, a3 - b3)的操作。

• 点积（或内积或标量）C = b. b = a1 b1 + a2 b2 + a3 b3。

• 向量C = A x B的叉积（或外积），它垂直于由A和B定义的平面，并具有元素（a2b3 - a3b2, a3b1 - a1b3, a1b2 – a2b1）。

• 向量A和B的逐元素或 Hadamard 乘积是向量C，其元素是A和B的元素的成对乘积；即C = (a1 b1, a2 b2, a3 b3)*。

我们可以在 Python 代码中定义并使用前述公式如下：

import numpy as np
A = np.array([1,2,3]) # create vector A
B = np.array([4,5,6]) # create vector B

然后，要找到A和B的和，输入以下代码：

A + B

输出如下：

array([5, 7, 9])

要计算差异，输入以下代码：

A - B

输出如下：

array([-3, -3, -3])

要找到逐元素乘积，输入以下代码：

A*B

输出如下：

array([ 4, 10, 18])

要找到点积，使用以下代码：

A.dot(B)

输出如下：

32

最后，叉积可以计算如下：

np.cross(A,B)

输出如下：

array([-3,  6, -3])

请注意，向量加法、减法和点积是可结合和可交换的操作，而叉积是可结合但不可交换的。换句话说，a x b 不等于 b x a，而是 b x a，这就是为什么它被称为反交换。

此外，向量A可以乘以标量λ。在这种情况下，您只需将每个向量元素乘以相同的数字，即标量：λ A = λ (a1, a2, a3) = (λ a1, λ a2, λ a3)

向量之间的另一个重要操作是点积，因为它可以说是数学、计算机科学及其应用中最常见的操作。点积是一种有趣的操作，它在实数领域中没有类似物。实际上，它需要两个向量作为输入，以产生单个标量作为输出。这意味着操作的结果（标量）与其成分（向量）的类型不同，因此通常不存在逆操作（点除法）。

根据定义，它如下所示：

图 5.18：θ角的图形表示

这可以用以下方程表示：

A.B = |A| |B| cos(θ)

这里，θ是向量A和B之间的角度。

让我们看一些典型的情况：

• 如果A和B是正交的，则点积消失：

如果且仅当θ = angle(A,B) = π/2时，A.B = 0，因为|A|和|B|不为零。

• 如果A和B共线且共方向，则θ = 0，cos(θ)=1和A.B = |A| |B|。如果它们共线且方向相反，则θ = π，cos(θ)=-1，A.B = -|A| |B|。

• 它遵循对向量与自身的点积的定义：A.A = |A| |A|或|A| = √(A.A)

• 它直接遵循A.B = |A| |B| cos(θ)，其中两个向量之间的角度如下给出：θ = arccos(A.B / |A| |B|)

这里，arccos是我们在前一节中看到的反余弦函数。

例如，我们可以编写一个 Python 程序，利用numpy和前面给出的关系计算任意两个给定向量之间的角度θ：

import numpy as np
from math import acos
A = np.array([2,10,0])
B = np.array([9,1,-1])

要找到每个向量的范数（大小），我们可以使用以下代码：

Amagn = np.sqrt(A.dot(A))
Bmagn = np.sqrt(B.dot(B))

作为替代方案，您还可以使用以下代码：

Amagn = np.linalg.norm(A)
Bmagn = np.linalg.norm(B)

打印它们的值：

print(Amagn, Bmagn)

您将得到以下输出：

10.198039027185569
9.1104335791443

这两种替代方案都会导致相同的结果，您可以通过再次打印Amagn和Bmagn来立即检查。

最后，我们可以按如下方式找到角度θ：

theta = acos(A.dot(B) / (Amagn * Bmagn))
print(theta)

输出如下：

1.2646655256233297

现在，让我们看一个练习，我们将执行刚学到的各种向量运算。

练习 5.07：可视化向量

在这个练习中，我们将编写一个在二维空间中绘制两个向量的函数。我们将不得不找到它们的和以及它们之间的角度。

执行以下步骤完成这个练习：

1. 导入必要的库，即numpy和matplotlib：

import numpy as np
import matplotlib.pyplot as plt

1. 创建一个函数，接受两个向量作为输入，每个向量作为一个列表，绘制它们，并可选择绘制它们的和向量：

def plot_vectors(vec1, vec2, isSum = False):

    label1 = "A"; label2 = "B"; label3 = "A+B"
    orig = [0.0, 0.0]  # position of origin of axes

vec1和vec2列表中分别包含两个实数。每对数字表示相应向量的端点（头部）坐标，而原点设置为（0,0）。标签设置为"A"、"B"和"A+B"，但您可以更改它们，甚至将它们设置为plot_vectors函数的变量（或不带）默认值。布尔变量isSum默认设置为False，和vec1+vec2将不会被绘制，除非显式设置为True。

1. 接下来，我们将坐标放在一个matplotlib.pyplot对象上：

    ax = plt.axes()
    ax.annotate(label1, [vec1[0]+0.5,vec1[1]+0.5] )   
    # shift position of label for better visibility
    ax.annotate(label2, [vec2[0]+0.5,vec2[1]+0.5] )
    if isSum: 
        vec3 = [vec1[0]+vec2[0], vec1[1]+vec2[1]]     
        # if isSum=True calculate the sum of the two vectors
        ax.annotate(label3, [vec3[0]+0.5,vec3[1]+0.5] )

    ax.arrow(*orig, *vec1, head_width=0.4, head_length=0.65)
    ax.arrow(*orig, *vec2, head_width=0.4, head_length=0.65, \
             ec='blue')
    if isSum:
        ax.arrow(*orig, *vec3, head_width=0.2, \
                 head_length=0.25, ec='yellow')
        # plot the vector sum as well

    plt.grid()
    e=3 
    # shift limits by e for better visibility
    plt.xlim(min(vec1[0],vec2[0],0)-e, max(vec1[0],\
                 vec2[0],0)+e) 
    # set plot limits to the min/max of coordinates
    plt.ylim(min(vec1[1],vec2[1],0)-e, max(vec1[1],\
                 vec2[1],0)+e) 
    # so that all vectors are inside the plot area

在这里，我们使用 annotate 方法向我们的向量添加标签，以及 arrow 方法来创建我们的向量。星号运算符*用于解包列表orig和vec1、vec2中的参数，以便它们可以从arrow()方法中正确读取。plt.grid()在绘图的背景上创建一个网格，以引导眼睛，这是可选的。添加e参数是为了使绘图限制足够宽，绘图可读。

1. 接下来，给我们的图表加上标题并绘制它：

    plt.title('Vector sum',fontsize=14)
    plt.show()
    plt.close()

1. 现在，我们将编写一个函数，计算两个输入向量之间的角度，如前所述，借助点（内）积：

def find_angle(vec1, vec2, isRadians = True, isSum = False):
    vec1 = np.array(vec1)
    vec2 = np.array(vec2)

    product12 = np.dot(vec1,vec2)
    cos_theta = product12/(np.dot(vec1,vec1)**0.5 * \
                           np.dot(vec2,vec2)**0.5 )
    cos_theta = round(cos_theta, 12)
    theta = np.arccos(cos_theta)

    plot_vectors(vec1, vec2, isSum=isSum)
    if isRadians:
        return theta
    else:
        return 180*theta/np.pi

首先，我们将我们的输入列表映射到numpy数组，以便我们可以使用这个模块的方法。我们计算点积（命名为product12），然后将其除以vec1的大小与vec2的大小的乘积。回想一下，向量的大小由其与自身的点积的平方根（或**0.5）给出。根据点积的定义，我们知道这个数量是两个向量之间角度的余弦。最后，在四舍五入余弦以避免输入错误后，利用numpy的arccos方法计算角度。

1. 我们希望将我们编写的两个函数find_angle和plot_vectors结合起来，并在后者内部调用前者。我们还希望给用户选择以弧度（isRadians=True）或度数（isRadians=False）打印角度结果的选项。现在我们准备尝试我们的函数。首先，让我们尝试使用两个垂直向量：

ve1 = [1,5]
ve2 = [5,-1]
find_angle(ve1, ve2, isRadians = False, isSum = True)

输出如下：

图 5.19：两个垂直向量的绘图

图看起来不错，结果是 90 度，如预期。

1. 现在，让我们尝试使用相同的函数创建两个共线向量：

ve1 = [1,5]
ve2 = [0.5,2.5]
find_angle(ve1, ve2, isRadians = False, isSum = True)

输出如下：

图 5.20：两个共线向量的绘图

输出为 0 度，如预期。

1. 最后，再次使用相同的函数，让我们创建两个通用向量：

ve1 = [1,5]
ve2 = [-3,-5]
find_angle(ve1, ve2, isRadians = False, isSum = True)

输出如下：

图 5.21：两个通用向量的绘图

总之，我们已经学习了向量作为生活在向量空间中的数学对象。我们学会了如何在 Python 中构造和表示向量以及如何可视化它们。向量遵循一些简单的规则，并且可以进行操作。在处理实数时，加法和减法遵循完全相同的逻辑。乘法有些更复杂，并且定义了不同类型的乘积。最常见的乘积是内积或点积，由于其简单的几何表示，在数学和物理界广受欢迎。我们学会了如何在 Python 中计算任意两个向量的点积，并且利用我们对点积的知识（以及一些 NumPy 方法）找到了这对向量之间的角度。简而言之，在二维空间中，向量是一对形成有趣属性的数字。

注意

要访问此特定部分的源代码，请参阅packt.live/2Zxu7n5.

您还可以在packt.live/2YPntJQ.上在线运行此示例

接下来，我们将学习如何将一对两个数字组合成一个更令人兴奋的对象，即复数。

复数

自古代数字系统以来，数学思想一直在发展，关于数字及其关系的数学思想也在历史上从具体到抽象不断演变。例如，自然数集合的概念是为了让我们周围世界中的所有物体直接对应于该集合中的某个数字。随着算术和代数的发展，人们意识到除了自然数或整数之外，还需要小数和有理数，因此引入了小数和有理数。同样，在毕达哥拉斯时代，人们发现有理数无法解决我们当时所知的几何构造的所有数学问题。这是因为引入了无理数——从其他数字的平方根得出并且没有比率表示的数字。

复数是实数的扩展，并包括一些特殊的数字，可以解决一些实数无法解决的方程。

这样的数字实际上是存在的，并且有符号i。它被称为虚数或虚数单位，尽管它并不虚构；它和我们见过的所有其他数字一样真实，并且正如我们将看到的那样，具有一些非常美丽的性质。

复数的基本定义

我们将虚数i定义如下：

i2 = -1

任何由实数和虚数（部分）组成的数字都称为复数。例如，考虑以下数字：

z = 3 – i

z = 14/11 + i 3

z = -√5 – i 2.1

所有前述的数字都是复数。它们的实部表示为Re(z)，虚部表示为Im(z)。对于前述的例子，我们得到以下结果：

Re(z) = 3 , Im(z) = -1

Re(z) = 14/11 , Im(z) = 3

Re(z) = -√5 , Im(z) = -2.1

让我们看一些使用代码的例子。在 Python 中，虚数单位用字母j表示，复数写成如下形式：

c = <real> + <imag>*1j,

这里，<real>和<imag>都是实数。同样，复数可以定义如下：

c = complex(<real>, <imag>).

在代码中，它变成了如下形式：

a = 1
b = -3
z = complex(a, b)
print(z)

输出如下：

(1-3j)

我们还可以使用real和imag函数来分离任何复数z的实部和虚部。首先，让我们使用real函数：

print(z.real)

输出如下：

1.0

现在，使用imag函数：

print(z.imag)

输出如下：

-3.0

换句话说，任何复数都可以被分解并写成z=Re(z) + i Im(z)。因此，一个复数是两个实数的一对，并且可以被视为生活在二维空间中的向量。因此，向量的几何和代数，如前一节所讨论的，也可以在这里应用。

接受复数作为输入的方法和函数可以在cmath模块中找到。该模块包含复数的数学函数。那里的函数接受整数、浮点数或复数作为输入参数。

复共轭被定义为与复数z具有相同实部和相反虚部的复数z**（也z̄）；也就是说，如果z = x+iy，那么z* = x -iy。注意，乘积zz**是实数x2+y2，这给出了z的模的平方：

zz = zz = |z|2

复数被绘制在复平面上，类似于向量（如下图所示）。这是由实部在x轴上和虚部在y轴上形成的平面。复共轭只是相对于实轴的向量的反射：

图 5.22：复数的绘图

复数z可以被视为具有坐标(x, y)的向量。或者，我们可以将其写成具有极坐标(r, φ)的向量。复共轭z**或z̄是一个与z相同的向量，但相对于x*轴反射。

如果一个复数的实部和虚部都是零，那么这个复数就是零。可以对两个复数z = x+iy和w = u+iv执行以下操作：

• 加法：z+w = (x+u) + i(y+v)

• 减法：z-w = (x-u) + i(y-v)

• 乘法：z w = (x+iy)(u+iv) = (xu-yv) + i(xv + yu)

• 除法：z/w = (x+iy)/(u+iv) = (ux+vy)+i(uy-xv) / (u2+v2)

极坐标表示和欧拉公式

复数很容易被视为复平面上的向量。因此，它有一个大小，由向量的大小确定，以及一个方向，由与x（实）轴形成的角度φ确定。要确定这两个数，我们需要找到z=x+iy的绝对值（或模），r：

r = |z| = √x2+y2

它的角度（也称为参数arg或相位），φ，如下：

φ = arg(z) = arctan(x+iy) = arctan(y/x)

这两个关系都源自复向量的几何关系。第一个关系简单地是勾股定理的应用，而第二个来自对角度φ应用正切关系。

通过检查向量的图形表示（见前面的图），我们可以看到以下内容：

cos(φ) = x/r 和

sin(φ) = y/r

或

x = r cos(φ) 和

y = r sin(φ)

通过用z = x+iy替换这些，我们得到以下结果：

z = r (cos(φ) + i sin(φ))

我们可以在 Python 中编写一些代码，一旦给定了(x, y)(笛卡尔坐标)，就可以找到(r, φ)(极坐标)，反之亦然：

def find_polar(z):
    from math import asin
    x = z.real
    y = z.imag
    r = (x**2 + y**2)**0.5
    phi = asin(y/r)
    return r, phi
find_polar(1-3j)

输出如下：

(3.1622776601683795, -1.2490457723982544)

同样，我们可以使用cmath模块中的polar方法：

import cmath
z = 1-3j
cmath.polar(z)

输出如下：

(3.1622776601683795, -1.2490457723982544)

注意

不允许输入(0,0)，因为这会导致除以零。

因此，复数可以用它的模，r，和相位，φ，来表示，而不是用它的横坐标（x，实部）和纵坐标（y，虚部）。模，r，是一个实数，非负数，相位，φ，在区间[-π,π]内：对于纯实数，它是0和π，对于纯虚数，它是π/2或-π/2。后者的表示被称为极坐标，而前者被称为矩形或笛卡尔；它们是等价的。也可以用以下表示：

z = r eiφ = r (cos(φ) + i sin(φ))

这是自然对数的底。这被称为欧拉公式。特殊情况φ=π给出了以下结果：

eiπ + 1 = 0

这就是欧拉恒等式。

使用欧拉公式的好处在于，复数乘法和除法获得了简单的几何表示。要乘（除）两个复数z1 和z2，我们只需将它们各自的模相乘（除）并加上（减去）它们的幅角：

z1 ** z2 = r eiφ = r1 ** r2 ei(φ1+φ2)

现在，让我们在 Python 中使用一些复数进行数学运算。我们将编写两个复数的加法、减法、乘法和除法：

def complex_operations2(c1, c2):
    print('Addition =', c1 + c2)
    print('Subtraction =', c1 - c2)
    print('Multiplication =', c1 * c2)
    print('Division =', c1 / c2)

现在，让我们尝试这些函数对一对通用的复数c1=10+2j/3和c2=2.9+1j/3：

complex_operations2(10+2j/3, 2.9+1j/3)

输出如下：

Addition = (12.9+1j)
Subtraction = (7.1+0.3333333333333333j)
Multiplication = (28.77777777777778+5.266666666666666j)
Division = (3.429391054896336-0.16429782240187768j)

我们可以对纯实数和纯虚数做同样的操作：

complex_operations2(1, 1j)

输出如下：

Addition = (1+1j)
Subtraction = (1-1j)
Multiplication = 1j
Division = -1j

从最后一行，我们很容易看出1/i = -i，这与虚数单位的定义一致。cmath库还为复数提供了有用的函数，如phase和polar，以及带有复数参数的三角函数：

import cmath
def complex_operations1(c):
    modulus = abs(c)
    phase = cmath.phase(c)
    polar = cmath.polar(c)
    print('Modulus =', modulus)
    print('Phase =', phase)
    print('Polar Coordinates =', polar)
    print('Conjugate =',c.conjugate())
    print('Rectangular Coordinates =', \
           cmath.rect(modulus, phase))
complex_operations1(3+4j)

输出如下：

Modulus = 5.0
Phase = 0.9272952180016122
Polar Coordinates = (5.0, 0.9272952180016122)
Conjugate = (3-4j)
Rectangular Coordinates = (3.0000000000000004+3.9999999999999996j)

因此，计算给定复数的模、相位或共轭变得非常简单。请注意，最后一行给出了复数的直角（或笛卡尔）形式，给定其模和相位。

现在我们已经了解了复数的算术和表示方式，让我们继续看一个涉及逻辑的练习，并结合我们在前几节中使用和学习的内容。

练习 5.08：复数的条件乘法

在这个练习中，您将编写一个函数，该函数读取一个复数c，如果复数的幅角大于零，则将其乘以自身，如果幅角小于零，则取c的平方根，如果幅角等于零，则不执行任何操作。绘制并讨论您的发现：

1. 导入必要的库，并且可以选择地抑制任何警告（这不是必需的，但如果您希望保持输出整洁，这是有帮助的，因为警告取决于您使用的库的版本）：

import cmath
from matplotlib import pyplot as plt
import warnings
warnings.filterwarnings("ignore")

1. 现在，定义一个函数，该函数使用 Matplotlib 的pyplot函数来绘制输入复数c的向量：

def plot_complex(c, color='b', label=None):

    ax = plt.axes()
    ax.arrow(0, 0, c.real, c.imag, head_width=0.2, \
             head_length=0.3, color=color)
    ax.annotate(label, xy=(0.6*c.real, 1.15*c.imag))
    plt.xlim(-3,3)
    plt.ylim(-3,3)
    plt.grid(b=True, which='major') #<-- plot grid lines

1. 现在，创建一个函数，读取输入复数c，通过调用之前定义的函数绘制它，然后根据输入的相位研究不同的情况。绘制操作前后的相位以及结果，以便将结果向量与输入向量进行比较：

def mult_complex(c, label1='old', label2='new'):

    phase = cmath.phase(c)
    plot_complex(c, label=label1)

    if phase == 0:
        result = -1
    elif phase < 0:
        print('old phase:', phase)
        result = cmath.sqrt(c)
        print('new phase:', cmath.phase(result))
        plot_complex(result, 'red', label=label2)
    elif phase > 0:
        print('old phase:', phase)
        result = c*c
        print('new phase:', cmath.phase(result))
        plot_complex(result, 'red', label=label2)
    return result

请注意，对于负相位，我们取c的平方根（使用math.sqrt()方法），而对于正相位，我们取c的平方。

1. 现在，转换一个位于复平面上半部的数字：

mult_complex(1 + 1.2j)

输出如下：

图 5.23：位于复平面上半部的数字的绘图

这里，一个具有正幅角φ（蓝色向量）的复数正在被转换（或映射）为一个具有更大模的新复数（红色向量），并且新的幅角是先前值的两倍。这是预期的：记住欧拉公式c=r eiφ的极坐标表示？很明显，平方c2 是一个具有原始幅角φ和模r2 两倍的数字。

1. 接下来，转换一个位于复平面下半部的数字：

mult_complex(1-1.2j)

输出如下：

图 5.24：位于复平面下半部的数字的绘图

在这种情况下，计算了平方根。与第一个例子类似，新转换的向量的模是原始向量的模的平方根，幅角是原始幅角的一半。

注意

有趣的事实：在两种情况下，向量都是逆时针旋转的。

1. 编写一个while迭代，调用mult_complex()函数n次，以检查如果我们保持向量旋转会发生什么：

c0 = 1+1.2j
n = 0
while n < 6:
    c0 = mult_complex(c0, None, str(n))
    n+=1

输出如下：

图 5.25：旋转向量的绘图

通过这样，我们看到了向量和向量代数如何用于可视化几何运算。特别是，除法和乘法复数的结果是获得几何表示，这在处理大量数据和可视化时非常有用。

注意

要访问此特定部分的源代码，请参阅packt.live/31yU8W1。

您也可以在packt.live/2BXWJOw上在线运行此示例。

活动 5.01：使用级数计算您的退休计划

在许多国家，一些雇主提供退休计划（也称为 401(k)）。这些计划允许您直接从工资中捐款，因此是一种轻松有效的储蓄和投资退休的方式。您的任务是编写一些代码，根据捐款金额和持续时间计算和绘制您的月度回报。

退休计划会像等比数列一样随时间累积。这是一种投资：您每月存钱，以便以后每月收取，带有附加值或利息。计算退休回报的主要要素是您当前的余额，每月的捐款，雇主匹配（雇主的捐款），退休年龄，回报率（您从 401(k)投资中预期的平均年回报），预期寿命和任何其他费用。在现实情况下，会引入上限：雇主匹配（通常在 50%至 100%之间）不能超过您年薪的 6%。同样，员工的捐款在一年内不能超过一定金额（通常为 18K），无论工资有多高。

执行以下步骤以完成此活动：

1. 确定我们问题的变量。这些将是我们函数的变量。确保仔细阅读活动描述，并内化已知和要计算的内容。

2. 识别序列并编写一个函数，计算某一年n的退休计划价值。该函数应接受当前余额、年薪、年份n等输入，并返回每年的捐款、雇主匹配和退休总价值的元组。

3. 识别级数并编写一个函数，计算n年后退休计划的累积价值。当前函数应读取输入，调用前一个函数计算每年计划的价值，并对所有（每年）储蓄进行求和。为了可视化目的，应将捐款（每年）、雇主匹配（每年）和总价值（每年）作为列表返回到一个元组中。

4. 对各种选择的值运行函数，并确保其正常运行。

5. 使用 Matplotlib 绘制结果。

注意

此活动的解决方案可在第 672 页找到。

摘要

在本章中，您已经得到了关于序列、级数、三角学、向量和复数等最核心数学概念的一般性和有益的阐述，更重要的是，我们使用具体而简短的例子在 Python 中实现了它们。作为一个现实生活的例子，我们研究了一个退休计划以及我们储蓄随时间的变化。然而，许多其他情况可以模拟成序列或级数，并通过应用向量或复数进行研究。这些概念和方法在物理学、工程学、数据科学等领域被广泛应用。线性代数，即对向量、矩阵和张量的研究，严重依赖于对几何和向量概念的理解，并几乎无处不在地出现在数据科学和神经网络的研究中。另一方面，几何和三角学明确地用于模拟物理运动（例如在视频游戏中）以及地理空间应用中的更高级概念。然而，对这些概念有背景知识使得使用和应用数据科学方法更加具体和可理解。

在下一章中，我们将讨论矩阵以及如何将它们应用于解决现实世界的问题。我们还将研究马尔可夫链，它们用于将与概率、矩阵和极限相关的概念联系在一起。

NDN74

ETB65

第六章：使用 Python 进行矩阵和马尔可夫链

概述

在本章中，我们打算使用 Python 深入研究矩阵及其应用。我们将探讨不同的矩阵操作技术，这将帮助我们有效地在实际应用中使用它们来构建有用的工具。

在本章结束时，您将了解矩阵并能对其进行操作。您将使用转移矩阵实现矩阵的一种应用，称为马尔可夫链，然后使用马尔可夫链和马尔可夫性质来解决现实世界问题。

介绍

矩阵是数字或表达式按行和列排列并被视为单个实体的矩形数组。由于我们将矩阵视为单个对象，如果我们对其进行操作，它将应用于其中的每个元素：

图 6.1：一个简单的 m×n 矩阵，其中 m 行 n 列

一个简单的线性单维数组很少能满足我们生活中与空间和时间相关的几乎所有属性都需要超过一个维度。紧凑性是使用矩阵的主要原因之一。当矩阵是封闭和有界的，或者简单地其点彼此之间的距离固定时，矩阵是紧凑的。主要基于这两个原因，矩阵在几乎每个领域都有应用，包括从图论、线性变换和概率论到物理学的不同分支，如量子力学和电磁学等基本数学概念。

马尔可夫链模型及其变体是这样一个应用，将矩阵、极限和概率的概念联系在一起，以产生在不确定性占主导地位的现实世界问题中的结果。在数学空间中，每当存在不确定性时，决策都是基于概率的；这构成了马尔可夫链的基础。这些使用一种特定类型的矩阵，称为转移矩阵，来构建状态图。马尔可夫链实际上是一个无记忆过程，主要基于当前状态来决定下一个状态的结果。我们在一些非常重要的用例中发现它们的应用，包括页面排名算法、自动完成应用程序和文本生成器。我们将在本章后面更详细地研究这些概念，为此，我们首先需要了解矩阵。

单矩阵上的矩阵操作

在本章中，我们将学习不同的矩阵操作方式以及如何在 Python 中实现它们。广义上理解矩阵的工作意味着理解二维或多维数组的基本原理。一旦我们对二维矩阵的基础有了很好的理解，有兴趣的人可以深入研究矩阵的高级研究，其中包括稀疏矩阵、向量空间、特征值和特征向量等特殊类型的矩阵，这可能涉及超过两个维度。

Python 中的矩阵可以使用列表或数组来实现。在 Python 中，嵌套列表可以很好地工作，但 Python 有一个强大的包，使矩阵实现变得更容易，称为 NumPy。SciPy 是另一个帮助进行矩阵操作的包，但通常更适合于更大的矩阵计算。在本章中，我们将在整个章节中使用这两个模块。

矩阵的基本操作

在这一点上，假设您已经安装了 Python 及其默认库以运行基本的 Python 程序。

一旦您安装好了这个包，让我们定义我们的第一个矩阵：

# Importing Numpy package
import numpy as np
# Initializing and printing matrix z
x = np.array([[1, 2], [3, 4]])
print(x)

这个矩阵与下面的矩阵z相同，只是更好地表示，并且在可能的情况下，作为一个良好的实践是可取的：

# Initializing and printing matrix z
z = np.array([[1, 2],\
              [3, 4]])
print(type(z))

请注意，在前面的代码中，我们已经打印了变量z的类型。你能猜到输出是什么吗？输出应该如下所示：

<type 'numpy.ndarray'>

ndarray是 NumPy 使用的标准数组格式。数组对象是同质的、多维的，并且具有与其分配相关联的数据类型对象。

例如，让我们取矩阵z中的一个元素，这个矩阵我们之前定义过：

# Printing the data types for matrix-z
print(z)
print(type(z))
print(z[0][1])
print(type(z[0][1]))

这将产生以下输出：

[[1 2]
 [3 4]]
[[5 6]
 [7 8]]
<type 'numpy.ndarray'>
2
<type 'numpy.int64'>
[Finished in 0.221s]

我们发现给定矩阵的元素是int64类型，即 64 位整数类型。其他数据类型包括np.float32、np.complex、np.bool、np.object、np.string_和np.unicode_。

目前，知道我们几乎每个数据结构都是使用 Python 版本 3.8 和 NumPy 版本 1.17 就足够了。截至出版日期，NumPy 有一个特殊的类叫做matrix类，它几乎做了与ndarray相同的事情。唯一的区别是matrix类保持其 2D 性质，并且内置了一些操作符，比如*表示乘法，**表示幂。虽然matrix类可能会派上用场并且可以被探索，但官方的 NumPy 文档建议使用普通的ndarray而不是np.matrix，因为后者可能会在将来被移除。因此，值得注意的是，在这个上下文中，术语ndarray可以被视为与术语matrix同义，并且在本章中可以互换使用。

让我们继续使用ndarray。假设我们有一个单独的矩阵，我们将看到一些我们可以对其进行的简单操作。我们可以使用之前定义的相同矩阵z。

让我们打印元素的总和：

# Sum of all elements of the matrix
print(z)
print(np.sum(z))

这将产生以下输出：

[[1 2]
 [3 4]]
10
[Finished in 0.237s]

这很简单。现在让我们看看我们可以做的其他事情。

让我们找到z矩阵的最大值、最小值、平均值和标准差：

# Value of the largest integer in the matrix
print("Max ", np.max(z))
# Value of the smallest integer in the matrix
print("Min ", np.min(z))
# Mean of elements in the matrix
print("Mean ", np.mean(z))
# Standard deviation
print("Standard deviation ", np.std(z))

这将产生以下输出：

('Max ', 4)
('Min ', 1)
('Mean ', 2.5)
('Standard deviation ', 1.1180339887498949)
[Finished in 0.207s]

还有许多其他操作可以在ndarray上执行，包括常见的数学函数，如 sin、cos、log 和平方根，以及统计函数，如寻找相关系数和累积和，我们很快将使用其中一些。

检查矩阵

现在，我们将处理一些有用的函数，这些函数可以帮助我们更多地了解我们正在使用的任何数组。让我们继续使用到目前为止一直在使用的相同矩阵/ndarray z：

1. 让我们打印一个矩阵的信息：

# 1\. Information about a matrix
print("Information: ")
print(np.info(z))

输出将如下所示：

Information: 
class:  ndarray
shape:  (2, 2)
strides:  (16, 8)
itemsize:  8
aligned:  True
contiguous:  True
fortran:  False
data pointer: 0x7ff57665fef0
byteorder:  little
byteswap:  False
type: int64
None

1. 现在，为了确定矩阵的形状，请写下以下代码：

# 2\. Gives the shape of the matrix
print("Shape: ")
print(np.shape(z))

输出将如下所示：

Shape: 
(2, 2)

1. 要检查矩阵是 2D 还是 3D 矩阵，请写下以下代码：

# 3\. Dimensions of the matrix
print("Dimensions: ")
print(np.ndim(z))

输出将如下所示：

Dimensions: 
2

1. 要打印矩阵的数据类型，请使用以下代码：

# 4\. Data type of the matrix
print("Data type of elements")
print(z.dtype.name)

输出将如下所示：

Data type of elements
int64

1. 要打印矩阵的长度，请使用以下代码：

print("Length of the ndarray: ")
print(len(z))

输出将如下所示：

Length of the ndarray: 
2

正如我们所看到的，info函数已经显示了我们调用的另外两个函数的值，即 shape 和 type。然而，这些函数只能起到有限的作用，有时这就是所需的。正如我们所知，多维ndarray是一个数组的数组，而 NumPy 数组的len函数将始终是第一维的长度。如果z是一个 2D 矩阵，那么len(z)将是z中的行数。

在接下来的练习中，我们将创建一个矩阵。我们可以用嵌套列表创建一个矩阵，但是这个问题将详细说明矩阵在现实世界中是如何打包和利用的。

练习 6.01：计算阳光到达地球所需的时间

在这个练习中，我们将计算阳光到达地球所需的时间。

正如我们所知，地球围绕太阳在椭圆轨道上运行。因此，地球和太阳之间的距离会发生变化，这将改变光到达地球所需的时间。我们可以使用三个主要方程来处理这个问题。

计算时间的数学公式如下：

图 6.2：计算时间的公式

我们需要计算地球和太阳之间的距离r：

图 6.3：计算距离的公式

在前述方程中，a的值为 149,600,000 千米，这是半长轴距离，e为 0.0167，这是地球轨道的离心率，θ是从近日点开始的角度。

需要计算前述方程中所需的因变量θ，计算如下：

图 6.4：计算因变量的公式

请注意这里n是从 1 月 3 日近日点开始的天数。为了简化问题，我们将其视为一年的开始。

不要被这些方程式所困扰，它们只是简单的常数数学乘法，可以轻松通过一个称为math的巧妙 Python 库解决。

现在让我们开始练习：

1. 首先，导入math和numpy库：

import math
import numpy as np

我们稍后将使用这些库。

1. 接下来，定义两个常数并使用大写字母，这是命名这些常数的标准 Python 做法：

def earth_sun_distance():
    # Semi-major axis between earth and sun
    A = 149600000
    # Eccentricity of earth
    E = 0.0167
    l = []

A在这里是地球和太阳之间的半长轴距离。

E被称为地球的离心率。

l是我们初始化用于存储值的列表。

1. 让我们跳到代码的主要部分。对于 365 天中的每一天，计算theta，因为每年的每一天它都不同。然后，计算地球到太阳的距离，最后将该距离附加到列表中：

    # Calculating the distance between earth and sun
    for i in range(365):
        theta = (2 * math.pi * i) / 365.25
        r = A*(1 - E**2) / (1 - (E * math.cos(theta)))
        l.append(r)
    return l

请注意我们之前导入的math库中的math.pi和math.cos函数的使用。

1. 计算所需的时间（以秒为单位），假设光速为恒定值 299,792 千米/秒：

# Calculating the time taken
S = 299792
t = np.divide(l, S)

在这里，我们首先利用 NumPy 数组的功能，使用divide函数，它将值应用于列表的所有成员，而无需使用循环。我们将其值存储在t中，它会自动转换为 NumPy 数组。

1. 最后，我们在这里做了两件事。首先，我们使用另一个有用的 Python 函数zip()，它将两个列表的相应元素绑在一起，然后我们使用np.asarray()函数，将元组列表转换为 NumPy 数组：

sunny = np.asarray(list(zip(l, t)))
print("Earth sun distance: \n", sunny)

运行程序以查看输出：

Earth sun distance:
[[  1.52098320e+08   5.07346160e+02]
 [  1.52097938e+08   5.07344885e+02]
 [  1.52096791e+08   5.07341061e+02]
 [  1.52094881e+08   5.07334688e+02]
 [  1.52092207e+08   5.07325770e+02]
 [  1.52088771e+08   5.07314309e+02]
 ...
 [  1.52072354e+08   5.07259546e+02]
 [  1.52078259e+08   5.07279242e+02]
 [  1.52083406e+08   5.07296411e+02]
 [  1.52087793e+08   5.07311046e+02]
 [  1.52091420e+08   5.07323143e+02]
 [  1.52094284e+08   5.07332697e+02]
 [  1.52096385e+08   5.07339707e+02]
 [  1.52097723e+08   5.07344168e+02]]
[Finished in 0.197s]

现在，我们以系统化的表格格式拥有了地球和太阳之间的距离值以及光线到达地球所需的时间。我们可以继续向我们的矩阵添加其他参数，这就是使用矩阵和 NumPy 数组的灵活性。

请注意，这些值绝不准确，我们为简单起见做了一些安全假设，但这仍然是矩阵可以用于几乎任何事情的很好的例证。还要注意，这里反映的值是 Python 中使用的科学计数法格式，并且可以根据需要轻松转换为浮点数或其他类型。左边的值是以千米为单位，右边的值是以 507.346...秒的形式。

1. 将结果附加如下：

d = []
for i in range(1,len(l) - 1):
    d.append(l[i]-l[i-1])
print(d)

输出的一部分如下：

[-382.2014582455158, -1146.4797523021698, -1910.3842301666737,
 -2673.6658524870872, -3436.075836390257, -4197.365758448839,
 -4957.287656396627, -5715.5941315591335, -6472.038449823856,
 -7226.374643236399, -7978.357610076666, -8727.743215203285,
 -9474.288]

注意

要访问此特定部分的源代码，请参阅packt.live/3irS3Bk.

您也可以在packt.live/3abV9pe上在线运行此示例。

矩阵中的运算和乘法

现在我们了解了如何执行简单的操作，让我们执行一系列操作，例如调整大小、重塑和转置，形成一个新的矩阵。

当矩阵中的行和列的索引交换时，沿对角线翻转它们，这被称为矩阵的转置。现在让我们来看看如何转置矩阵。这可以通过三种不同的方式完成，如下所示：

print("matrix z: ")
print(z)
# Transpose matrix
# Method 1
print("new matrix: ")
print(np.transpose(z))
# Method 2
print(z.transpose())
# Method 3
t = z.transpose()
print(t)

如果您运行此代码，输出将如下所示：

matrix z: 
[[1 2]
 [3 4]]
new matrix: 
[[1 3]
 [2 4]]
[[1 3]
 [2 4]]
[[1 3]
 [2 4]]
[Finished in 0.207s]

在方法 3 中，我们将转置矩阵的值分配给一个新变量。

我们现在将看到的函数是在执行矩阵操作时最常用的函数之一。

我们将处理的第一个函数是展平。将矩阵转换为单行的过程称为矩阵的展平：

# Flatten the array
y = z.ravel()
print(y)

这将产生以下输出：

[1 2 3 4]

现在让我们看一下各种比较运算符：

# Comparison operators on matrix
print(z == 3)

在这种情况下，矩阵中的所有值都与基本值（在本例中为3）进行比较，并且布尔结果针对矩阵中的相应索引显示。输出如下：

[[False False]
 [ True False]]

要检查z的值是否小于3，请使用以下代码：

print(z < 3 )

输出如下：

[[ True  True]
 [False False]]

reshape是一个函数，用于根据函数内部传递的行和列的值来改变矩阵的维度。要重塑矩阵，请使用以下代码：

# Reshaping the matrix
r = z.reshape(4,1)
print(r)

这将产生以下输出：

[[1]
 [2]
 [3]
 [4]]

要调整矩阵的大小，请使用以下代码：

# Resizing the matrix
resize = np.resize(z,(3,3))
print(resize)

这将产生以下输出：

 [[1 2 3]
 [4 1 2]
 [3 4 1]]
[Finished in 0.205s]

请注意，当我们使用resize函数时，值会被迭代重复，直到满足大小，即使可能并未添加原始矩阵中的所有值。还要注意，通常使用reshape函数代替ravel函数来展平矩阵。

矩阵中的轴

这个相对简单的主题易于理解，但同样容易误解，因此我们需要独立处理它。对于 Python 中的数组，为具有多个维度的任何矩阵或数组定义轴。在处理复杂的数据科学和数据操作问题时，我们经常需要处理超过两个维度的数据，这很难可视化并且容易混淆。为了简化这一点，矩阵中的维度由轴表示。

简单地说，2D 矩阵将有两个轴，水平和垂直，但在这种情况下，它们将以数字形式表示或命名。第一个轴称为轴 0，沿着行向下垂直运行，第二个称为轴 1，沿着列水平运行。

我们之前使用的一组函数可以用于沿着单个轴运行，这在大型数据集的情况下减少了计算的开销。让我们处理一些例子。为了清晰起见，我们将创建一个稍大一些的矩阵：

import numpy as np
z = np.array([[1, 5, 9, 4],\
              [8, 3, 7, 6]])
print(z.max(axis = 0))
print(z.max(axis = 1))

这将产生以下输出：

[8 5 9 6]
[9 8]
[Finished in 0.198s]

这里发生的是沿着每个轴计算最大值。在返回的第一个数组中，比较是1和8，5和3，9和7，以及4和6，因为这些是沿着轴 0 的唯一两个元素。类似地，在轴 1 的情况下，比较是沿着子数组的四个元素进行的，并返回最大元素。

让我们再举一个例子：

print(z.sum(axis = 1))

你能猜到结果吗？让我们看一下输出：

[19 24]
[Finished in 0.255s]

现在让我们看一个最后的、更复杂的例子：

print(np.concatenate((z[1], z[0]), axis=0))

这将产生以下输出：

[8 3 7 6 1 5 9 4]
[Finished in 0.252s]

我们所做的首先是使用一个接受两个数组的连接函数。取出的两个数组分别是数组z的第一个和第二个元素，分别是[8 3 7 6]和[1 5 9 4]。由于这两个数组各自具有单个维度，我们沿着轴 0 取它们。如果我们在这里输入轴 1，NumPy 将抛出AxisError，如下所示：

print(np.concatenate((z[1], z[0]), axis=1))

这将产生以下输出：

Traceback (most recent call last):
  File "/matrix.py", line 9, in <module>
    print(np.concatenate((z[1], z[0]), axis=1))
numpy.core._internal.AxisError: axis 1 is out of bounds for array of dimension 1

练习 6.02：矩阵搜索

在这个练习中，我们将在按升序排序的矩阵中搜索给定的输入值，无论是按行还是按列。这将帮助我们理解在矩阵内部遍历的一般规则，特别是如果我们不使用 NumPy 数组。

为了给出一个提示，我们将在矩阵上实现二分搜索。即使您以前没有处理过二分搜索，这也足够容易理解。

目标是根据值是否存在于矩阵中返回True或False值：

1. 让我们定义我们要搜索的输入矩阵：

matrix = [[7, 10, 15, 18],\
          [25, 29, 35, 47],\
           [56, 78, 85, 104]]

1. 现在，让我们定义并编写一个名为matrixsearch()的函数，它将以此矩阵作为输入，以及我们要搜索的值。我们首先将处理边缘情况，即矩阵为空或目标值为非零的情况：

def matrixsearch(matrix, value):
    # Check for edge cases
    if value is None or not matrix:
         return False

1. 接下来，我们将定义四个变量：

# Initialize the variables
    row = len(matrix)
    col = len(matrix[0])
    start = 0
    end   = row * col - 1

请注意这里row和column变量的初始化方式。在任何矩阵中，它们的初始化方式都是一样的，值得理解。start和end变量被初始化为矩阵中的第一个和最后一个值，因为矩阵已经排序，可以被视为一个单一的列表，从起始到对角线的另一端。

1. 现在是程序的实际逻辑，我们将把它分解成几个步骤来帮助理解。在从起始到结束的循环中，首先我们找到矩阵的中点（将矩阵视为一个列表）：

    while start <= end:
        mid = int((start + end) / 2)

1. 然后，我们定义一个名为pointer的变量，它由我们找到的这个中间值的值初始化：

        pointer = matrix[int(mid/col)][int(mid%col)]
        print(int(mid/col), int(mid%col))

请注意，/用于除法，%在这里用作模数。因此，在第一次迭代中，它们的值将分别为（1,1）。

1. 现在，我们来到二分查找的核心部分，我们通过与我们拥有的值进行比较来增加或减少我们的指针。如果我们找到该值，我们返回True，否则我们将继续循环，直到我们找到或者在最后找不到任何东西时返回False：

        if pointer == value:
            return True
        elif pointer < value:
            start = mid + 1
        else:
            end = mid - 1
    return False
sol = matrixsearch(matrix, 78)
print(sol)

输出将如下所示：

1 1
2 0
2 2
2 1
True

在这个练习中，我们使用 NumPy 实现了对矩阵的二分查找，并且根据矩阵的值，我们的代码返回了True。

注意

要访问此特定部分的源代码，请参阅packt.live/3eVd0Ch。

您也可以在packt.live/2ZusZkj上在线运行此示例。

多个矩阵

到目前为止，我们已经学会了在拥有单个矩阵时执行操作和操作。接下来，我们将处理多个矩阵。Python 中矩阵的组合在当今的数据科学中最常用，因为它需要存储和处理大型数组。让我们从一个简单的例子开始。我们将取两个矩阵z和x，并将值相乘如下：

import numpy as np
z = np.array([[1, 2],\
              [3, 4]])
x = np.array([[4, 5],\
              [7, 8]])
print(np.multiply(x,z))
print(np.multiply(z,x))

如果运行上述代码，您将看到以下输出：

[[ 4 10]
 [21 32]]
[[ 4 10]
 [21 32]]
[Finished in 0.206s]

输出显示，直观地，两个矩阵的相应元素相乘得到一个乘积值。这只是元素级的乘法，或者在数学上称为 Hadamard 乘积。

现在让我们稍微改变矩阵x：

import numpy as np
z = np.array([[1, 2],\
              [3, 4]])
x = np.array([[4, 5, 6],\
              [7, 8, 9]])
print("Multiplying with a number: ")
print(np.multiply(3,x))
print(np.multiply(x,3))
print("Multiplication between matrices of different sizes: ")
print(np.multiply(x,z))
print(np.multiply(z,x))

现在，输出将如下所示：

Multiplying with a number: 
[[12 15 18]
 [21 24 27]]
[[12 15 18]
 [21 24 27]]
Multiplication between matrices of different sizes: 
Traceback (most recent call last):
    File "/Users/…/matrix operations.py", line 52, in <module>
    print(np.multiply(x,z))
ValueError: operands could not be broadcast together with shapes (2,3) (2,2) 
[Finished in 0.204s]

我们得到的是ValueError，这是由于 NumPy 中数组的广播属性导致的。

广播

重要的是要理解广播的概念，以便我们知道在使用数组进行矩阵运算时允许什么，不允许什么。

简而言之，广播是 NumPy 处理两个具有不同形状的数组的方式。一般规则是，在两者中较小的数组将以某种方式在较大的数组上进行广播，使它们兼容。根据官方文档，广播的一般规则如下：

• 它从尾部维度开始向前工作。

• 当它们中的一个为 1 时，或者当它们两个相等时，比较的两个维度是相等的。

注意

您也可以参考官方文档docs.scipy.org/doc/numpy/user/basics.broadcasting.html。

因此，正如我们在前面的例子中看到的，当乘以相同维度的矩阵和标量变量时，乘法运算完美地进行。另一方面，如果两个矩阵的维度不同，将抛出ValueError，因为 NumPy 无法有效地将较小矩阵的值广播到较大矩阵上。这种广播主要是在内部完成的，以使数组更快速和更节省内存。它提供了一种将数组向量化以在 C 中实现循环而不是 Python 的方法，这有效地使其更快。在这里需要记住的一件重要的事情是，在一对 NumPy 数组的情况下，操作是基于逐元素的。为了克服维度的问题，采用的两种主要方法是reshape和newaxis。在我们结束之前，让我们看一下广播概念的另一个变化：

import numpy as np
z = np.array([[1,2],\
              [3]])
x = np.array([[4,5]])
print(np.multiply(x,z))

猜猜输出会是什么样子？让我们来看看：

[[list([1, 2, 1, 2, 1, 2, 1, 2]) list([3, 3, 3, 3, 3])]]
[Finished in 0.244s]

由于这里的数组z不是一个常规的方形数组，NumPy 在内部不将其解释为矩阵，而是将其视为常规对象行，并对其进行逐元素乘法。因此，z[0]乘以x[0]，z[1]乘以x[1]，以产生在这种情况下恰好是一个列表的对象。

多个矩阵的操作

我们现在将执行两个或多个矩阵之间的操作，并查看将帮助我们实现这一目标的函数。我们将涵盖如何编写矩阵的逆、逻辑运算符、点积、特征值和特征向量、外积以及矩阵的行列式。

应该注意到，矩阵还有许多其他用途，官方的 NumPy 文档是一个非常好的资源，可以根据用户的需求引用信息。我们将要涵盖的大部分主题都是 NumPy 库的线性代数包的一部分。对于我们将要涵盖的每个主题，物理学和数学中都有更广泛的应用，但应该足够了解它们在理解矩阵方面起着非常重要的作用。

注意

有关 NumPy 库的更多信息，请参阅docs.scipy.org/doc/numpy/reference/。

单位矩阵

单位矩阵在对角线上有 1，在其他地方都是 0。我们将使用 NumPy 的linalg函数创建单位矩阵：

import numpy as np
from numpy.linalg import inv
def identity():
    print(np.identity(3))

输出将如下所示：

[[ 1\.  0\.  0.]
 [ 0\.  1\.  0.]
 [ 0\.  0\.  1.]]
[Finished in 0.206s]

单位矩阵

eye函数类似于单位矩阵，唯一的区别是可以偏移矩阵的值。这意味着它将从第 k 行开始创建一个单位矩阵，如下所示：

def eye():
    print(np.eye(3, k = 1))

输出将如下所示：

[[ 0\.  1\.  0.]
 [ 0\.  0\.  1.]
 [ 0\.  0\.  0.]]
[Finished in 0.277s]

矩阵的逆

逆或乘法逆是当你将其与原始矩阵相乘时产生单位矩阵的矩阵。矩阵的逆在应用于 3D 几何和图形时最常用：

def inverse():
    z = np.array([[1,2],\
                  [3,4]])
    z_inv = inv(z)
    product = np.dot(z, z_inv)
    print(z_inv)
    print(product)

这将产生以下输出：

# Output of print(z_inv)
[[-2\.   1\. ]
 [ 1.5 -0.5]]
# Output of print(product)
 [[  1.00000000e+00   0.00000000e+00]
 [  8.88178420e-16   1.00000000e+00]]
[Finished in 0.202s]

我们这里有两个输出。第一个是所谓的矩阵的逆，第二个是我们将逆乘以原始矩阵，使用dot函数产生单位矩阵。显示的值是浮点数，不应该成为一个问题。

逻辑运算符

我们将在这里创建两个列表，包含二进制的True（1）或False（0）值。我们将使用 NumPy 的内置logical_and()函数对它们进行AND操作的输出：

def and_op():
    m1 = [True, False, False]
    m2 = [True, False, True]
    print(np.logical_and(m1, m2))

输出将如下所示：

[ True False False]
[Finished in 0.253s]

非常简单。你也可以使用 1 和 0 代替True和False，它仍然会给出相同的结果。实际上，只要不是 0，它就被视为True。让我们看一个使用 1 和 0 的例子：

def and_op():
    m1 = [0, 1, 0]
    m2 = [1, 1, 0]
    print(np.logical_and(m1, m2))

输出将如下所示：

[False  True False]
[Finished in 0.192s]

使用logical_or、logical_not和logical_xor函数也可以对其他逻辑函数进行同样的操作。

外函数或向量积

np.outer是用于生成向量或两个矩阵的叉积的函数：

def vector():
    horizontal = np.array([[1,3,2]])
    vertical = np.array([[2],\
                         [0],\
                         [1]])
    print("Output for dimension 1 x 1: ")
    print(horizontal.dot(vertical))
    print("Output for dimension 3 x 3: ")
    print(vertical.dot(horizontal))
    print("Output using outer for getting cross product: ")
    print(np.outer(vertical.ravel(), horizontal.ravel()))
    print(np.outer(horizontal.ravel(), vertical.ravel()))

这将产生以下输出：

Output for dimension 1 x 1: 
[[4]]
Output for dimension 3 x 3: 
[[2 6 4]
 [0 0 0]
 [1 3 2]]
Output using outer for getting cross product: 
[[2 6 4]
 [0 0 0]
 [1 3 2]]
[[2 0 1]
 [6 0 3]
 [4 0 2]]
[Finished in 0.289s]

到目前为止，我们已经学会了使用矩阵的各种方法。我们使用的方法列表绝对不是限制性的，随时需要进行某种操作时，深入探索库始终是一个好习惯。值得再次提到的是，根据工作领域的要求，有几种特定类型的矩阵具有有限的使用情况。

使用矩阵解决线性方程

线性方程是代数的基本组成部分，任何学过基本初等数学的人都知道它们是如何工作的。让我们简要地介绍一下它们，然后我们可以看看如何在 Python 中使用矩阵轻松解决它们。

线性方程通常以以下形式出现：

图 6.5：计算线性方程的公式

这里，a1、a2、..、an 是系数，x1、x2、.. xn 是变量。

这些具有两个变量的线性方程可以在二维空间图中表示，其中x是水平维度，y是垂直维度。

让我们快速举一个具有两个变量的线性方程的例子。假设方程是y = 2x + 6。这种表示被称为斜率-截距形式，格式为y = mx + c，其中m是斜率，c是方程的y截距。

这里，m=2，c=6，并且可以通过绘制不同的x和y值在图表上画出一条线：

图 6.6：在二维空间中表示 y = 2x + 6

不详细讨论，我们可以想象平面上可能有另一条线，它要么与该线平行，要么与该线相交。线性方程的目的是找到这些线的交点，并根据交点的值找到变量x和y的值。随着维度的增加，可视化变得困难，但基本上概念是相同的。矩阵极大地简化了解这些方程的过程。通常有两个矩阵，一个包含x和y的系数，另一个包含变量x和y。它们的点积产生了结果矩阵，即先前提到的常数或y截距。一旦我们快速进行了一次练习，就很容易理解。

练习 6.03：使用矩阵执行线性方程

我们现在将使用矩阵来解决一个线性方程问题。

约翰出差三天，心情不错，准备花光所有的钱。他随身携带了三种货币面额。第一天，约翰在他喜欢的最新电子平板上花了 435 美元，其中a面额的37，b面额的20，c面额的12。第二天，他去跳伞，分别花了a、b和c面额的15、32和4，总共 178 美元。第三天，他决定去看戏，花了 70 美元，分别是a、b和c面额的5、40和2。你能告诉各个面额的值是多少吗？

从问题中我们可以看出，有三个方程和三个未知变量需要计算。

1. 让我们使用 NumPy 数组将三天的已知值放入矩阵中：

import numpy as np
# Input
z = np.array([[37, 20, 12],\
              [15, 32, 4],\
              [5,  40, 2]])

我们现在有了需要处理的矩阵。有几种方法可以解决这个问题。本质上，这就是我们需要做的：

Ax = b

其中A是我们知道值的矩阵，x是具有未知变量的矩阵，b是结果矩阵。

1. 结果的b矩阵如下。这些是约翰在三天内的花费金额：

r = np.array([[435],[178],[70]])

有几种方法可以在 Python 中解决这个问题：

方法 1：通过计算 x = A-1b 来找到x：

1. 让我们首先使用我们之前学到的函数来计算矩阵 A 的逆：

print(np.linalg.inv(z))

注意

我们将使用矩阵的点积而不是纯乘法，因为这些不是标量变量。

输出如下：

[[-0.06282723  0.28795812 -0.19895288]
 [-0.0065445   0.0091623   0.02094241]
 [ 0.28795812 -0.90314136  0.57853403]]

这里不需要理解这个矩阵，因为这只是一个中间步骤。

1. 然后，我们取两个矩阵的点积来产生一个矩阵X，这是我们的输出：

X = np.linalg.inv(z).dot(r)
print(X)

输出如下：

[[10\.  ]
 [ 0.25]
 [ 5\.  ]]

方法 2：使用linalg包中的内置方法：

1. 甚至可以使用另一个名为solve()的 NumPy 函数更轻松地完成相同的事情。让我们在这里将输出变量命名为y。

y = np.linalg.solve(z,r)
print(y)

这产生了以下输出：

[[10\.  ]
 [ 0.25]
 [ 5\.  ]]

而且在一行代码中，我们就能够用 Python 解决线性方程。我们可以推断和理解，使用 Python 可以轻松解决具有大量未知变量的类似方程。

因此，我们可以看到使用方法 1 和方法 2 后的输出是相同的，即$10、25 美分和$5，这些是我们试图确定的相应面额。

如果我们接收约翰的开销信息是迭代的而不是一次性的呢？

1. 让我们首先添加我们收到的有关约翰开销的信息：

a = np.array([[37, 20, 12]])

1. 然后，让我们也添加有关约翰其他两笔开销的信息：

b = np.array([[15, 32, 4]])
c = np.array([[5,  40, 2]])

1. 我们可以轻松地使用concat()函数将这些数组绑定在一起形成一个矩阵：

u = np.concatenate((a, b, c), axis=0)
print(u)

这产生了以下输出：

[[37 20 12]
 [15 32  4]
 [ 5 40  2]]
[Finished in 0.188s]

这是我们用于前面程序的相同输入。

同样，如果我们有更多这样的数据，我们可以应用循环来形成一个更大的矩阵，然后用它来解决问题。

注意

要访问此特定部分的源代码，请参考packt.live/3eStF9N.

您也可以在packt.live/38rZ6Fl上在线运行此示例。

转移矩阵和马尔可夫链

现在，我们将看一下矩阵的一个应用，这是一个独立的研究领域。马尔可夫链利用转移矩阵、概率和极限来解决现实世界的问题。现实世界很少像我们为解决它们创建的数学模型那样完美。一辆汽车可能想要从 A 点到 B 点，但在现实中，距离和速度可能是不够的参数。一只过马路的猫可能会完全改变计算出的汽车行驶时间。股市可能在几天内似乎遵循可预测的模式，但在夜间发生了一个事件，完全崩溃了。这个事件可能是一些全球性事件、政治声明或公司报告的发布。当然，我们在数学和计算模型方面的发展还没有达到我们可以预测每一个这些事件结果的地步，但我们可以尝试并确定某些事件发生的概率比其他事件更高。以前面的例子为例，如果公司报告将在特定日期发布，那么我们可以预期某只股票会受到影响，并且我们可以根据对公司进行的市场分析来对此进行建模。

马尔可夫链就是这样一种模型，其中依赖于马尔可夫性质的变量只考虑当前状态来预测下一个状态的结果。因此，从本质上讲，马尔可夫链是一个无记忆过程。它们利用转换状态图和转换矩阵进行表示，用于映射在给定当前事件的情况下发生某一事件的概率。当我们称之为无记忆过程时，很容易将其与与过去事件无关的东西混淆，但实际情况并非如此。当我们举例说明它们是如何工作时，事情就容易理解得多了。在我们开始使用马尔可夫链之前，让我们首先深入了解转换状态和矩阵的工作原理，以及它们为什么被使用。

马尔可夫链的基础知识

为了保持简单，让我们将概念分解成片段，并从我们已有的信息中逐步学习它们，然后才能将它们放在一起理解这些概念。

随机模型与确定性模型

当我们试图解决现实世界中的问题时，我们经常遇到超出我们控制范围的情况，很难制定。模型旨在模拟给定系统的功能方式。虽然我们可以在我们的模型中考虑系统的大部分元素，但许多方面无法被确定，因此根据它们发生的可能性进行模拟。这就是概率介入的地方。我们试图找到在一组情况下发生特定事件的概率。我们使用的模型有两种主要类型，确定性和随机性。确定性模型是具有一组参数值和函数的模型，并且可以形成可预测的数学方程，并提供可预测的输出。随机模型包括随机性，即使它们具有初始值和方程，它们也提供可能性的定量值。在随机模型中，不会有固定的答案，而是某些事件发生的可能性大于其他事件。线性规划是确定性模型的一个很好的例子，而天气预测是一个随机模型。

过渡状态图

广义上，这些构成了面向对象编程的基础，您可以根据给定的事件和条件描述对象可能具有的所有可能状态。状态是在满足某个先前条件时对象在某一特定时刻的状态。让我们举个例子来说明这一点：

图 6.7：风扇的状态转移图

这是台式风扇调节器的状态转移图，通常情况下，每次我们顺时针旋转它，它都会改变状态，直到它转回“关闭”位置。

在这里，台式风扇的状态在速度方面发生变化，而动作是扭转。在这种情况下，它是基于事件的，而在其他一些情况下，它将基于满足的条件。

让我们以使用与我们将要实施的内容一致的马尔可夫链进行文本生成的示例。我们将回顾一个儿歌的前两行：

汉普蒂·韦尔壁上坐着，

汉普蒂·韦尔摔了个大跟头。

首先，让我们准备一个包含句子中所有单词频率的频率表：

图 6.8：韵文中单词的频率表

标记是存在的单词数，而键是唯一单词。因此，值将如下所示：

标记 = 12

键=9

我们甚至可能不需要我们在这里学到的一切，但一旦您决定实施更复杂的问题，这将是重要的。每个转换图都有一个起始状态和结束状态，因此我们将在这里添加这两个状态作为关键：

图 6.9：起始和结束状态的频率表

然后我们准备一个状态图表来显示从一个状态到下一个状态的转换。在这种情况下，需要显示哪个词将跟随当前词。因此，我们将形成如下的成对关系：

图 6.10：词对

如果我们按照关键词而不是标记来压缩这个，我们会看到一些关键词有多个转换，如下所示：

图 6.11：一些关键词有多个转换

这不仅是为了减少状态转换，还为我们正在做的事情增加意义，我们很快就会看到。整个目的是确定哪些词可以与其他词配对。我们现在准备绘制我们的状态转换图。

我们将所有唯一的关键字添加为状态，并显示这些单词可以转换到哪些状态：

图 6.12：状态转换图

如果我们看前面的图表，我们可以根据给定的条件跟随任何单词来完成韵脚。剩下的是关键词在给定单词之后出现的概率。为此，请看下面的图表，我们可以很直观地看到概率如何根据它们的频率在关键词之间分配。请注意，Humpty 总是跟在 Dumpty 后面，因此概率为 1：

图 6.13：带概率的状态转换图

现在我们已经讨论了状态转换图，我们将继续绘制转换矩阵。

转换矩阵

在马尔可夫过程中，我们需要以数学格式显示状态转换的概率，为此使用转换矩阵。行和列简单地是转换图的状态。转换矩阵中的每个值显示从一个状态到另一个状态的转换概率。可以想象，这种矩阵中的许多值将为 0。对于前面讨论的问题，转换矩阵将有 9 个状态和许多 0。我们将以一个更简单的转换图为例，并找到其相应的矩阵：

图 6.14：具有状态 1、2 和 3 的状态图

当我们看这个图表时，我们看到了三个转换状态。请注意，我们在这里没有明确包括起始和结束状态，但在某些情况下可能是必要的。外向箭头表示从一个状态到下一个状态的转换。一旦我们有了图表，就很容易绘制矩阵。写行和列等于图表的状态。在这种情况下，将是 3。然后，第 0 行将显示第 1 个状态的转换，第 1 行将显示第二个状态，依此类推。一般来说，矩阵中的每一行表示一个状态的转换概率。

让我们先看矩阵，然后我们可以讨论一些其他事情：

图 6.15：转换矩阵

除了行的属性，我们还可以观察到另一件事。给定行中所有概率的总和将始终等于 1。在这里，第一行的总和将是 1/5 + 2/5 + 2/5 = 5/5 = 1。这是因为这些状态是穷尽的。

如果两个给定状态之间没有转换，那么矩阵中的状态值将为 0。我们可以通过比较矩阵中的值的数量和图表中可以看到的状态转换的数量来验证这一点。在这种情况下，它们都等于 7。

练习 6.04：查找状态转换的概率

给定一个包含四个状态A、B、C和D的数组，这些状态是随机生成的，让我们找到这四个状态之间的转移概率。我们将找到每个状态转移的概率，并从中形成一个转移矩阵。

让我们从给定的输入数组在 Python 中生成一个转移矩阵。我们将在未来推断相同的概念，用于创建马尔可夫链。

1. 使用 Python 中的random包生成一个包含字符A、B、C和D的随机状态数组。然后，我们将通过创建一个名为LEN_STR的常量来定义我们想要的元素数量，在这种情况下，我们将其设置为50：

# Generate random letters from 4 states A B C D
import random
tokens = []
LEN_STR = 50
for i in range(LEN_STR):
    tokens.append(random.choice("ABCD"))
print(tokens)
LEN_TOKENS = len("ABCD")

这产生了以下输出：

['C', 'A', 'A', 'B', 'A', 'A', 'D', 'C', 'B', 'A', 'B',  'A', 'A', 'D', 'A', 'A', 'C', 'B', 'C', 'D', 'D', 'C',  'C', 'B', 'A', 'D', 'D', 'C', 'A', 'A', 'D', 'C', 'A',  'D', 'A', 'A', 'A', 'C', 'B', 'D', 'D', 'C', 'A', 'A',  'B', 'A', 'C', 'A', 'D', 'D']

注意

我们从字符串的长度生成的另一个常量LEN_TOKENS的使用将指示问题中将存在的状态的数量。

1. 接下来，我们将找到字母的相对值并将它们转换为整数，主要是因为整数更容易进行计算：

# Finding the relative values with ordinal values of 
# ASCII characters
relative_value = [(ord(t) - ord('A')) for t in tokens]
print(relative_value)

这产生了以下输出：

[2, 0, 0, 1, 0, 0, 3, 2, 1, 0, 1, 0, 0, 3, 0, 0, 2, 1, 
 2, 3, 3, 2, 2, 1, 0, 3, 3, 2, 0, 0, 3, 2, 0, 3, 0, 0, 
 0, 2, 1, 3, 3, 2, 0, 0, 1, 0, 2, 0, 3, 3]

这里我们使用基数值是为了方便，但我们也可以使用字典或其他方法来做到这一点。如果你不知道，ord()函数在这里返回字符串中字符的 ASCII 值。例如，A和D的 ASCII 值分别为65和68。

1. 现在，找到这些 ASCII 值之间的差异，并将它们放入一个名为ti的列表中。我们也可以直接更新令牌列表，但为了清晰起见，我们将它们分开保持在不同的列表中：

#create Matrix of zeros
m = [[0]*LEN_TOKENS for j in range(LEN_TOKENS)]
print(m)
# Building the frequency table(matrix) from the given data
for (i,j) in zip(relative_value, relative_value [1:]):
    m[i][j] += 1
print(list(zip(relative_value, relative_value [1:])))
print(m)

这产生了以下输出：

[[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]
[(2, 0), (0, 0), (0, 1), (1, 0), (0, 0), (0, 3), (3, 2),  (2, 1), (1, 0), (0, 1), (1, 0), (0, 0), (0, 3), (3, 0),  (0, 0), (0, 2), (2, 1), (1, 2), (2, 3), (3, 3), (3, 2),  (2, 2), (2, 1), (1, 0), (0, 3), (3, 3), (3, 2), (2, 0),  (0, 0), (0, 3), (3, 2), (2, 0), (0, 3), (3, 0), (0, 0),  (0, 0), (0, 2), (2, 1), (1, 3), (3, 3), (3, 2), (2, 0),  (0, 0), (0, 1), (1, 0), (0, 2), (2, 0), (0, 3), (3, 3)]
[[8, 3, 3, 6], [5, 0, 1, 1], [5, 4, 1, 1], [2, 0, 5, 4]]

我们现在根据我们之前生成的LEN_TOKENS常量的大小初始化了一个零矩阵，并用它来构建一个零矩阵。

在第二部分中，我们正在创建成对的元组，就像我们在之前的问题中所做的那样，并根据两个状态之间的转换次数更新转移矩阵的频率。这个的输出是前面代码中的最后一行。

注意

我们选择迭代地更新矩阵m的值，而不是创建新的矩阵。

1. 我们现在将生成概率，这仅仅是给定行中的相对频率。就像在第一行中，从 A 到 A 的转换是 8，从 A 到任何状态的总转换次数是 20。所以，在这种情况下，概率将是8/20 = 0.4：

# Finding the Probability
for state in m:
    total = sum(state)
    if total > 0:
        state[:] = [float(f)/sum(state) for f in state]

代码对每一行都是这样的，如果sum函数大于0，我们就找到概率。请注意这里使用float函数是为了避免在某些 Python 版本中进行类型转换为int。另外，请注意使用state[:]，它创建了一个浅拷贝，从而防止内部类型转换冲突。

1. 现在，让我们通过添加以下代码来打印state对象：

for state in m:
      print(state)

在这里，我们遍历矩阵中的行并打印出值，这就是我们的转移矩阵。

这产生了以下输出：

[0.4, 0.15, 0.15, 0.3]
[0.7142857142857143, 0.0, 0.14285714285714285,  0.14285714285714285]
[0.45454545454545453, 0.36363636363636365,  0.09090909090909091, 0.09090909090909091]
[0.18181818181818182, 0.0, 0.45454545454545453,  0.36363636363636365]

因此，我们能够构建一个描述状态转移的转移矩阵，它显示了从一个状态转移到下一个状态的概率。因此，A找到A作为下一个字母的可能性是0.4，A转到B的概率是0.15，依此类推。

注意

要访问此特定部分的源代码，请参阅packt.live/31Ejr9c.

您也可以在packt.live/3imNsAb.上在线运行此示例。

马尔可夫链和马尔可夫性质

转换状态和矩阵基本上涵盖了大部分马尔可夫链。此外，还有一些值得理解的东西。如前所述，当变量仅依赖于当前状态来确定其下一个状态时，马尔可夫性质适用。形成的概率模型可能确定当前状态的结果的可能性，但过去的状态被视为独立的，不会影响结果。以抛硬币为例；我们可以创建一个关于正面或反面概率的图表，但这不会决定结果。

马尔可夫性质本质上应该满足两个标准：

• 它应该只依赖于当前状态。

• 它应该是特定的离散时间。

在不至于混淆的情况下，模型中考虑的时间要么是离散的，要么是连续的。抛硬币可以被视为离散时间事件，因为它有明确的结果，比如正面或反面。另一方面，天气模式或股票市场是连续时间事件；例如，天气在一天中是变化的，没有开始和结束时间来衡量其变化。为了处理这样的连续事件，我们需要使用分箱等技术使它们变得离散。简而言之，分箱意味着根据数量或时间将数据分组成固定数量。由于马尔可夫链是无记忆的，它本质上成为了一个离散时间和状态空间过程。

有一些专门用于特定目的的特殊矩阵。例如，稀疏矩阵在数据科学中被广泛使用，因为它们在内存和计算上都很高效。我们没有过多地处理矩阵内的元素操作，因为这本质上就像处理一个列表的列表，但在未来花一些时间进行这方面的工作是值得的。

除了马尔可夫链，还有一些随机过程的模型。这些包括自回归模型、泊松模型、高斯模型、移动平均模型等。每种模型都以不同的方式处理随机性，Python 中几乎所有模型都有支持库。即使在马尔可夫链中，还涉及到涉及多维矩阵或二阶矩阵、隐马尔可夫模型、MCMC 或马尔可夫链蒙特卡洛方法等复杂的主题。你必须自己决定你想深入研究的程度。

活动 6.01：使用马尔可夫链构建文本预测器

这项活动的目标是基于我们所学的知识构建我们自己的文本预测器。我们将获取一位著名领导人的演讲文本，并基于演讲内容使用马尔可夫链模型和状态转换构建文本预测器。执行以下步骤来实现这个目标：

1. 首先，找到一位著名人物的演讲的合适且足够大的文本，比如科学家或政治或精神领袖。为了开始，样本文本文件churchill.txt已经添加到 GitHub 仓库中，网址为packt.live/38rZy6v。

2. 生成一个描述状态转换的列表，显示给定单词与其后跟随的单词之间的相关性。

3. 整理你制作的列表，并通过将跟随给定单词的单词分组在不同位置形成哈希表。例如，这些后续单词将分组形成John: [cannot, completely, thought, …]：

John cannot…, John completely…, and John thought ..,

1. 使用随机单词生成器生成并为第一个单词分配一个值。

2. 最后，创建一个随机生成器，根据我们声明的转换状态创建一个句子。

提示

这项活动需要一些你应该熟悉的 Python 方法。

为了开始，你可以使用open('churchill.txt').read()从文本文件中读取演讲文本，然后使用split()将其分割成单词列表。

然后，您可以遍历列表，并将元素附加到新列表中，该列表将存储关键字和其后的单词。

然后，使用字典为新列表中的每个元组形成键值对。

可以使用np.random()函数生成随机词语语料库。

句子的形成来自于连接我们生成的列表的元素。

注意

此活动的解决方案可在第 677 页找到。

通过这种方式，我们制作了自己的文本预测器。这样的文本预测器可以被认为是在文本生成领域的广阔和快速增长的领域中的基础步骤。它远非完美；这有几个原因，其中一些如下：

• 我们选择的文本样本通常比我们选择的文本大得多。我们的文本包含约 22,000 个单词，而在实践中，数百万个单词被作为数据输入。

• 在停用词、标点和句子的开头/结尾方面，使用 NLP 的适当规则进行了更好的调节。

• 我们使用简单的随机生成器来选择我们的单词，而实际模型使用概率和统计模型来生成更准确的输出。

话虽如此，我们已经完成了我们的第一个文本预测器，更复杂的文本预测器基本上是基于我们所描述的方式。

尽管这绝对不能被认为是流畅的，但我们仍然只用了几行代码编写了我们的第一个文本预测器，这是一个很好的开始。

摘要

在本章中，我们能够涵盖矩阵的主题，这对于许多主题都是基础，无论是在数学中还是在使用 Python 中。今天的数据科学主要基于矩阵的有效使用。我们研究了它们在马尔可夫链的形式中的应用，尽管这是一个重要的主题，但在应用矩阵的领域下，还有许多可以探索的主题。

接下来，我们将深入了解统计学的世界。我们将首先使用更正式和系统化的方法来理解统计学的基本知识，然后理解概率和变量的作用，最后将这些概念联系起来以实施统计建模。

YKA34

LHL39