srcmini本文概述
当你想用Python进行科学工作时, 你可以使用的第一个库是SciPy。正如你将在本教程中看到的那样, SciPy不仅是一个库, 而且是一个完整的库生态系统, 这些生态系统协同工作以帮助你快速而可靠地完成复杂的科学任务。
在本教程中, 你将学习如何：
- 查找有关使用SciPy可以做的所有事情的信息
- 在计算机上安装SciPy
- 使用SciPy通过多个变量对数据集进行聚类
- 使用SciPy寻找最佳功能
让我们深入了解SciPy的美好世界！
免费奖金：单击此处可访问Conda速查表, 并提供用于管理Python环境和软件包的便捷用法示例。
移除广告
区分科学生态系统和科学库
当你想将Python用于科学计算任务时, 建议你使用多个库, 其中包括：
- NumPy
- 科学
- Matplotlib
- IPython的
- symbol
- pandas
这些库共同构成了SciPy生态系统, 并旨在协同工作。其中许多直接依赖于NumPy数组进行计算。本教程希望你熟悉创建NumPy数组并对其进行操作。
注意：如果需要NumPy上的快速入门或进阶知识, 则可以查看以下教程：
- 看起来不错, 没有麻烦：使用NumPy进行数组编程
- NumPy arange()：如何使用np.arange()
- MATLAB vs.Python：基本数组操作概述
在本教程中, 你将学习SciPy库, 这是SciPy生态系统的核心组件之一。 SciPy库是Python中科学计算的基础库。它为诸如数字积分, 优化, 信号处理, 线性代数等任务提供了许多高效且用户友好的界面。
了解SciPy模块
SciPy库由许多模块组成, 这些模块将库分成不同的功能单元。如果你想了解SciPy包括的不同模块, 则可以在scipy上运行help(), 如下所示：
>>>
>>> import scipy
>>> help(scipy)
这会为整个SciPy库产生一些帮助输出, 下面显示了其中的一部分：
Subpackages
-----------
Using any of these subpackages requires an explicit import. For example, ``import scipy.cluster``.
::
cluster --- Vector Quantization / Kmeans
fft --- Discrete Fourier transforms
fftpack --- Legacy discrete Fourier transforms
integrate --- Integration routines
...
此代码块显示了帮助输出的“子包”部分, 该列表是SciPy中可用于计算的所有可用模块的列表。
请注意本节顶部的文本, 其中指出：“使用这些子软件包中的任何一个都需要显式导入。”如果要使用SciPy中模块的功能, 则需要导入要专门使用的模块。在本教程的稍后部分, 你将看到一些示例, SciPy文档中显示了从SciPy导入库的准则。
一旦确定了要使用的模块, 就可以签出SciPy API参考, 其中包含SciPy中每个模块的所有详细信息。如果你需要更多说明, SciPy讲义是深入了解许多SciPy模块的好资源。
在本教程的后面, 你将学习集群和优化, 这是SciPy库中的两个模块。但是首先, 你需要在计算机上安装SciPy。
在计算机上安装SciPy
与大多数Python软件包一样, 有两种主要方法可以在计算机上安装SciPy：
- 水蟒
- PyPI和pip
在这里, 你将学习如何使用这两种方法来安装库。 SciPy唯一的直接依赖项是NumPy软件包。如有必要, 这两种安装方法都会自动安装除SciPy之外的NumPy。
水蟒
Anaconda是Python的流行发行版, 主要是因为它包含适用于Windows, macOS和Linux的最受欢迎的科学Python软件包的预构建版本。如果你尚未在计算机上安装Python, 那么Anaconda是入门的绝佳选择。 Anaconda预先安装了SciPy及其必需的依赖项, 因此, 一旦安装了Anaconda, 就无需执行其他任何操作！
你可以从他们的下载页面下载并安装Anaconda。确保下载最新的Python 3版本。将安装程序安装到计算机上之后, 即可根据应用程序的平台执行应用程序的默认设置过程。
注意：确保将Anaconda安装在不需要管理员权限即可修改的目录中。这是安装程序中的默认设置。
如果你已经安装了Anaconda, 但是你想安装或更新SciPy, 那么你也可以这样做。在macOS或Linux上打开终端应用程序, 或在Windows上打开Anaconda Prompt, 然后键入以下代码行之一：
$ conda install scipy
$ conda update scipy
如果需要安装SciPy, 则应使用第一行；如果只想更新SciPy, 则应使用第二行。为了确保已安装SciPy, 请在终端中运行Python并尝试导入SciPy：
>>>
>>> import scipy
>>> print(scipy.__file__)
/.../lib/python3.7/site-packages/scipy/__init__.py
在此代码中, 你已导入scipy并打印了从中加载scipy的文件的位置。上面的示例适用于macOS。你的计算机可能会显示其他位置。现在, 你已经在计算机上安装了SciPy, 可以使用了。你可以跳到下一部分, 开始使用SciPy！
移除广告
点子
如果你已经安装了不是Anaconda的Python版本, 或者你不想使用Anaconda, 那么你将使用pip来安装SciPy。要了解有关什么是点的更多信息, 请查看什么是点？新Pythonista指南。
注意：pip使用称为wheel的格式安装软件包。采用轮式格式, 代码在发送到你的计算机之前会先进行编译。尽管轮格式文件与Anaconda格式略有不同, 并且两者不可互换, 但这几乎与Anaconda所采用的方法相同。
要使用pip安装SciPy, 请打开终端应用程序, 然后键入以下代码行：
$ python -m pip install -U scipy
如果尚未安装该代码, 则将安装SciPy；如果已安装, 则将升级SciPy。为了确保已安装SciPy, 请在终端中运行Python并尝试导入SciPy：
>>>
>>> import scipy
>>> print(scipy.__file__)
/.../lib/python3.7/site-packages/scipy/__init__.py
在此代码中, 你已导入scipy并打印了从中加载scipy的文件的位置。上面的示例适用于使用pyenv的macOS。你的计算机可能会显示其他位置。现在, 你已经在计算机上安装了SciPy。让我们看看如何使用SciPy解决你可能遇到的几个问题！
在SciPy中使用群集模块
聚类是一种流行的技术, 它通过将数据关联到组来对数据进行分类。 SciPy库包含k-means聚类算法以及几种分层聚类算法的实现。在此示例中, 你将在scipy.cluster.vq中使用k-means算法, 其中vq表示矢量量化。
首先, 你应该看一下将用于此示例的数据集。数据集包含4827个真实邮件和747个垃圾短信(或SMS)消息。原始数据集可以在UCI机器学习存储库或作者的网页上找到。
注意：这些数据是由Tiago A. Almeida和JoséMaríaGómezHidalgo收集的, 并发表在2011年ACM文档工程研讨会(DOCENG)上的标题为“ SMS垃圾邮件过滤研究的贡献：新收集和结果”的文章中’11)于2011年在美国加利福尼亚山景城举办。
在数据集中, 每个消息都具有两个标签之一：
- 火腿提供合法信息
- 垃圾邮件
全文消息与每个标签关联。扫描数据时, 你可能会发现垃圾邮件中往往包含很多数字。它们通常包括电话号码或奖金。让我们根据邮件中的位数来预测邮件是否为垃圾邮件。为此, 你将根据消息中显示的位数将数据分为三组：
- 不是垃圾邮件：预计位数最少的邮件不是垃圾邮件。
- 未知：具有中间位数的消息是未知的, 需要由更高级的算法处理。
- 垃圾邮件：位数最多的邮件预计为垃圾邮件。
让我们开始对短信进行聚类。首先, 你应该导入在此示例中使用的库：
1 from pathlib import Path
2 import numpy as np
3 from scipy.cluster.vq import whiten, kmeans, vq
你会看到你正在从scipy.cluster.vq导入三个函数。这些函数中的每一个都接受NumPy数组作为输入。这些数组应在列中具有数据集的特征, 在行中具有观测值。
要素是你感兴趣的变量, 而每次你记录每个要素时都会创建一个观察值。在此示例中, 数据集中有5, 574个观测值或单个消息。此外, 你会发现有两个功能：
- 短信中的位数
- 该位数出现在整个数据集中的次数
接下来, 你应该从UCI数据库加载数据文件。数据以文本文件形式出现, 其中消息的类别通过制表符与消息分开, 并且每条消息都位于单独的行上。你应该使用pathlib.Path将数据读入列表：
4 data = Path("SMSSpamCollection").read_text()
5 data = data.strip()
6 data = data.split("\n")
在此代码中, 你使用pathlib.Path.read_text()将文件读取为字符串。然后, 使用.strip()删除所有尾随空格, 并使用.split()将字符串拆分为一个列表。
接下来, 你可以开始分析数据。你需要计算每条短信中出现的位数。 Python在标准库中包含collections.Counter, 用于收集类字典结构的对象计数。但是, 由于scipy.cluster.vq中的所有函数都希望NumPy数组作为输入, 因此在此示例中, 你不能使用collections.Counter。相反, 你可以使用NumPy数组并手动实现计数。
同样, 你对给定SMS消息中的位数以及有多少位数的SMS消息感兴趣。首先, 你应该创建一个NumPy数组, 该数组将给定消息中的位数与该消息的结果相关联, 无论是火腿还是垃圾邮件：
7 digit_counts = np.empty((len(data), 2), dtype=int)
在此代码中, 你将创建一个空的NumPy数组digit_counts, 其中包含两列和5, 574行。行数等于数据集中的消息数。你将使用digit_counts将邮件中的位数与邮件是否为垃圾邮件相关联。
你应该在进入循环之前创建阵列, 这样就不必在阵列扩展时分配新的内存。这样可以提高代码的效率。接下来, 你应该处理数据以记录数字位数和消息状态：
8 for i, line in enumerate(data):
9 case, message = line.split("\t")
10 num_digits = sum(c.isdigit() for c in message)
11 digit_counts[i, 0] = 0 if case == "ham" else 1
12 digit_counts[i, 1] = num_digits
以下是此代码的工作原理的逐行细分：
- 第8行：循环遍历数据。你可以使用enumerate()将列表中的值放入行中, 并为此列表创建索引i。要了解有关enumerate()的更多信息, 请签出使用enumerate()保持运行索引。
- 第9行：在制表符上拆分行以创建大小写和消息。 case是说明邮件是垃圾邮件还是垃圾邮件的字符串, 而message是包含邮件文本的字符串。
- 第10行：使用理解值的sum()计算消息中的位数。在理解中, 你使用isdigit()检查消息中的每个字符, 如果元素是数字, 则返回True, 否则返回False。 sum()然后将每个True结果视为1, 将每个False视为0。因此, 此理解中sum()的结果是isdigit()返回True的字符数。
- 第11行：将值分配给digit_counts。如果邮件是合法的(火腿), 则将i行的第一列分配为0；如果邮件是垃圾邮件, 则将其分配为1。
- 第12行：将值分配给digit_counts。你将i行的第二列指定为消息中的位数。
现在, 你有一个NumPy数组, 其中包含每个消息中的位数。但是, 你想将聚类算法应用于具有一定位数的消息数量的数组。换句话说, 你需要创建一个数组, 其中第一列具有消息中的位数, 而第二列是具有该位数的消息数。查看以下代码：
13 unique_counts = np.unique(digit_counts[:, 1], return_counts=True)
np.unique()将数组作为第一个参数, 并返回另一个具有该参数唯一元素的数组。它还需要几个可选参数。在这里, 你使用return_counts = True来指示np.unique()还返回一个数组, 其中包含输入数组中每个唯一元素出现的次数。这两个输出将作为元组返回, 并存储在unique_counts中。
接下来, 你需要将unique_counts转换为适合聚类的形状：
14 unique_counts = np.transpose(np.vstack(unique_counts))
你可以使用np.vstack()将来自np.unique()的两个1xN输出合并为一个2xN数组, 然后将它们转置为Nx2数组。此格式是你在聚类功能中使用的格式。现在, unique_counts中的每一行都有两个元素：
- 消息中的位数
- 具有该数字位数的邮件数
这两个操作的输出的子集如下所示：
[[ 0 4110]
[ 1 486]
[ 2 160]
...
[ 40 4]
[ 41 2]
[ 47 1]]
在数据集中, 有4110条没有数字的消息, 486条没有数字的消息, 依此类推。现在, 你应该将k-means聚类算法应用于此数组：
15 whitened_counts = whiten(unique_counts)
16 codebook, _ = kmeans(whitened_counts, 3)
你可以使用whiten()将每个特征归一化以具有单位方差, 从而改善kmeans()的结果。然后, kmeans()将变白的数据和要创建的簇数作为参数。在此示例中, 你要创建3个群集, 分别是绝对火腿, 绝对垃圾邮件和未知簇。 kmeans()返回两个值：
- 具有三行两列的数组, 代表每个组的质心：kmeans()算法通过最小化观测值到每个质心的距离来计算每个聚类的质心的最佳位置。该数组分配给代码本。
- 从观测值到质心的平均欧几里得距离：在本示例的其余部分中, 你将不需要该值, 因此可以将其分配给_。
接下来, 你应该使用vq()确定每个观察值属于哪个群集：
17 codes, _ = vq(unique_counts, codebook)
vq()将代码簿中的代码分配给每个观察值。它返回两个值：
- 第一个值是与unique_counts长度相同的数组, 其中每个元素的值是一个整数, 表示观察值分配给的簇。由于在此示例中使用了三个聚类, 因此将每个观察值分配给聚类0、1或2。
- 第二个值是每个观测值与其质心之间的欧几里得距离的数组。
现在你已经将数据聚类了, 你应该使用它来对SMS消息进行预测。你可以检查计数, 以确定聚类算法将多少位数划定为绝对垃圾邮件和未知垃圾之间以及未知垃圾邮件和垃圾邮件之间的界限。
聚类算法会为每个聚类随机分配代码0、1或2, 因此你需要确定哪个是哪个。你可以使用以下代码查找与每个集群关联的代码：
18 ham_code = codes[0]
19 spam_code = codes[-1]
20 unknown_code = list(set(range(3)) ^ set((ham_code, spam_code)))[0]
在此代码中, 第一行找到与火腿消息关联的代码。根据上面的假设, 火腿消息的位数最少, 并且数字数组的编号从最小到最大。因此, 火腿消息簇从代码的开头开始。
同样, 垃圾邮件消息的位数最多, 并且形成代码中的最后一个簇。因此, 垃圾邮件的代码将等于代码的最后一个元素。最后, 你需要找到未知消息的代码。由于该代码只有3个选项, 并且你已经标识了其中两个, 因此可以在Python集上使用symmetric_difference运算符来确定最后一个代码值。然后, 你可以打印与每种消息类型关联的群集：
21 print("definitely ham:", unique_counts[codes == ham_code][-1])
22 print("definitely spam:", unique_counts[codes == spam_code][-1])
23 print("unknown:", unique_counts[codes == unknown_code][-1])
在此代码中, 每一行都获得了unique_counts中的行, 其中vq()分配了不同的代码值。由于该操作返回一个数组, 因此你应该获取数组的最后一行以确定分配给每个群集的最高位数。输出如下所示：
definitely ham: [0 4110]
definitely spam: [47 1]
unknown: [20 18]
在此输出中, 你看到肯定是火腿消息是消息中数字为零的消息, 未知消息是1到20位之间的所有内容, 垃圾邮件肯定是21到47位之间的所有内容, 这是最大数量数据集中的数字。
现在, 你应该检查预测在此数据集上的准确性。首先, 为digit_counts创建一些掩码, 以便你可以轻松获取邮件的垃圾邮件或垃圾邮件状态：
24 digits = digit_counts[:, 1]
25 predicted_hams = digits == 0
26 predicted_spams = digits > 20
27 predicted_unknowns = np.logical_and(digits > 0, digits <= 20)
在此代码中, 你正在创建预测的_hams掩码, 其中消息中没有数字。然后, 为所有超过20位数字的邮件创建预测的垃圾邮件掩码。最后, 中间的消息是predicted_unknowns。
接下来, 将这些掩码应用于实际数字计数以检索预测：
28 spam_cluster = digit_counts[predicted_spams]
29 ham_cluster = digit_counts[predicted_hams]
30 unk_cluster = digit_counts[predicted_unknowns]
在这里, 你要将在最后一个代码块中创建的掩码应用于digit_counts数组。这将创建三个新阵列, 其中仅包含已集群到每个组中的消息。最后, 你可以看到每种群集中有多少种消息类型：
31 print("hams:", np.unique(ham_cluster[:, 0], return_counts=True))
32 print("spams:", np.unique(spam_cluster[:, 0], return_counts=True))
33 print("unknowns:", np.unique(unk_cluster[:, 0], return_counts=True))
此代码显示集群中每个唯一值的计数。请记住, 0表示邮件为垃圾邮件, 1表示邮件为垃圾邮件。结果如下所示：
hams: (array([0, 1]), array([4071, 39]))
spams: (array([0, 1]), array([ 1, 232]))
unknowns: (array([0, 1]), array([755, 476]))
从此输出中, 你可以看到4110条消息属于绝对火腿组, 其中4071条实际上是火腿, 只有39条是垃圾邮件。相反, 在属于绝对垃圾邮件组的233条消息中, 只有1条实际上是垃圾邮件, 其余均为垃圾邮件。
当然, 超过1200条消息属于未知类别, 因此需要进行更高级的分析以对这些消息进行分类。你可能希望研究自然语言处理等方法来帮助提高预测的准确性, 并且可以使用Python和Keras来提供帮助。
移除广告
在SciPy中使用优化模块
当你需要优化函数的输入参数时, scipy.optimize包含许多用于优化各种函数的有用方法：
- minimal_scalar()和minimum()分别最小化一个变量和许多变量的函数
- curve_fit() to fit a function to a set of data
- root_scalar()和root()分别查找一个变量和多个变量的函数的零
- linprog() to minimize a linear objective function with linear inequality and equality constraints
实际上, 所有这些功能都在执行一种或另一种优化。在本节中, 你将学习两个最小化函数, minimum_scalar()和minimum()。
用一个变量最小化一个函数
接受一个数字并产生一个输出的数学函数称为标量函数。通常将它与接受多个数字并产生多个输出的多元函数形成对比。在下一部分中, 你将看到优化多元函数的示例。
对于此部分, 标量函数将是四次多项式, 并且你的目标是找到函数的最小值。该函数为y =3x⁴-2x +1。在下图中从0到1的x范围绘制函数。
在图中, 你可以看到此函数的最小值约为x = 0.55。你可以使用minimum_scalar()确定最小值的确切x和y坐标。首先, 从scipy.optimize导入minimal_scalar()。然后, 你需要定义要最小化的目标函数：
1 from scipy.optimize import minimize_scalar
2
3 def objective_function(x):
4 return 3 * x ** 4 - 2 * x + 1
Objective_function()接受输入x并对其进行必要的数学运算, 然后返回结果。在函数定义中, 可以使用所需的任何数学函数。唯一的限制是该函数必须在最后返回一个数字。
接下来, 使用minimum_scalar()查找此函数的最小值。 minimal_scalar()仅具有一个必需的输入, 这是目标函数定义的名称：
5 res = minimize_scalar(objective_function)
minimal_scalar()的输出是OptimizeResult的一个实例。该类从优化器的运行中收集了许多相关的详细信息, 包括优化是否成功以及成功的最终结果。该函数的minimum_scalar()的输出如下所示：
fun: 0.17451818777634331
nfev: 16
nit: 12
success: True
x: 0.5503212087491959
这些结果都是OptimizeResult的所有属性。 success是一个布尔值, 指示优化是否成功完成。如果优化成功, 则有趣的是目标函数在最佳值x处的值。从输出中可以看到, 正如预期的那样, 此函数的最佳值接近x = 0.55。
注意：你可能知道, 并非每个功能都有最小值。例如, 尝试看看如果你的目标函数为y =x³会发生什么。对于minimal_scalar(), 没有最小值的目标函数通常会导致OverflowError, 因为优化器最终会尝试使用太大的数字以致无法由计算机计算。
在没有最小值的函数的另一侧是具有多个最小值的函数。在这些情况下, 不能保证minimal_scalar()会找到函数的全局最小值。但是, minimum_scalar()具有方法关键字参数, 你可以指定该参数来控制用于优化的求解器。 SciPy库具有三种内置的标量最小化方法：
布伦特
是一个实现
布伦特算法
。此方法是默认方法。
金色的
是
黄金分割搜索
。的
文件资料
注意到布伦特的方法通常更好。
有界的
是布伦特算法的有界实现。当最小值在已知范围内时, 限制搜索区域很有用。
当方法为brent或golden时, minimum_scalar()将使用另一个称为方括号的参数。这是一个由两个或三个元素组成的序列, 可为具有最小值的区域边界提供初始猜测。但是, 这些求解器不能保证找到的最小值在此范围内。
另一方面, 当有界方法时, minimum_scalar()将使用另一个称为界的参数。这是两个元素的序列, 这些元素严格限制了搜索区域的最小值。尝试使用函数y =x⁴-x²的有界方法。该功能如下图所示：
使用前面的示例代码, 你可以像这样重新定义objective_function()：
7 def objective_function(x):
8 return x ** 4 - x ** 2
首先, 尝试默认的brent方法：
9 res = minimize_scalar(objective_function)
在这段代码中, 你没有传递方法的值, 因此默认情况下, minimum_scalar()使用了brent方法。输出是这样的：
fun: -0.24999999999999994
nfev: 15
nit: 11
success: True
x: 0.7071067853059209
你可以看到优化成功。发现x = 0.707和y = -1/4附近的最优值。如果你通过解析求出方程的最小值, 则将在x = 1 /√2处找到最小值, 该最小值与最小化函数找到的答案非常接近。但是, 如果要在x = -1 /√2处找到对称最小值, 该怎么办？你可以通过向brent方法提供方括号参数来返回相同的结果：
10 res = minimize_scalar(objective_function, bracket=(-1, 0))
在此代码中, 你提供了序列(-1, 0)括起来以在-1和0之间的区域中开始搜索。由于目标函数关于y轴对称, 因此你期望在该区域中有一个最小值。 。但是, 即使使用括号, brent方法仍会在x = + 1 /√2时返回最小值。要在x = -1 /√2处找到最小值, 可以将bounded方法与bounds一起使用：
11 res = minimize_scalar(objective_function, method='bounded', bounds=(-1, 0))
在此代码中, 将方法和界限添加为minimal_scalar()的参数, 并将界限设置为-1到0。此方法的输出如下：
fun: -0.24999999999998732
message: 'Solution found.'
nfev: 10
status: 0
success: True
x: -0.707106701474177
如预期的那样, 最小值在x = -1 /√2处发现。请注意此方法的附加输出, 其中包括res中的message属性。此字段通常用于某些最小化求解器的更详细的输出。
移除广告
最小化具有多个变量的函数
scipy.optimize还包括更通用的minimum()。此函数可以处理多变量输入和输出, 并且具有更复杂的优化算法才能处理此问题。此外, maximum()可以处理对你的问题的解决方案的约束。你可以指定三种约束类型：
- LinearConstraint：通过将解决方案x值的内积与用户输入数组进行比较, 并将结果与上下边界进行比较, 来约束解决方案。
- NonlinearConstraint：通过对解决方案x值应用用户提供的函数并将返回值与上下限进行比较, 来约束解决方案。
- 界限：解决方案x值被约束在上下限之间。
使用这些约束时, 它可能会限制你可以使用的优化方法的特定选择, 因为并非所有可用的方法都以这种方式支持约束。
让我们尝试演示如何使用minimal()。假设你是一位股票经纪人, 并且希望通过出售固定数量的股票来获得总收入。你已经确定了一组特定的买家, 并且对于每个买家, 你都知道他们要支付的价格以及他们手中有多少现金。
你可以将此问题表述为约束优化问题。目标功能是你要最大化你的收入。但是, maximum()会找到函数的最小值, 因此你需要将目标函数乘以-1才能找到产生最大负数的x值。
这个问题有一个限制, 那就是买方购买的总股份总数不超过你手头的股份数目。每个解决方案变量也都有界限, 因为每个购买者都有可用现金上限, 而下限为零。解决方案的x值负值表示你要向买家付款！
试试下面的代码来解决这个问题。首先, 导入所需的模块, 然后设置变量以确定市场上的购买者数量和要出售的股票数量：
1 import numpy as np
2 from scipy.optimize import minimize, LinearConstraint
3
4 n_buyers = 10
5 n_shares = 15
在此代码中, 你从scipy.optimize导入numpy, minimum()和LinearConstraint。然后, 你设定一个由10个买家组成的市场, 这些买家将总共向你购买15股。
接下来, 给定前两个数组, 创建数组以存储每个买方支付的价格, 他们可以负担的最大金额以及每个买方可以负担的最大股票数量。对于此示例, 可以在np.random中使用随机数生成来生成数组：
6 np.random.seed(10)
7 prices = np.random.random(n_buyers)
8 money_available = np.random.randint(1, 4, n_buyers)
在此代码中, 你为NumPy的随机数生成器设置了种子。此功能可确保你每次运行此代码时, 都会获得相同的随机数集。这里是为了确保你的输出与用于比较的教程相同。
在第7行中, 你生成了买家将支付的价格数组。 np.random.random()在半开间隔[0, 1)上创建一个随机数数组。数组中元素的数量由参数的值确定, 在这种情况下, 参数的数量是购买者的数量。
在第8行中, 你从[1, 4)的半开间隔生成了一个整数数组, 同样具有购买者数量的大小。此数组代表每个买家可用的总现金。现在, 你需要计算每个买家可以购买的最大股票数量：
9 n_shares_per_buyer = money_available / prices
10 print(prices, money_available, n_shares_per_buyer, sep="\n")
在第9行中, 你使用money_available与价格的比率来确定每个买家可以购买的最大股票数量。最后, 打印每个用换行符分隔的数组。输出如下所示：
[0.77132064 0.02075195 0.63364823 0.74880388 0.49850701 0.22479665
0.19806286 0.76053071 0.16911084 0.08833981]
[1 1 1 3 1 3 3 2 1 1]
[ 1.29647768 48.18824404 1.57816269 4.00638948 2.00598984 13.34539487
15.14670609 2.62974258 5.91328161 11.3199242 ]
第一行是价格数组, 它们是0到1之间的浮点数。此行后是1到4之间的整数形式的最大可用现金量。最后, 你看到每个买家可以购买的股票数量。
现在, 你需要为求解器创建约束和界限。约束条件是购买的股份总数之和不能超过可用股份总数。这是一个约束而非约束, 因为它涉及多个解决方案变量。
为了用数学表示这一点, 你可以说x [0] + x [1] + … + x [n] = n_shares, 其中n是购买者的总数。更简洁地说, 你可以将矢量的点或内积与解值相乘, 并将其约束为等于n_shares。请记住, LinearConstraint将输入数组的点积与解值相乘, 并将其与上下限进行比较。你可以使用它在n_shares上设置约束：
11 constraint = LinearConstraint(np.ones(n_buyers), lb=n_shares, ub=n_shares)
在此代码中, 你将创建一个长度为n_buyers的数组, 并将其作为第一个参数传递给LinearConstraint。由于LinearConstraint将此参数与解向量的点积相乘, 因此将得出已购买股票的总和。
然后将此结果约束为位于其他两个参数之间：
- 下界磅
- 上限ub
由于lb = ub = n_shares, 因此这是一个相等约束, 因为值的总和必须等于lb和ub。如果lb与ub不同, 那么它将是一个不平等约束。
接下来, 为解决方案变量创建界限。边界将购买的股票数量限制为下限为0, 上限为n_shares_per_buyer。 minimal()期望的边界格式是上下限元组的序列：
12 bounds = [(0, n) for n in n_shares_per_buyer]
在此代码中, 你将使用一个理解为每个购买者生成一个元组列表。运行优化之前的最后一步是定义目标函数。回想一下, 你正在尝试使收入最大化。同样, 你希望使收入的负数尽可能大地成为负数。
你从每次销售中获得的收入是买方支付的价格乘以他们所购买的股票数量。数学上, 你可以将其写为prices [0] * x [0] + prices [1] * x [1] + … + prices [n] * x [n], 其中n再次是购买者的总数。
再一次, 你可以使用内积或x.dot(prices)来更简洁地表示这一点。这意味着你的目标函数应将当前解决方案值x和价格数组作为参数：
13 def objective_function(x, prices):
14 return -x.dot(prices)
在此代码中, 你定义了objective_function()以接受两个参数。然后, 将x的点乘积与价格相乘, 然后返回该值的负数。请记住, 你必须返回负数, 因为你正试图使该数字尽可能小或尽可能接近负无穷大。最后, 你可以调用minimum()：
15 res = minimize(
16 objective_function, 17 x0=10 * np.random.random(n_buyers), 18 args=(prices, ), 19 constraints=constraint, 20 bounds=bounds, 21 )
在此代码中, res是OptimizeResult的一个实例, 就像minimal_scalar()一样。就像你将看到的, 即使问题完全不同, 也存在许多相同的字段。在对minimum()的调用中, 你传递了五个参数：
- Objective_function：第一个位置参数必须是你要优化的函数。
- x0：下一个参数是对解决方案值的初始猜测。在这种情况下, 你只是提供了一个介于0到10之间且长度为n_buyers的随机值数组。对于某些算法或某些问题, 选择适当的初始猜测可能很重要。但是, 对于此示例, 它似乎并不重要。
- args：下一个参数是必须传递给目标函数的其他参数的元组。 minimal()总是将解x的当前值传递到目标函数中, 因此此参数用作收集任何其他必要输入的地方。在此示例中, 你需要将价格传递给objective_function(), 这样就可以了。
- 约束：下一个参数是对问题的一系列约束。你传递的是你先前生成的关于可用份额数的约束。
- bounds：最后一个参数是你先前生成的解决方案变量的界限序列。
一旦求解器运行, 你应该通过打印来检查res：
fun: -8.783020157087478
jac: array([-0.7713207 , -0.02075195, -0.63364828, -0.74880385, -0.49850702, -0.22479665, -0.1980629 , -0.76053071, -0.16911089, -0.08833981])
message: 'Optimization terminated successfully.'
nfev: 204
nit: 17
njev: 17
status: 0
success: True
x: array([1.29647768e+00, 2.78286565e-13, 1.57816269e+00, 4.00638948e+00, 2.00598984e+00, 3.48323773e+00, 3.19744231e-14, 2.62974258e+00, 2.38121197e-14, 8.84962214e-14])
在此输出中, 你可以看到消息和状态, 指示优化的最终状态。对于此优化器, 状态为0表示优化已成功终止, 你也可以在消息中看到该状态。由于优化成功, 因此fun以优化的解决方案值显示了目标函数的值。你将从此次销售中获得$ 8.78的收入。
你可以在res.x中看到优化函数的x值。在这种情况下, 结果是你应该向第一个买家出售约1.3股, 向第二个买家出售约0股, 向第三个买家出售1.6股, 向第四个买家出售4.0股, 依此类推。
你还应该检查并确保满足你设置的约束和界限。你可以使用以下代码执行此操作：
22 print("The total number of shares is:", sum(res.x))
23 print("Leftover money for each buyer:" money_available - res.x * prices)
在此代码中, 你将打印每个买家购买的股票总数, 该总数应等于n_shares。然后, 你打印每个买家的手头现金与他们花费的金额之间的差额。这些值均应为正。这些检查的输出如下所示：
The total number of shares is: 15.000000000000002
Leftover money for each buyer: [3.08642001e-14 1.00000000e+00
3.09752224e-14 6.48370246e-14 3.28626015e-14 2.21697984e+00
3.00000000e+00 6.46149800e-14 1.00000000e+00 1.00000000e+00]
如你所见, 解决方案上的所有约束和界限都得到了满足。现在, 你应该尝试解决问题, 以便解决者无法找到解决方案。将n_shares的值更改为1000, 以便你尝试向这些相同的买方出售1000股。当运行minimal()时, 你将发现结果如下所示：
fun: nan
jac: array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan])
message: 'Iteration limit exceeded'
nfev: 2182
nit: 101
njev: 100
status: 9
success: False
x: array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan])
请注意, 现在status属性的值为9, 并且该消息指出已超过迭代限制。考虑到每个买家的钱数和市场上的买家数量, 无法出售1000股股票。但是, minimum()并不会引发错误, 而是返回OptimizeResult实例。在进行进一步的计算之前, 你需要确保检查状态码。
移除广告
结论
在本教程中, 你了解了SciPy生态系统及其与SciPy库的区别。你了解了SciPy中可用的一些模块, 并了解了如何使用Anaconda或pip安装SciPy。然后, 你重点研究一些使用SciPy中的群集和优化功能的示例。
在集群示例中, 你开发了一种算法, 可以根据合法邮件对垃圾邮件进行分类。使用kmeans(), 你发现超过20位数字的邮件极有可能是垃圾邮件！
在优化示例中, 你首先在只有一个变量的数学上清晰的函数中找到了最小值。然后, 你解决了更复杂的问题, 即通过卖出股票获得最大利润。使用minimal(), 你发现了要出售给一组购买者的最佳股票数量, 并获得了8.79美元的利润！
SciPy是一个巨大的库, 还有许多其他模块可供学习。有了你现在所掌握的知识, 你就可以开始探索！
评论前必须登录！
注册