原文：
zh.annas-archive.org/md5/123a7612a4e578f6816d36f968cfec22译者：飞龙
第八章：回归和预测
统计学家或数据科学家最重要的任务之一是生成对两组数据之间关系的系统性理解。这可能意味着两组数据之间的“连续”关系，其中一个值直接取决于另一个变量的值。或者，它可能意味着分类关系，其中一个值根据另一个值进行分类。处理这些问题的工具是回归。在其最基本的形式中，回归涉及将一条直线拟合到两组数据的散点图中，并进行一些分析，以查看这条直线如何“拟合”数据。当然，我们通常需要更复杂的东西来模拟现实世界中存在的更复杂的关系。
时间序列代表这些回归类型问题的一种专门类别，其中一个值在一段时间内发展。与更简单的问题不同，时间序列数据通常在连续值之间有复杂的依赖关系；例如，一个值可能依赖于前两个值，甚至可能依赖于前一个“噪音”。时间序列建模在科学和经济学中非常重要，有各种工具可用于建模时间序列数据。处理时间序列数据的基本技术称为ARIMA，它代表自回归综合移动平均。该模型包括两个基本组件，一个自回归（AR）组件和一个移动平均**（MA）组件，用于构建观察数据的模型。
在本章中，我们将学习如何建立两组数据之间的关系模型，量化这种关系的强度，并对其他值（未来）生成预测。然后，我们将学习如何使用逻辑回归，在分类问题中，这是简单线性模型的一种变体。最后，我们将使用 ARIMA 为时间序列数据构建模型，并基于这些模型构建不同类型的数据。我们将通过使用一个名为 Prophet 的库来自动生成时间序列数据模型来结束本章。
在本章中，我们将涵盖以下内容：
使用基本线性回归
使用多元线性回归
使用对数回归进行分类
使用 ARMA 对时间序列数据进行建模
使用 ARIMA 从时间序列数据进行预测
使用 ARIMA 对季节性数据进行预测
使用 Prophet 对时间序列进行建模
让我们开始吧！
技术要求
在本章中，像往常一样，我们需要导入 NumPy 包并使用别名np，导入 Matplotlib pyplot模块并使用plt别名，以及导入 Pandas 包并使用pd别名。我们可以使用以下命令来实现：
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
在本章中，我们还需要一些新的包。statsmodels 包用于回归和时间序列分析，scikit-learn包（sklearn）提供通用数据科学和机器学习工具，Prophet 包（fbprophet）用于自动生成时间序列数据模型。这些包可以使用您喜欢的包管理器（如pip）进行安装：
python3.8 -m pip install statsmodels sklearn fbprophet
Prophet 包可能在某些操作系统上安装起来比较困难，因为它的依赖关系。如果安装fbprophet出现问题，您可能希望尝试使用 Python 的 Anaconda 发行版及其包管理器conda，它可以更严格地处理依赖关系：
conda install fbprophet
最后，我们还需要一个名为tsdata的小模块，该模块包含在本章的存储库中。该模块包含一系列用于生成样本时间序列数据的实用程序。
本章的代码可以在 GitHub 存储库的Chapter 07文件夹中找到：github.com/PacktPublishing/Applying-Math-with-Python/tree/master/Chapter%2007。
查看以下视频以查看代码实际操作：bit.ly/2Ct8m0B。
使用基本线性回归
线性回归是一种建模两组数据之间依赖关系的工具，这样我们最终可以使用这个模型进行预测。名称来源于我们基于第二组数据形成一个线性模型（直线）。在文献中，我们希望建模的变量通常被称为响应变量，而我们在这个模型中使用的变量是预测变量。
在这个步骤中，我们将学习如何使用 statsmodels 包执行简单的线性回归，以建模两组数据之间的关系。
准备工作
对于这个步骤，我们需要导入 statsmodels 的api模块并使用别名sm，导入 NumPy 包并使用别名np，导入 Matplotlib 的pyplot模块并使用别名plt，以及一个 NumPy 默认随机数生成器的实例。所有这些都可以通过以下命令实现：
import statsmodels.api as sm
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import default_rng
rng = default_rng(12345)
如何做到这一点…
以下步骤概述了如何使用 statsmodels 包对两组数据执行简单线性回归：
- 首先，我们生成一些示例数据进行分析。我们将生成两组数据，这将说明一个良好的拟合和一个不太好的拟合：
x = np.linspace(0, 5, 25)
rng.shuffle(x)
trend = 2.0
shift = 5.0
y1 = trend*x + shift + rng.normal(0, 0.5, size=25)
y2 = trend*x + shift + rng.normal(0, 5, size=25)
- 执行回归分析的一个很好的第一步是创建数据集的散点图。我们将在同一组坐标轴上完成这一步：
fig, ax = plt.subplots()
ax.scatter(x, y1, c="b", label="Good correlation")
ax.scatter(x, y2, c="r", label="Bad correlation")
ax.legend()
ax.set_xlabel("X"),
ax.set_ylabel("Y")
ax.set_title("Scatter plot of data with best fit lines")
- 我们需要使用
sm.add_constant实用程序例程，以便建模步骤将包括一个常数值：
pred_x = sm.add_constant(x)
- 现在，我们可以为我们的第一组数据创建一个
OLS模型，并使用fit方法来拟合模型。然后，我们使用summary方法打印数据的摘要：
model1 = sm.OLS(y1, pred_x).fit()
print(model1.summary())
- 我们重复第二组数据的模型拟合并打印摘要：
model2 = sm.OLS(y2, pred_x).fit()
print(model2.summary())
- 现在，我们使用
linspace创建一个新的x值范围，我们可以用它来在散点图上绘制趋势线。我们需要添加constant列以与我们创建的模型进行交互：
model_x = sm.add_constant(np.linspace(0, 5))
- 接下来，我们在模型对象上使用
predict方法，这样我们就可以使用模型在前一步生成的每个x值上预测响应值：
model_y1 = model1.predict(model_x)
model_y2 = model2.predict(model_x)
- 最后，我们将在散点图上绘制在前两个步骤中计算的模型数据：
ax.plot(model_x[:, 1], model_y1, 'b')
ax.plot(model_x[:, 1], model_y2, 'r')
散点图以及我们添加的最佳拟合线（模型）可以在下图中看到：
图 7.1：使用最小二乘法回归计算的数据散点图和最佳拟合线。
工作原理…
基本数学告诉我们，一条直线的方程如下所示：
在这里，c是直线与y轴相交的值，通常称为y截距，m是直线的斜率。在线性回归的背景下，我们试图找到响应变量Y和预测变量X之间的关系，使其形式为一条直线，使以下情况发生：
在这里，c和m现在是要找到的参数。我们可以用另一种方式来表示这一点，如下所示：
这里，E是一个误差项，一般来说，它取决于X。为了找到“最佳”模型，我们需要找到E误差被最小化的c和m参数值（在适当的意义上）。找到使这个误差最小化的参数值的基本方法是最小二乘法，这给了这里使用的回归类型以它的名字：普通最小二乘法。一旦我们使用这种方法建立了响应变量和预测变量之间的某种关系，我们的下一个任务就是评估这个模型实际上如何代表这种关系。为此，我们形成了以下方程给出的残差：
我们对每个数据点*X[i]和Y[i]*进行这样的操作。为了对我们对数据之间关系建模的准确性进行严格的统计分析，我们需要残差满足某些假设。首先，我们需要它们在概率意义上是独立的。其次，我们需要它们围绕 0 呈正态分布，并且具有相同的方差。（在实践中，我们可以稍微放松这些条件，仍然可以对模型的准确性做出合理的评论。）
在这个配方中，我们使用线性关系从预测数据中生成响应数据。我们创建的两个响应数据集之间的差异在于每个值的误差的“大小”。对于第一个数据集y1，残差呈正态分布，标准差为 0.5，而对于第二个数据集y2，残差的标准差为 5.0。我们可以在图 7.1中的散点图中看到这种变异性，其中y1的数据通常非常接近最佳拟合线 - 这与用于生成数据的实际关系非常接近 - 而y2数据则远离最佳拟合线。
来自 statsmodels 包的OLS对象是普通最小二乘回归的主要接口。我们将响应数据和预测数据作为数组提供。为了在模型中有一个常数项，我们需要在预测数据中添加一列 1。sm.add_constant例程是一个简单的实用程序，用于添加这个常数列。OLS类的fit方法计算模型的参数，并返回一个包含最佳拟合模型参数的结果对象（model1和model2）。summary方法创建一个包含有关模型和拟合优度的各种统计信息的字符串。predict方法将模型应用于新数据。顾名思义，它可以用于使用模型进行预测。
除了参数值本身之外，摘要中报告了另外两个统计量。第一个是R²值，或者调整后的版本，它衡量了模型解释的变异性与总变异性之间的关系。这个值将在 0 和 1 之间。较高的值表示拟合效果更好。第二个是 F 统计量的 p 值，它表示模型的整体显著性。与 ANOVA 测试一样，较小的 F 统计量表明模型是显著的，这意味着模型更有可能准确地对数据进行建模。
在这个配方中，第一个模型model1的调整后的R²值为 0.986，表明该模型非常紧密地拟合了数据，p 值为 6.43e-19，表明具有很高的显著性。第二个模型的调整后的R²值为 0.361，表明该模型与数据的拟合程度较低，p 值为 0.000893，也表明具有很高的显著性。尽管第二个模型与数据的拟合程度较低，但从统计学的角度来看，并不意味着它没有用处。该模型仍然具有显著性，尽管不如第一个模型显著，但它并没有解释所有的变异性（或者至少是数据中的一个显著部分）。这可能表明数据中存在额外的（非线性）结构，或者数据之间的相关性较低，这意味着响应和预测数据之间的关系较弱（由于我们构造数据的方式，我们知道后者是真实的）。
还有更多…
简单线性回归是统计学家工具包中一个很好的通用工具。它非常适合找到两组已知（或被怀疑）以某种方式相互关联的数据之间关系的性质。衡量一个数据集依赖于另一个数据集的程度的统计测量称为相关性。我们可以使用相关系数来衡量相关性，例如Spearman 秩相关系数。高正相关系数表示数据之间存在强烈的正相关关系，就像在这个示例中看到的那样，而高负相关系数表示强烈的负相关关系，其中通过数据的最佳拟合线的斜率为负。相关系数为 0 意味着数据没有相关性：数据之间没有关系。
如果数据集之间明显相关，但不是线性（直线）关系，那么它可能遵循一个多项式关系，例如，一个值与另一个值的平方有关。有时，您可以对一个数据集应用转换，例如对数，然后使用线性回归来拟合转换后的数据。当两组数据之间存在幂律关系时，对数特别有用。
使用多元线性回归
简单线性回归，如前面的示例所示，非常适合产生一个响应变量和一个预测变量之间关系的简单模型。不幸的是，有一个单一的响应变量依赖于许多预测变量更为常见。此外，我们可能不知道从一个集合中选择哪些变量作为良好的预测变量。对于这个任务，我们需要多元线性回归。
在这个示例中，我们将学习如何使用多元线性回归来探索响应变量和几个预测变量之间的关系。
准备就绪
对于这个示例，我们将需要导入 NumPy 包作为np，导入 Matplotlib 的pyplot模块作为plt，导入 Pandas 包作为pd，并创建 NumPy 默认随机数生成器的实例，使用以下命令：
from numpy.random import default_rng
rng = default_rng(12345)
我们还需要导入 statsmodels 的api模块作为sm，可以使用以下命令导入：
import statsmodels.api as sm
如何做…
以下步骤向您展示了如何使用多元线性回归来探索几个预测变量和一个响应变量之间的关系：
- 首先，我们需要创建要分析的预测数据。这将采用 Pandas 的
DataFrame形式，有四个项。我们将通过添加一个包含 1 的列来添加常数项：
p_vars = pd.DataFrame({
"const": np.ones((100,)),
"X1": rng.uniform(0, 15, size=100),
"X2": rng.uniform(0, 25, size=100),
"X3": rng.uniform(5, 25, size=100)
})
- 接下来，我们将仅使用前两个变量生成响应数据：
residuals = rng.normal(0.0, 12.0, size=100)
Y = -10.0 + 5.0*p_vars["X1"] - 2.0*p_vars["X2"] + residuals
- 现在，我们将生成响应数据与每个预测变量的散点图：
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True,
tight_layout=True)
ax1.scatter(p_vars["X1"], Y)
ax2.scatter(p_vars["X2"], Y)
ax3.scatter(p_vars["X3"], Y)
- 然后，我们将为每个散点图添加轴标签和标题，因为这是一个很好的做法：
ax1.set_title("Y against X1")
ax1.set_xlabel("X1")
ax1.set_ylabel("Y")
ax2.set_title("Y against X2")
ax2.set_xlabel("X2")
ax3.set_title("Y against X3")
ax3.set_xlabel("X3")
The resulting plots can be seen in the following figure:
Figure 7.2: Scatter plots of the response data against each of the predictor variables
As we can see, there appears to be some correlation between the response data and the first two predictor columns, `X1` and `X2`. This is what we expect, given how we generated the data.
5. We use the same `OLS` class to perform multilinear regression; that is, providing the response array and the predictor `DataFrame`:
model = sm.OLS(Y, p_vars).fit()
print(model.summary())
The output of the `print` statement is as follows:
OLS 回归结果
==================================================================
因变量：y R 平方：0.770
模型：OLS 调整 R 平方：0.762
方法：最小二乘法 F 统计量：106.8
日期：2020 年 4 月 23 日星期四 概率（F 统计量）：1.77e-30
时间：12:47:30 对数似然：-389.38
观测数量：100 AIC：786.8
残差自由度：96 BIC：797.2
模型自由度：3
协方差类型：非鲁棒
===================================================================
coef std err t P>|t| [0.025 0.975]
常数 -9.8676 4.028 -2.450 0.016 -17.863 -1.872
X1 4.7234 0.303 15.602 0.000 4.122 5.324
X2 -1.8945 0.166 -11.413 0.000 -2.224 -1.565
X3 -0.0910 0.206 -0.441 0.660 -0.500 0.318
===================================================================
Omnibus: 0.296 Durbin-Watson: 1.881
Prob(Omnibus): 0.862 Jarque-Bera (JB): 0.292
偏度：0.123 Prob(JB): 0.864
峰度：2.904 条件数 72.9
===================================================================
In the summary data, we can see that the `X3` variable is not significant since it has a p value of 0.66.
6. Since the third predictor variable is not significant, we eliminate this column and perform the regression again:
second_model = sm.OLS(Y, p_vars.loc[:, “const”:“X2”]).fit()
print(second_model.summary())
This results in a small increase in the goodness of fit statistics.
How it works...
Multilinear regression works in much the same way as simple linear regression. We follow the same procedure here as in the previous recipe, where we use the statsmodels package to fit a multilinear model to our data. Of course, there are some differences behind the scenes. The model we produce using multilinear regression is very similar in form to the simple linear model from the previous recipe. It has the following form:
Here, *Y* is the response variable, *X[i]* represents the predictor variables, *E* is the error term, and *β*[*i*] is the parameters to be computed. The same requirements are also necessary for this context: residuals must be independent and normally distributed with a mean of 0 and a common standard deviation.
In this recipe, we provided our predictor data as a Pandas `DataFrame` rather than a plain NumPy array. Notice that the names of the columns have been adopted in the summary data that we printed. Unlike the first recipe, *Using basic linear regression*, we included the constant column in this `DataFrame`, rather than using the `add_constant` utility from statsmodels.
In the output of the first regression, we can see that the model is a reasonably good fit with an adjusted *R²* value of 0.762, and is highly significant (we can see this by looking at the regression F statistic p value). However, looking closer at the individual parameters, we can see that both of the first two predictor values are significant, but the constant and the third predictor are less so. In particular, the third predictor parameter, `X3`, is not significantly different from 0 and has a p value of 0.66\. Given that our response data was constructed without using this variable, this shouldn't come as a surprise. In the final step of the analysis, we repeat the regression without the predictor variable, `X3`, which is a mild improvement to the fit.
Classifying using logarithmic regression
Logarithmic regression solves a different problem to ordinary linear regression. It is commonly used for classification problems where, typically, we wish to classify data into two distinct groups, according to a number of predictor variables. Underlying this technique is a transformation that's performed using logarithms. The original classification problem is transformed into a problem of constructing a model for the **log-odds***.* This model can be completed with simple linear regression. We apply the inverse transformation to the linear model, which leaves us with a model of the probability that the desired outcome will occur, given the predictor data. The transform we apply here is called the **logistic function**, which gives its name to the method. The probability we obtain can then be used in the classification problem we originally aimed to solve.
In this recipe, we will learn how to perform logistic regression and use this technique in classification problems.
Getting ready
For this recipe, we will need the NumPy package imported as `np`, the Matplotlib `pyplot`module imported as `plt`, the Pandas package imported as `pd`, and an instance of the NumPy default random number generator to be created using the following commands:
从 numpy.random 导入 default_rng
rng = default_rng(12345)
We also need several components from the `scikit-learn` package to perform logistic regression. These can be imported as follows:
从 sklearn.linear_model 导入 LogisticRegression
从 sklearn.metrics 导入 classification_report
How to do it...
Follow these steps to use logistic regression to solve a simple classification problem:
1. First, we need to create some sample data that we can use to demonstrate how to use logistic regression. We start by creating the predictor variables:
df = pd.DataFrame({
“var1”: np.concatenate([rng.normal(3.0, 1.5, size=50),
rng.normal(-4.0, 2.0, size=50)]),
“var2”: rng.uniform(size=100),
“var3”: np.concatenate([rng.normal(-2.0, 2.0, size=50),
rng.normal(1.5, 0.8, size=50)])
})
2. Now, we use two of our three predictor variables to create our response variable as a series of Boolean values:
分数 = 4.0 + df[“var1”] - df[“var3”]
Y = score >= 0
3. Next, we scatter plot the points, styled according to the response variable, of the `var3` data against the `var1` data, which are the variables used to construct the response variable:
fig1, ax1 = plt.subplots()
ax1.plot(df.loc[Y, “var1”], df.loc[Y, “var3”], “bo”, label="True
数据")
ax1.plot(df.loc[~Y, “var1”], df.loc[~Y, “var3”], “rx”, label="False
数据")
ax1.legend()
ax1.set_xlabel(“var1”)
ax1.set_ylabel(“var3”)
ax1.set_title(“Scatter plot of var3 against var1”)
The resulting plot can be seen in the following figure:
Figure 7.3: Scatter plot of the var3 data against var1, with classification marked
4. Next, we create a `LogisticRegression` object from the `scikit-learn` package and fit the model to our data:
model = LogisticRegression()
model.fit(df, Y)
5. Next, we prepare some extra data, different from what we used to fit the model, to test the accuracy of our model:
test_df = pd.DataFrame({
“var1”: np.concatenate([rng.normal(3.0, 1.5, size=50),
rng.normal(-4.0, 2.0, size=50)]),
“var2”: rng.uniform(size=100),
“var3”: np.concatenate([rng.normal(-2.0, 2.0, size=50),
rng.normal(1.5, 0.8, size=50)])
})
test_scores = 4.0 + test_df[“var1”] - test_df[“var3”]
test_Y = test_scores >= 0
6. Then, we generate predicted results based on our logistic regression model:
test_predicts = model.predict(test_df)
7. Finally, we use the `classification_report` utility from `scikit-learn` to print a summary of predicted classification against known response values to test the accuracy of the model. We print this summary to the Terminal:
print(classification_report(test_Y, test_predicts))
The report that's generated by this routine looks as follows:
precision recall f1-score support
False 1.00 1.00 1.00 18
True 1.00 1.00 1.00 32
accuracy 1.00 50
macro avg 1.00 1.00 1.00 50
weighted avg 1.00 1.00 1.00 50
How it works...
Logistic regression works by forming a linear model of the *log odds* ratio *(*or *logit*), which, for a single predictor variable, *x*, has the following form:
Here, *p*(*x*) represents the probability of a true outcome in response to the given the predictor, *x*. Rearranging this gives a variation of the logistic function for the probability:
The parameters for the log odds are estimated using a maximum likelihood method.
The `LogisticRegression` class from the `linear_model` module in `scikit-learn` is an implementation of logistic regression that is very easy to use. First, we create a new model instance of this class, with any custom parameters that we need, and then use the `fit` method on this object to fit (or train) the model to the sample data. Once this fitting is done, we can access the parameters that have been estimated using the `get_params` method.
The `predict` method on the fitted model allows us to pass in new (unseen) data and make predictions about the classification of each sample. We could also get the probability estimates that are actually given by the logistic function using the `predict_proba` method.
Once we have built a model for predicting the classification of data, we need to validate the model. This means we have to test the model with some previously unseen data and check whether it correctly classifies the new data. For this, we can use `classification_report`, which takes a new set of data and the predictions generated by the model and computes the proportion of the data that was correctly predicted by the model. This is the *precision* of the model.
The classification report we generated using the `scikit-learn` utility performs a comparison between the predicted results and the known response values. This is a common method for validating a model before using it to make actual predictions. In this recipe, we saw that the reported precision for each of the categories (`True` and `False`) was 1.00, indicating that the model performed perfectly in predicting the classification with this data. In practice, it is unlikely that the precision of a model will be 100%.
There's more...
There are lots of packages that offer tools for using logistic regression for classification problems. The statsmodels package has the `Logit` class for creating logistic regression models. We used the `scikit-learn` package in this recipe, which has a similar interface. `Scikit-learn` is a general-purpose machine learning library and has a variety of other tools for classification problems.
Modeling time series data with ARMA
Time series, as the name suggests, tracks a value over a sequence of distinct time intervals. They are particularly important in the finance industry, where stock values are tracked over time and used to make predictions – known as forecasting – of the value at some future time. Good predictions coming from such data can be used to make better investments. Time series also appear in many other common situations, such as weather monitoring, medicine, and any places where data is derived from sensors over time.
Time series, unlike other types of data, do not usually have independent data points. This means that the methods that we use for modeling independent data will not be particularly effective. Thus, we need to use alternative techniques to model data with this property. There are two ways in which a value in a time series can depend on previous values. The first is where there is a direct relationship between the value and one or more previous values. This is the *autocorrelation* property and is modeled by an *autoregressive* model. The second is where the noise that's added to the value depends on one or more previous noise terms. This is modeled by a *moving average* model. The number of terms involved in either of these models is called the *order* of the model.
In this recipe, we will learn how to create a model for stationary time series data with ARMA terms.
Getting ready
For this recipe, we need the Matplotlib `pyplot` module imported as `plt` and the statsmodels package `api` module imported as `sm`. We also need to import the `generate_sample_data` routine from the `tsdata` package from this book's repository, which uses NumPy and Pandas to generate sample data for analysis:
from tsdata import generate_sample_data
How to do it...
Follow these steps to create an autoregressive moving average model for stationary time series data:
1. First, we need to generate the sample data that we will analyze:
sample_ts, _ = generate_sample_data()
2. As always, the first step in the analysis is to produce a plot of the data so that we can visually identify any structure:
ts_fig, ts_ax = plt.subplots()
sample_ts.plot(ax=ts_ax, label=“Observed”)
ts_ax.set_title(“Time series data”)
ts_ax.set_xlabel(“Date”)
ts_ax.set_ylabel(“Value”)
The resulting plot can be seen in the following figure. Here, we can see that there doesn't appear to be an underlying trend, which means that the data is likely to be stationary:
Figure 7.4: Plot of the time series data that we will analyze. There doesn't appear to be a trend in this data
3. Next, we compute the augmented Dickey-Fuller test. The null hypothesis is that the time series is not stationary:
adf_results = sm.tsa.adfuller(sample_ts)
adf_pvalue = adf_results[1]
print(“Augmented Dickey-Fuller test:\nP-value:”, adf_pvalue)
The reported p value is 0.000376 in this case, so we reject the null hypothesis and conclude that the series is stationary.
4. Next, we need to determine the order of the model that we should fit. For this, we'll plot the **autocorrelation function** (**ACF**) and the **partial autocorrelation function**(**PACF**) for the time series:
ap_fig, (acf_ax, pacf_ax) = plt.subplots(2, 1, sharex=True,
tight_layout=True)
sm.graphics.tsa.plot_acf(sample_ts, ax=acf_ax,
title=“Observed autocorrelation”)
sm.graphics.tsa.plot_pacf(sample_ts, ax=pacf_ax,
title=“Observed partial autocorrelation”)
pacf_ax.set_xlabel(“Lags”)
pacf_ax.set_ylabel(“Value”)
acf_ax.set_ylabel(“Value”)
The plots of the ACF and PACF for our time series can be seen in the following figure. These plots suggest the existence of both autoregressive and moving average processes:
Figure 7.5: ACF and PACF for the sample time series data
5. Next, we create an ARMA model for the data, using the `ARMA` class from statsmodels, `tsa` module. This model will have an order 1 AR and an order 1 MA:
arma_model = sm.tsa.ARMA(sample_ts, order=(1, 1))
6. Now, we fit the model to the data and get the resulting model. We print a summary of these results to the Terminal:
arma_results = arma_model.fit()
print(arma_results.summary())
The summary data given for the fitted model is as follows:
ARMA Model Results
===================================================================
Dep. Variable: y No. Observations: 366
Model: ARMA(1, 1) Log Likelihood -513.038
方法：css-mle 创新的标准偏差 0.982
日期：2020 年 5 月 1 日 AIC 1034.077
时间：12:40:00 BIC 1049.687
Sample: 01-01-2020 HQIC 1040.280
- 12-31-2020
===================================================================
coef std err z P>|z| [0.025 0.975]
const -0.0242 0.143 -0.169 0.866 -0.305 0.256
ar.L1.y 0.8292 0.057 14.562 0.000 0.718 0.941
ma.L1.y -0.5189 0.090 -5.792 0.000 -0.695 -0.343
Roots
===================================================================
实部 虚部 模数 频率
AR.1 1.2059 +0.0000j 1.2059 0.0000
MA.1 1.9271 +0.0000j 1.9271 0.0000
Here, we can see that both of the estimated parameters for the AR and MA components are significantly different from 0\. This is because the value in the `P >|z|` column is 0 to 3 decimal places.
7. Next, we need to verify that there is no additional structure remaining in the residuals (error) of the predictions from our model. For this, we plot the ACF and PACF of the residuals:
residuals = arma_results.resid
rap_fig, (racf_ax, rpacf_ax) = plt.subplots(2, 1,
sharex=True, tight_layout=True)
sm.graphics.tsa.plot_acf(residuals, ax=racf_ax,
title=“Residual autocorrelation”)
sm.graphics.tsa.plot_pacf(residuals, ax=rpacf_ax,
标题=“残差部分自相关”)
rpacf_ax.set_xlabel(“滞后”)
rpacf_ax.set_ylabel(“值”)
racf_ax.set_ylabel(“值”)
The ACF and PACF of the residuals can be seen in the following figure. Here, we can see that there are no significant spikes at lags other than 0, so we conclude that there is no structure remaining in the residuals:
Figure 7.6: ACF and PACF for the residuals from our model
8. Now that we have verified that our model is not missing any structure, we plot the values that are fitted to each data point on top of the actual time series data to see whether the model is a good fit for the data. We plot this model in the plot we created in *step 2*:
fitted = arma_results.fittedvalues
fitted.plot(c=“r”, ax=ts_ax, 标签=“拟合”)
ts_ax.legend()
The updated plot can be seen in the following figure:
Figure 7.7: Plot of the fitted time series data over the observed time series data
The fitted values give a reasonable approximation of the behavior of the time series, but reduce the noise from the underlying structure.
How it works...
A time series is stationary if it does not have a trend. They usually have a tendency to move in one direction rather than another. Stationary processes are important because we can usually remove the trend from an arbitrary time series and model the underlying stationary series. The ARMA model that we used in this recipe is a basic means of modeling the behavior of stationary time series. The two parts of an ARMA model are the autoregressive and moving average parts, which model the dependence of the terms and noise, respectively, on previous terms and noise.
An order 1 autoregressive model has the following form:
Here, *φ[i]* represents the parameters and *ε[t]* is the noise at a given step. The noise is usually assumed to be normally distributed with a mean of 0 and a standard deviation that is constant across all the time steps. The *Y[t]* value represents the value of the time series at the time step, *t*. In this model, each value depends on the previous value, though it can also depend on some constants and some noise. The model will give rise to a stationary time series precisely when the *φ[1]* parameter lies strictly between -1 and 1.
An order 1 moving average model is very similar to an autoregressive model and is given by the following equation:
Here, the variants of *θ[i]* are parameters. Putting these two models together gives us an ARMA(1, 1) model, which has the following form:
In general, we can have an ARMA(p, q) model that has an order *p* AR component and an order q MA component. We usually refer to the quantities, *p* and *q*, as the orders of the model.
Determining the orders of the AR and MA components is the most tricky aspect of constructing an ARMA model. The ACF and PACF give some information toward this, but even then, it can be quite difficult. For example, an autoregressive process will show some kind of decay or oscillating pattern on the ACF as lag increases, and a small number of peaks on the PACF and values that are not significantly different from 0 beyond that. The number of peaks that appear on the PAF plot can be taken as the order of the process. For a moving average process, the reverse is true. There are usually a small number of significant peaks on the ACF plot, and a decay or oscillating pattern on the PACF plot. Of course, sometimes, this isn't obvious.
In this recipe, we plotted the ACF and PACF for our sample time series data. In the autocorrelation plot in *Figure 7.5* (top), we can see that the peaks decay rapidly until they lie within the confidence interval of zero (meaning they are not significant). This suggests the presence of an autoregressive component. On the partial autocorrelation plot in *Figure 7.5* (bottom), we can see that there are only two peaks that can be considered not zero, which suggests an autoregressive process of order 1 or 2\. You should try to keep the order of the model as small as possible. Due to this, we chose an order 1 autoregressive component. With this assumption, the second peak on the partial autocorrelation plot is indicative of decay (rather than an isolated peak), which suggests the presence of a moving average process. To keep the model simple, we try an order 1 moving average process. This is how the model that we used in this recipe was decided on. Notice that this is not an exact process, and you might have decided differently.
We use the augmented Dickey-Fuller test to test the likelihood that the time series that we have observed is stationary. This is a statistical test, such as those seen in Chapter 6, *Working with Data and Statistics*, that generates a test statistic from the data. This test statistic, in turn, determines a p-value that is used to determine whether to accept or reject the null hypothesis. For this test, the null hypothesis is that a unit root is present in the time series that's been sampled. The alternative hypothesis – the one we are really interested in – is that the observed time series is (trend) stationary. If the p-value is sufficiently small, then we can conclude with the specified confidence that the observed time series is stationary. In this recipe, the p-value was 0.000 to 3 decimal places, which indicates a strong likelihood that the series is stationary. Stationarity is an essential assumption for using the ARMA model for the data.
Once we have determined that the series is stationary, and also decided on the orders of the model, we have to fit the model to the sample data that we have. The parameters of the model are estimated using a maximum likelihood estimator. In this recipe, the learning of the parameters is done using the `fit` method, in *step 6*.
The statsmodels package provides various tools for working with time series, including utilities for calculating – and plotting –ACF and PACF of time series data, various test statistics, and creating ARMA models for time series. There are also some tools for automatically estimating the order of the model.
We can use the **Akaike information criterion** (**AIC**), **Bayesian information criterion** (**BIC**), and **Hannan-Quinn Information Criterion** (**HQIC**) quantities to compare this model to other models to see which model best describes the data. A smaller value is better in each case.
When using ARMA to model time series data, as in all kinds of mathematical modeling tasks, it is best to pick the simplest model that describes the data to the extent that is needed. For ARMA models, this usually means picking the smallest order model that describes the structure of the observed data.
There's more...
Finding the best combination of orders for an ARMA model can be quite difficult. Often, the best way to fit a model is to test multiple different configurations and pick the order that produces the best fit. For example, we could have tried ARMA(0, 1) or ARMA(1, 0) in this recipe, and compared it to the ARMA(1, 1) model we used to see which produced the best fit by considering the **Akaike Information Criteria** (**AIC**) statistic reported in the summary. In fact, if we build these models, we will see that the AIC value for ARMA(1, 1) – the model we used in this recipe – is the "best" of these three models.
Forecasting from time series data using ARIMA
In the previous recipe, we generated a model for a stationary time series using an ARMA model, which consists of an **autoregressive** (**AR**) component and an **m****oving average** (**MA**) component. Unfortunately, this model cannot accommodate time series that have some underlying trend; that is, they are not stationary time series. We can often get around this by *differencing* the observed time series one or more times until we obtain a stationary time series that can be modeled using ARMA. The incorporation of differencing into an ARMA model is called an ARIMA model, which stands for **Autoregressive** (**AR**) **Integrated** (**I**) **Moving Average** (**MA**).
Differencing is the process of computing the difference of consecutive terms in a sequence of data. So, applying first-order differencing amounts to subtracting the value at the current step from the value at the next step (*t[i+1] - t[i]*). This has the effect of removing the underlying upward or downward linear trend from the data. This helps to reduce an arbitrary time series to a stationary time series that can be modeled using ARMA. Higher-order differencing can remove higher-order trends to achieve similar effects.
An ARIMA model has three parameters, usually labeled *p*, *d*, and *q*. The *p* and *q* order parameters are the order of the autoregressive component and the moving average component, respectively, just as they are for the ARMA model. The third order parameter, *d*, is the order of differencing to be applied. An ARIMA model with these orders is usually written as ARIMA (*p*, *d*, *q*). Of course, we will need to determine what order differencing should be included before we start fitting the model.
In this recipe, we will learn how to fit an ARIMA model to a non-stationary time series and use this model to make forecasts about future values.
Getting ready
For this recipe, we will need the NumPy package imported as `np`, the Pandas package imported as `pd`, the Matplotlib `pyplot` module as `plt`, and the statsmodels `api` module imported as `sm`. We will also need the utility for creating sample time series data from the `tsdata` module, which is included in this book's repository:
从 tsdata 导入生成样本数据
How to do it...
The following steps show you how to construct an ARIMA model for time series data and use this model to make forecasts:
1. First, we load the sample data using the `generate_sample_data` routine:
sample_ts，test_ts = generate_sample_data(trend=0.2, undiff=True)
2. As usual, the next step is to plot the time series so that we can visually identify the trend of the data:
ts_fig, ts_ax = plt.subplots(tight_layout=True)
sample_ts.plot(ax=ts_ax, c=“b”, 标签=“观察到的”)
ts_ax.set_title(“训练时间序列数据”)
ts_ax.set_xlabel(“日期”)
ts_ax.set_ylabel(“值”)
The resulting plot can be seen in the following figure. As we can see, there is a clear upward trend in the data, so the time series is certainly not stationary:
Figure 7.8: Plot of the sample time series. There is an obvious positive trend in the data.
3. Next, we difference the series to see if one level of differencing is sufficient to remove the trend:
差分=sample_ts.diff().dropna()
4. Now, we plot the ACF and PACF for the differenced time series:
ap_fig, (acf_ax, pacf_ax) = plt.subplots(1, 2,
tight_layout=True, sharex=True)
sm.graphics.tsa.plot_acf(diffs, ax=acf_ax)
sm.graphics.tsa.plot_pacf(diffs, ax=pacf_ax)
acf_ax.set_ylabel(“值”)
pacf_ax.set_xlabel(“滞后”)
pacf_ax.set_ylabel(“值”)
The ACF and PACF can be seen in the following figure. We can see that there does not appear to be any trends left in the data and that there appears to be both an autoregressive component and a moving average component:
Figure 7.9: ACF and PACF for the differenced time series
5. Now, we construct the ARIMA model with order 1 differencing, an autoregressive component, and a moving average component. We fit this to the observed time series and print a summary of the model:
model = sm.tsa.ARIMA(sample_ts, order=(1,1,1))
fitted = model.fit(trend=“c”)
print(fitted.summary())
The summary information that's printed looks as follows:
ARIMA 模型结果
==================================================================
因变量：D.y 观察次数：365
模型：ARIMA(1, 1, 1) 对数似然-512.905
方法：css-mle 创新的标准差 0.986
日期：星期六，2020 年 5 月 2 日 AIC 1033.810
时间：14:47:25 BIC 1049.409
样本：2020 年 01 月 02 日 HQIC 1040.009
- 12-31-2020
==================================================================
coef std err z P>|z| [0.025 0.975]
const 0.9548 0.148 6.464 0.000 0.665 1.244
ar.L1.D.y 0.8342 0.056 14.992 0.000 0.725 0.943
ma.L1.D.y -0.5204 0.088 -5.903 0.000 -0.693 -0.348
根
==================================================================
实际 虚数 模数 频率
AR.1 1.1987 +0.0000j 1.1987 0.0000
MA.1 1.9216 +0.0000j 1.9216 0.0000
Here, we can see that all three of our estimated coefficients are significantly different from 0 due to the fact that all three have 0 to 3 decimal places in the `P>|z|` column.
6. Now, we can use the `forecast` method to generate predictions of future values. This also returns the standard error and confidence intervals for predictions:
forecast, std_err, fc_ci = fitted.forecast(steps=50)
forecast_dates = pd.date_range(“2021-01-01”, periods=50)
forecast = pd.Series(forecast, index=forecast_dates)
7. Next, we plot the forecast values and their confidence intervals on the figure containing the time series data:
forecast.plot(ax=ts_ax, c=“g”, 标签=“预测”)
ts_ax.fill_between(forecast_dates, fc_ci[:, 0], fc_ci[:, 1],
颜色=“r”, alpha=0.4)
8. Finally, we add the actual future values to generate, along with the sample in *step 1*, to the plot (it might be easier if you repeat the plot commands from *step 1* to regenerate the whole plot here):
test_ts.plot(ax=ts_ax, c=“k”, 标签=“实际”)
ts_ax.legend()
The final plot containing the time series with the forecast and the actual future values can be seen in the following figure:
Figure 7.10: Plot of the sample time series with forecast values and actual future values for comparison
Here, we can see that the actual future values are within the confidence interval for the forecast values.
How it works...
The ARIMA model – with orders *p*, *d*, and *q –* is simply an ARMA (*p*, *q*) model that's applied to a time series. This is obtained by applying differencing of order *d* to the original time series data. It is a fairly simple way to generate a model for time series data. The statsmodels `ARIMA` class handles the creation of a model, while the `fit` method fits this model to the data. We passed the `trend="c"` keyword argument because we know, from *Figure 7.9*, that the time series has a constant trend.
The model is fit to the data using a maximum likelihood method and the final estimates for the parameters – in this case, one parameter for the autoregressive component, one for the moving average component, the constant trend parameter, and the variance of the noise. These parameters are reported in the summary. From this output, we can see that the estimates for the AR coefficient (0.8342) and the MA constant (-0.5204) are very good approximations of the true estimates that were used to generate the data, which were 0.8 for the AR coefficient and -0.5 for the MA coefficient. These parameters are set in the `generate_sample_data` routine from the `tsdata.py` file in the code repository for this chapter. This generates the sample data in *step 1*. You might have noticed that the constant parameter (0.9548) is not 0.2, as specified in the `generate_sample_data` call in *step 1*. In fact, it is not so far from the actual drift of the time series.
The `forecast` method on the fitted model (the output of the `fit` method) uses the model to make predictions about the value after a given number of steps. In this recipe, we forecast for up to 50 time steps beyond the range of the sample time series. The output of the `forecast` method is a tuple containing the forecast values, the standard error for the forecasts, and the confidence interval (by default, 95% confidence) for the forecasts. Since we provided the time series as a Pandas series, these are returned as `Series` objects (the confidence interval is a `DataFrame`).
When you construct an ARIMA model for time series data, you need to make sure you use the smallest order differencing that removes the underlying trend. Applying more differencing than is necessary is called *overdifferencing* and can lead to problems with the model.
Forecasting seasonal data using ARIMA
Time series often display periodic behavior so that peaks or dips in the value appear at regular intervals. This behavior is called *seasonality* in the analysis of time series. The methods we have used to far in this chapter to model time series data obviously do not account for seasonality. Fortunately, it is relatively easy to adapt the standard ARIMA model to incorporate seasonality, resulting in what is sometimes called a SARIMA model.
In this recipe, we will learn how to model time series data that includes seasonal behavior and use this model to produce forecasts.
Getting ready
For this recipe, we will need the NumPy package imported as `np`, the Pandas package imported as `pd`, the Matplotlib `pyplot`module as `plt`, and the statsmodels `api`module imported as `sm`. We will also need the utility for creating sample time series data from the `tsdata`module, which is included in this book's repository:
从 tsdata 导入生成样本数据
How to do it...
Follow these steps to produce a seasonal ARIMA model for sample time series data and use this model to produce forecasts:
1. First, we use the `generate_sample_data` routine to generate a sample time series to analyze:
sample_ts，test_ts = generate_sample_data(undiff=True，
季节性=True)
2. As usual, our first step is to visually inspect the data by producing a plot of the sample time series:
ts_fig, ts_ax = plt.subplots(tight_layout=True)
sample_ts.plot(ax=ts_ax, 标题=“时间序列”, 标签=“观察到的”)
ts_ax.set_xlabel(“日期”)
ts_ax.set_ylabel(“值”)
The plot of the sample time series data can be seen in the following figure. Here, we can see that there seem to be periodic peaks in the data:
Figure 7.11: Plot of the sample time series data
3. Next, we plot the ACF and PACF for the sample time series:
ap_fig，(acf_ax，pacf_ax) = plt.subplots(2, 1，
sharex=True, tight_layout=True)
sm.graphics.tsa.plot_acf(sample_ts, ax=acf_ax)
sm.graphics.tsa.plot_pacf(sample_ts, ax=pacf_ax)
pacf_ax.set_xlabel(“滞后”)
acf_ax.set_ylabel(“值”)
pacf_ax.set_ylabel(“值”)
The ACF and PACF for the sample time series can be seen in the following figure:
Figure 7.12: ACF and PACF for the sample time series
These plots possibly indicate the existence of autoregressive components, but also a significant spike on the PACF with lag 7.
4. Next, we difference the time series and produce plots of the ACF and PACF for the differenced series. This should make the order of the model clearer:
差分=sample_ts.diff().dropna()
dap_fig, (dacf_ax, dpacf_ax) = plt.subplots(2, 1, sharex=True,
tight_layout=True)
sm.graphics.tsa.plot_acf(diffs, ax=dacf_ax,
标题=“差分 ACF”)
sm.graphics.tsa.plot_pacf(diffs, ax=dpacf_ax,
标题=“差分 PACF”)
dpacf_ax.set_xlabel(“滞后”)
dacf_ax.set_ylabel(“值”)
dpacf_ax.set_ylabel(“值”)
The ACF and PACF for the differenced time series can be seen in the following figure. We can see that there is definitely a seasonal component with lag 7:
Figure 7.13: Plot of the ACF and PACF for the differenced time series
5. Now, we need to create a `SARIMAX` object that holds the model, with ARIMA order `(1, 1, 1)` and seasonal ARIMA order `(1, 0, 0, 7)`. We fit this model to the sample time series and print summary statistics. We plot the predicted values on top of the time series data:
model = sm.tsa.SARIMAX(sample_ts, order=(1, 1, 1),
季节性顺序=(1, 0, 0, 7))
fitted_seasonal = model.fit()
print(fitted_seasonal.summary())
fitted_seasonal.fittedvalues.plot(ax=ts_ax, c=“r”,
标签=“预测”)
The summary statistics that are printed to the Terminal look as follows:
SARIMAX 结果
===================================================================
因变量：y 观察次数：366
模型：SARIMAX(1, 1, 1)x(1, 0, [], 7) 对数似然-509.941
日期：星期一，2020 年 5 月 4 日 AIC 1027.881
时间：18:03:27 BIC 1043.481
样本：2020 年 01 月 01 日 HQIC 1034.081
- 12-31-2020
协方差类型： opg
===================================================================
coef std err z P>|z| [0.025 0.975]
ar.L1 0.7939 0.065 12.136 0.000 0.666 0.922
ma.L1 -0.4544 0.095 -4.793 0.000 -0.640 -0.269
ar.S.L7 0.7764 0.034 22.951 0.000 0.710 0.843
sigma2 0.9388 0.073 12.783 0.000 0.795 1.083
===================================================================
Ljung-Box (Q): 31.89 Jarque-Bera (JB): 0.47
Prob(Q): 0.82 Prob(JB): 0.79
异方差性（H）: 1.15 偏度: -0.03
Prob(H) (双侧): 0.43 峰度: 2.84
===================================================================
警告:
[1] 使用外积计算的协方差矩阵
梯度的数量（复杂步骤）。
6. This model appears to be a reasonable fit, so we move ahead and forecast `50` time steps into the future:
forecast_result = fitted_seasonal.get_forecast(steps=50)
forecast_index = pd.date_range(“2021-01-01”, periods=50)
预测 = 预测结果.预测均值
7. Finally, we add the forecast values to the plot of the sample time series, along with the confidence interval for these forecasts:
forecast.plot(ax=ts_ax, c=“g”, label=“预测”)
conf = forecast_result.conf_int()
ts_ax.fill_between(forecast_index, conf[“lower y”],
conf[“upper y”], color=“r”, alpha=0.4)
test_ts.plot(ax=ts_ax, color=“k”, label=“实际未来”)
ts_ax.legend()
The final plot of the time series, along with the predictions and the confidence interval for the forecasts, can be seen in the following figure:
Figure 7.14: Plot of the sample time series, along with the forecasts and confidence interval
How it works...
Adjusting an ARIMA model to incorporate seasonality is a relatively simple task. A seasonal component is similar to an autoregressive component, where the lag starts at some number larger than 1\. In this recipe, the time series exhibits seasonality with period 7 (weekly), which means that the model is approximately given by the following equation:
Here *φ[1]* and *Φ**[1]**are the parameters and *ε[t]* is the noise at time step *t*. The standard ARIMA model is easily adapted to include this additional lag term.* *The SARIMA model incorporates this additional seasonality into the ARIMA model. It has four additional order terms on top of the three for the underlying ARIMA model. These four additional parameters are the seasonal AR, differencing, and MA components, along with the period of the seasonality. In this recipe, we took the seasonal AR to be order 1, with no seasonal differencing or MA components (order 0), and a seasonal period of 7\. This gives us the additional parameters (1, 0, 0, 7) that we used in *step 5* of this recipe.
Seasonality is clearly important in modeling time series data that is measured over a period of time covering days, months, or years. It usually incorporates some kind of seasonal component based on the time frame that they occupy. For example, a time series of national power consumption measured hourly over several days would probably have a 24-hour seasonal component since power consumption will likely fall during the night hours.
Long-term seasonal patterns might be hidden if the time series data that you are analyzing does not cover a sufficiently large time period for the pattern to emerge. The same is true for trends in the data. This can lead to some interesting problems when trying to produce long-term forecasts from a relatively short period represented by observed data.
The `SARIMAX` class from the statsmodels package provides the means of modeling time series data using a seasonal ARIMA model. In fact, it can also model external factors that have an additional effect on the model, sometimes called *exogenous regressors*. (We will not cover these here.) This class works much like the `ARMA` and `ARIMA` classes that we used in the previous recipes. First, we create the model object by providing the data and orders for both the ARIMA process and the seasonal process, and then use the `fit` method on this object to create a fitted model object. We use the `get_forecasts` method to generate an object holding the forecasts and confidence interval data that we can then plot, thus producing the *Figure 7.14*.
## There's more...
There is a small difference in the interface between the `SARIMAX` class used in this recipe and the `ARIMA` class used in the previous recipe. At the time of writing, the statsmodels package (v0.11) includes a second `ARIMA` class that builds on top of the `SARIMAX` class, thus providing the same interface. However, at the time of writing, this new `ARIMA` class does not offer the same functionality as that used in this recipe.
# Using Prophet to model time series data
The tools we have seen so far for modeling time series data are very general and flexible methods, but they require some knowledge of time series analysis in order to be set up. The analysis needed to construct a good model that can be used to make reasonable predictions into the future can be intensive and time-consuming, and may not be viable for your application. The Prophet library is designed to automatically model time series data quickly, without the need for input from the user, and make predictions into the future.
In this recipe, we will learn how to use Prophet to produce forecasts from a sample time series.
## Getting ready
For this recipe, we will need the Pandas package imported as `pd`, the Matplotlib `pyplot` package imported as `plt`, and the `Prophet` object from the Prophet library, which can be imported using the following command:
from fbprophet import Prophet
We also need to import the `generate_sample_data` routine from the `tsdata` module, which is included in the code repository for this book:
from tsdata import generate_sample_data
## How to do it...
The following steps show you how to use the Prophet package to generate forecasts for a sample time series:
1. First, we use `generate_sample_data` to generate the sample time series data:
sample_ts, test_ts = generate_sample_data(undiff=True, trend=0.2)
2. We need to convert the sample data into a `DataFrame` that Prophet expects:
df_for_prophet = pd.DataFrame({
“ds”: sample_ts.index, # dates
“y”: sample_ts.values # values
})
3. Next, we make a model using the `Prophet` class and fit it to the sample time series:
model = Prophet()
model.fit(df_for_prophet)
4. Now, we create a new `DataFrame` that contains the time intervals for the original time series, plus the additional periods for the forecasts:
forecast_df = model.make_future_dataframe(periods=50)
5. Then, we use the `predict` method to produce the forecasts along the time periods we just created:
forecast = model.predict(forecast_df)
6. Finally, we plot the predictions on top of the sample time series data, along with the confidence interval and the true future values:
fig, ax = plt.subplots(tight_layout=True)
sample_ts.plot(ax=ax, label=“观察到的”, title=“预测”)
forecast.plot(x=“ds”, y=“yhat”, ax=ax, c=“r”,
label=“预测”)
ax.fill_between(forecast[“ds”].values, forecast[“yhat_lower”].values,
forecast[“yhat_upper”].values, color=“r”, alpha=0.4)
test_ts.plot(ax=ax, c=“k”, label=“未来”)
ax.legend()
ax.set_xlabel(“日期”)
ax.set_ylabel(“值”)
The plot of the time series, along with forecasts, can be seen in the following figure:
Figure 7.15: Plot of sample time series data, along with forecasts and a confidence interval
## How it works...
Prophet is a package that's used to automatically produce models for time series data based on sample data, with little extra input needed from the user. In practice, it is very easy to use; we just need to create an instance of the `Prophet` class, call the `fit` method, and then we are ready to produce forecasts and understand our data using the model.
The `Prophet` class expects the data in a specific format: a `DataFrame` with columns named `ds` for the date/time index, and `y` for the response data (the time series values). This `DataFrame` should have integer indices. Once the model has been fit, we use `make_future_dataframe` to create a `DataFrame` in the correct format, with appropriate date intervals, and with additional rows for future time intervals. The `predict` method then takes this `DataFrame` and produces values using the model to populate these time intervals with predicted values. We also get other information, such as the confidence intervals, in this forecast's `DataFrame`.
## There's more...
Prophet does a fairly good job of modeling time series data without any input from the user. However, the model can be customized using various methods from the `Prophet` class. For example, we could provide information about the seasonality of the data using the `add_seasonality` method of the `Prophet` class, prior to fitting the model.
There are alternative packages for automatically generating models for time series data. For example, popular machine learning libraries such as TensorFlow can be used to model time series data.
# Further reading
A good textbook on regression in statistics is the book *Probability and Statistics* by Mendenhall, Beaver, and Beaver, as mentioned in Chapter 6, *Working with Data and Statistics*. The following books provide a good introduction to classification and regression in modern data science:
* *James, G. and Witten, D., 2013\. An Introduction To Statistical Learning: With Applications In R. New York: Springer.*
* *Müller, A. and Guido, S., 2016\. Introduction To Machine Learning With Python. Sebastopol: O'Reilly Media.*
A good introduction to time series analysis can be found in the following book:
* *Cryer, J. and Chan, K., 2008\. Time Series Analysis. New York: Springer.**
```**
# 第九章：几何问题
本章描述了关于二维几何的几个问题的解决方案。几何是数学的一个分支，涉及点、线和其他图形（形状）的特征，这些图形之间的相互作用以及这些图形的变换。在本章中，我们将重点关注二维图形的特征以及这些对象之间的相互作用。
在 Python 中处理几何对象时，我们必须克服几个问题。最大的障碍是表示问题。大多数几何对象占据二维平面中的一个区域，因此不可能存储区域内的每个点。相反，我们必须找到一种更紧凑的方式来表示可以存储为相对较少的点的区域。例如，我们可以存储沿对象边界的一些点，从而可以重建边界和对象本身。此外，我们将几何问题重新表述为可以使用代表性数据回答的问题。
第二个最大的问题是将纯几何问题转化为可以使用软件理解和解决的形式。这可能相对简单-例如，找到两条直线相交的点是解决矩阵方程的问题-或者可能非常复杂，这取决于所提出的问题类型。解决这些问题的常见技术是使用更简单的对象表示所讨论的图形，并使用每个简单对象解决（希望）更容易的问题。然后，这应该给我们一个关于原始问题的解决方案的想法。
我们将首先向您展示如何可视化二维形状，然后学习如何确定一个点是否包含在另一个图形中。然后，我们将继续查看边缘检测、三角剖分和寻找凸包。最后，我们将通过构造贝塞尔曲线来结束本章。
本章包括以下教程：
+ 可视化二维几何形状
+ 寻找内部点
+ 在图像中查找边缘
+ 对平面图形进行三角剖分
+ 计算凸包
+ 构造贝塞尔曲线
让我们开始吧！
# 技术要求
对于本章，我们将像往常一样需要`numpy`包和`matplotlib`包。我们还需要 Shapely 包和`scikit-image`包，可以使用您喜欢的软件包管理器（如`pip`）进行安装：
```py
python3.8 -m pip install numpy matplotlib shapely scikit-image
该章节的代码可以在 GitHub 存储库的Chapter 08文件夹中找到：github.com/PacktPublishing/Applying-Math-with-Python/tree/master/Chapter%2008。
查看以下视频以查看代码实际操作：bit.ly/3hpeKEF。
可视化二维几何形状
本章的重点是二维几何，因此我们的第一个任务是学习如何可视化二维几何图形。这里提到的一些技术和工具可能适用于三维几何图形，但通常需要更专门的软件包和工具。
几何图形，至少在本书的上下文中，是指边界是一组线和曲线的任何点、线、曲线或封闭区域（包括边界）的集合。简单的例子包括点和线（显然）、矩形、多边形和圆。
在本教程中，我们将学习如何使用 Matplotlib 可视化几何图形。
准备工作
对于本教程，我们需要导入 NumPy 包作为np，导入 Matplotlib pyplot模块作为plt。我们还需要从 Matplotlib patches模块导入Circle类和从 Matplotlib collections模块导入PatchCollection类。可以使用以下命令完成这些操作：
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from matplotlib.collections import PatchCollection
我们还需要本章的代码库中的swisscheese-grid-10411.csv数据文件。
如何做…
以下步骤向您展示了如何可视化一个二维几何图形：
- 首先，我们从本书的代码库中的
swisscheese-grid-10411.csv文件加载数据：
data = np.loadtxt("swisscheese-grid-10411.csv")
- 我们创建一个代表绘图区域的新补丁对象。这将是一个圆（圆盘），其中心在原点，半径为
1。我们创建一个新的轴集，并将这个补丁添加到其中：
fig, ax = plt.subplots()
outer = Circle((0.0, 0.0), 1.0, zorder=0, fc="k")
ax.add_patch(outer)
- 接下来，我们从步骤 1中加载的数据创建一个
PatchCollection对象，其中包含了一些其他圆的中心和半径。然后我们将这个PatchCollection添加到步骤 2中创建的轴上：
col = PatchCollection(
(Circle((x, y), r) for x, y, r in data),
facecolor="white", zorder=1, linewidth=0.2,
ls="-", ec="k"
)
ax.add_collection(col)
- 最后，我们设置*x-和y-*轴的范围，以便整个图像都能显示出来，然后关闭轴线：
ax.set_xlim((-1.1, 1.1))
ax.set_ylim((-1.1, 1.1))
ax.set_axis_off()
结果图像是一个瑞士奶酪，如下所示：
图 8.1：瑞士奶酪的绘图
它是如何工作的…
这个食谱的关键是Circle和PatchCollection对象，它们代表了 Matplotlib Axes上的绘图区域的区域。在这种情况下，我们创建了一个大的圆形补丁，它位于原点，半径为1，具有黑色的面颜色，并使用zorder=0将其放在其他补丁的后面。这个补丁被添加到Axes对象中使用add_patch方法。
下一步是创建一个对象，它将呈现从 CSV 文件中加载的数据所代表的圆。这些数据包括中心（x，y）和半径r的值，用于表示个别圆的中心和半径（总共 10,411 个）。PatchCollection对象将一系列补丁组合成一个单一的对象，可以添加到Axes对象中。在这里，我们为我们的数据中的每一行添加了一个Circle，然后使用add_collection方法将其添加到Axes对象中。请注意，我们已经将面颜色应用到整个集合，而不是每个单独的Circle成员。我们将面颜色设置为白色（使用facecolor="w"参数），边缘颜色设置为黑色（使用ec="k"），边缘线宽设置为 0.2（使用linewidth=0.2），边缘样式设置为连续线。所有这些放在一起，就得到了我们的图像。
我们在这里创建的图像被称为“瑞士奶酪”。这些最初是由爱丽丝·罗斯在 1938 年在有理逼近理论中使用的；随后它们被重新发现，并且类似的构造自那时以来已经被多次使用。我们使用这个例子是因为它由一个大的个体部分和一个大量的较小的个体部分组成。罗斯的瑞士奶酪是平面上具有正面积但没有拓扑内部的一个集合的例子。（这样一个集合甚至能存在是相当惊人的！）更重要的是，有一些连续函数在这个瑞士奶酪上是不能被有理函数逼近的。这个特性使得类似的构造在均匀 代数理论中非常有用。
Circle类是更一般的Patch类的子类。还有许多其他Patch类，代表不同的平面图形，比如Polygon和PathPatch，它们代表了由路径（曲线或曲线集合）所界定的区域。这些可以用来生成可以在 Matplotlib 图中呈现的复杂补丁。集合可以用来同时应用设置到多个补丁对象上，这在本例中特别有用，因为你有大量的对象都将以相同的样式呈现。
还有更多…
Matplotlib 中有许多不同的补丁类型。在本教程中，我们使用了Circle补丁类，它表示坐标轴上的圆形区域。还有Polygon补丁类，它表示多边形（规则或其他）。还有PatchPath对象，它们是由不一定由直线段组成的曲线包围的区域。这类似于许多矢量图形软件包中可以构造阴影区域的方式。
除了 Matplotlib 中的单个补丁类型外，还有许多集合类型，它们将许多补丁聚集在一起，以便作为单个对象使用。在本教程中，我们使用了PatchCollection类来收集大量的Circle补丁。还有更多专门的补丁集合，可以用来自动生成这些内部补丁，而不是我们自己生成它们。
另请参阅
在数学中，可以在以下传记文章中找到关于瑞士奶酪的更详细的历史：Daepp,U., Gauthier, P., Gorkin, P. and Schmieder, G., 2005. Alice in Switzerland: The life and mathematics of Alice Roth. The Mathematical Intelligencer, 27(1), pp.41-54。
查找内部点
在编程环境中处理二维图形的一个问题是，您不可能存储所有位于图形内的点。相反，我们通常存储表示图形的远少于的点。在大多数情况下，这将是一些点（通过线连接）来描述图形的边界。这在内存方面是有效的，并且可以使用 MatplotlibPatches轻松在屏幕上可视化它们。但是，这种方法使确定点或其他图形是否位于给定图形内变得更加困难。这是许多几何问题中的一个关键问题。
在本教程中，我们将学习如何表示几何图形并确定点是否位于图形内。
做好准备
对于本教程，我们需要将matplotlib包（整体）导入为mpl，将pyplot模块导入为plt：
import matplotlib as mpl
import matplotlib.pyplot as plt
我们还需要从 Shapely 包的geometry模块中导入Point和Polygon对象。Shapely 包包含许多用于表示、操作和分析二维几何图形的例程和对象：
from shapely.geometry import Polygon, Point
如何操作…
以下步骤向您展示如何创建多边形的 Shapely 表示，然后测试点是否位于此多边形内：
- 创建一个样本多边形进行测试：
polygon = Polygon(
[(0, 2), (-1, 1), (-0.5, -1), (0.5, -1), (1, 1)],
)
- 接下来，我们在新图上绘制多边形。首先，我们需要将多边形转换为可以添加到图中的 Matplotlib
Polygon补丁：
fig, ax = plt.subplots()
poly_patch = mpl.patches.Polygon(polygon.exterior, ec="k",
lw="1", alpha=0.5)
ax.add_patch(poly_patch)
ax.set(xlim=(-1.05, 1.05), ylim=(-1.05, 2.05))
ax.set_axis_off()
- 现在，我们需要创建两个测试点，其中一个将位于多边形内，另一个将位于多边形外：
p1 = Point(0.0, 0.0)
p2 = Point(-1.0, -0.75)
- 我们在多边形上方绘制并注释这两个点，以显示它们的位置：
ax.plot(0.0, 0.0, "k*")
ax.annotate("p1", (0.0, 0.0), (0.05, 0.0))
ax.plot(-0.8, -0.75, "k*")
ax.annotate("p2", (-0.8, -0.75), (-0.8 + 0.05, -0.75))
- 最后，我们使用
contains方法测试每个点是否位于多边形内，然后将结果打印到终端：
print("p1 inside polygon?", polygon.contains(p1))
print("p2 inside polygon?", polygon.contains(p2))
结果显示，第一个点p1包含在多边形内，而第二个点p2不包含。这也可以在下图中看到，清楚地显示了一个点包含在阴影多边形内，而另一个点不包含：
图 8.2：多边形区域内外的点
工作原理…
ShapelyPolygon类是多边形的表示，它将其顶点存储为点。外边界包围的区域-存储顶点之间的五条直线-对我们来说是明显的，并且很容易被眼睛识别，但是“在”边界内的概念很难以一种计算机容易理解的方式定义。甚至很难给出关于“在”给定曲线内的含义的正式数学定义。
确定点是否位于简单闭合曲线内有两种主要方法 - 即从同一位置开始并结束且不包含任何自交点的曲线。第一种方法使用数学概念称为绕数，它计算曲线“绕”点的次数，以及射线交叉计数方法，其中我们计算从点到无穷远处的点的射线穿过曲线的次数。幸运的是，我们不需要自己计算这些数字，因为我们可以使用 Shapely 包中的工具来为我们执行这些计算。这就是多边形的contains方法所做的。（在底层，Shapely 使用 GEOS 库执行此计算。）
Shapely 的Polygon类可用于计算与这些平面图形相关的许多数量，包括周长和面积。contains方法用于确定点或一组点是否位于对象表示的多边形内（该类有关多边形的表示存在一些限制）。实际上，您可以使用相同的方法来确定一个多边形是否包含在另一个多边形内，因为正如我们在这个示例中看到的，多边形由一组简单的点表示。
在图像中找到边缘
在图像中找到边缘是减少包含大量噪音和干扰的复杂图像为包含最突出轮廓的非常简单图像的好方法。这可以作为我们分析过程的第一步，例如在图像分类中，或者作为将线轮廓导入计算机图形软件包的过程。
在这个示例中，我们将学习如何使用scikit-image包和 Canny 算法来找到复杂图像中的边缘。
准备工作
对于这个示例，我们需要将 Matplotlib 的pyplot模块导入为plt，从skimage.io模块导入imread例程，以及从skimage.feature模块导入canny例程：
import matplotlib.pyplot as plt
from skimage.io import imread
from skimage.feature import canny
如何做到…
按照以下步骤学习如何使用scikit-image包在图像中找到边缘：
- 从源文件加载图像数据。这可以在本章的 GitHub 存储库中找到。关键是，我们传入
as_gray=True以以灰度加载图像：
image = imread("mandelbrot.png", as_gray=True)
以下是原始图像，供参考。集合本身由白色区域显示，如您所见，边界由较暗的阴影表示，非常复杂：
图 8.3：使用 Python 生成的 Mandelbrot 集合的绘图
- 接下来，我们使用
canny例程，需要从scikit-image包的features模块导入。对于这个图像，sigma值设置为 0.5：
edges = canny(image, sigma=0.5)
- 最后，我们将
edges图像添加到一个新的图中，使用灰度（反转）色图：
fig, ax = plt.subplots()
ax.imshow(edges, cmap="gray_r")
ax.set_axis_off()
已检测到的边缘可以在以下图像中看到。边缘查找算法已经识别出 Mandelbrot 集合边界的大部分可见细节，尽管并不完美（毕竟这只是一个估计）：
图 8.4：使用 scikit-image 包的 Canny 边缘检测算法找到的 Mandelbrot 集合的边缘
它是如何工作的…
scikit-image包提供了各种用于操作和分析从图像中导出的数据的实用程序和类型。正如其名称所示，canny例程使用 Canny 边缘检测算法来找到图像中的边缘。该算法使用图像中的强度梯度来检测边缘，其中梯度较大。它还执行一些过滤以减少它找到的边缘中的噪音。
我们提供的sigma关键字值是应用于图像的高斯平滑的标准偏差，用于计算边缘检测。这有助于我们去除图像中的一些噪音。我们设置的值（0.5）小于默认值（1），但在这种情况下可以提供更好的分辨率。较大的值会遮盖 Mandelbrot 集边界的一些细节。
对平面图形进行三角剖分
正如我们在第三章中看到的，微积分和微分方程，我们经常需要将连续区域分解为更小、更简单的区域。在之前的示例中，我们将实数区间缩小为一系列长度较小的小区间。这个过程通常称为离散化。在本章中，我们正在处理二维图形，因此我们需要这个过程的二维版本。为此，我们将一个二维图形（在这个示例中是一个多边形）分解为一系列更小和更简单的多边形。所有多边形中最简单的是三角形，因此这是二维离散化的一个很好的起点。找到一组"铺砌"几何图形的三角形的过程称为三角剖分。
在这个示例中，我们将学习如何使用 Shapely 包对多边形（带有孔）进行三角剖分。
准备工作
在这个示例中，我们需要将 NumPy 包导入为np，将 Matplotlib 包导入为mpl，并将pyplot模块导入为plt。
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
我们还需要从 Shapely 包中获取以下项目：
from shapely.geometry import Polygon
from shapely.ops import triangulate
如何做…
以下步骤向您展示了如何使用 Shapely 包对带有孔的多边形进行三角剖分：
- 首先，我们需要创建一个代表我们希望进行三角剖分的图形的
Polygon对象：
polygon = Polygon(
[(2.0, 1.0), (2.0, 1.5), (-4.0, 1.5), (-4.0, 0.5),
(-3.0, -1.5), (0.0, -1.5), (1.0, -2.0), (1.0, -0.5),
(0.0, -1.0), (-0.5, -1.0), (-0.5, 1.0)],
holes=[np.array([[-1.5, -0.5], [-1.5, 0.5], [-2.5, 0.5],
[-2.5, -0.5]])]
)
- 现在，我们应该绘制图形，以便了解我们将在其中工作的区域：
fig, ax = plt.subplots()
plt_poly = mpl.patches.Polygon(polygon.exterior,
ec="k", lw="1", alpha=0.5, zorder=0)
ax.add_patch(plt_poly)
plt_hole = mpl.patches.Polygon(polygon.interiors[0],
ec="k", fc="w")
ax.add_patch(plt_hole)
ax.set(xlim=(-4.05, 2.05), ylim=(-2.05, 1.55))
ax.set_axis_off()
可以在下图中看到这个多边形。正如我们所看到的，这个图形中有一个"孔"，必须仔细考虑：
图 8.5：带有孔的示例多边形
- 我们使用
triangulate例程生成多边形的三角剖分。这个三角剖分包括外部边缘，这是我们在这个示例中不想要的：
triangles = triangulate(polygon)
- 为了去除位于原始多边形外部的三角形，我们需要使用内置的
filter例程，以及contains方法（在本章前面已经看到）：
filtered = filter(lambda p: polygon.contains(p), triangles)
- 将三角形绘制在原始多边形上，我们需要将 Shapely 三角形转换为 Matplotlib
Patch对象，然后将其存储在PatchCollection中：
patches = map(lambda p: mpl.patches.Polygon(p.exterior), filtered)
col = mpl.collections.PatchCollection(patches, fc="none", ec="k")
- 最后，我们将三角形补丁的集合添加到之前创建的图形中：
ax.add_collection(col)
在原始多边形上绘制的三角剖分可以在下图中看到。在这里，我们可以看到每个顶点都连接到另外两个顶点，形成了覆盖整个原始多边形的三角形系统：
图 8.6：带有孔的示例多边形的三角剖分
它是如何工作的…
triangulate例程使用一种称为Delaunay 三角剖分的技术将一组点连接到一组三角形中。在这种情况下，这组点是多边形的顶点。Delaunay 方法以这样一种方式找到这些三角形，即没有任何点包含在任何三角形的外接圆内。这是该方法的技术条件，但这意味着三角形被有效地选择，因为它避免了非常长、细的三角形。得到的三角剖分利用了原始多边形中存在的边缘，并连接了一些外部边缘。
为了去除原多边形外的三角形，我们使用内置的filter例程，它通过移除标准函数失败的项目来创建一个新的可迭代对象。这与 Shapely Polygon对象上的contains方法一起使用，以确定每个三角形是否位于原始图形内。正如我们之前提到的，我们需要将这些 Shapely 项目转换为 Matplotlib 补丁，然后才能将它们添加到图中。
还有更多…
三角剖分通常用于将复杂的几何图形简化为一组三角形，这些三角形对于某种计算任务来说要简单得多。然而，它们也有其他用途。三角剖分的一个特别有趣的应用是解决“艺术画廊问题”。这个问题涉及找到必要的“守卫”艺术画廊的最大数量。三角剖分是 Fisk 对艺术画廊定理的简单证明的重要部分，这个定理最初是由 Chvátal 证明的。
假设这个食谱中的多边形是一个艺术画廊的平面图，并且一些守卫需要放置在顶点上。一点工作就会表明，你需要在多边形的顶点处放置三个守卫，整个博物馆才能被覆盖。在下面的图像中，我们绘制了一个可能的布局：
图 8.7：在顶点上放置守卫的艺术画廊问题的一个可能解决方案。
点由点表示，并且它们相应的视野范围被阴影表示。
每个顶点都放置了一个守卫，并且他们的视野范围由相应的阴影区域表示。在这里，你可以看到整个多边形至少被一种颜色覆盖。艺术画廊问题的解决方案——实际上是原问题的一个变体——告诉我们，最多需要四名守卫。
另请参阅
关于艺术画廊问题的更多信息可以在 O’Rourke 的经典著作中找到：ORourke, J. (1987). Art gallery theorems and algorithms. New York: Oxford University Press.
计算凸包
如果图形内的每一对点都可以使用一条直线连接，并且这条直线也包含在图形内，那么几何图形被称为凸。凸体的简单例子包括点、直线、正方形、圆（圆盘）、正多边形等。图 8.5 中显示的几何图形不是凸的，因为孔的对面的点不能通过保持在图形内的直线连接起来。
从某种角度来看，凸图形是简单的，这意味着它们在各种应用中都很有用。一个特别的问题涉及找到包含一组点的最小凸集。这个最小凸集被称为这组点的凸包。
在这个食谱中，我们将学习如何使用 Shapely 包找到一组点的凸包。
准备工作
对于这个食谱，我们需要导入 NumPy 包作为np，导入 Matplotlib 包作为mpl，并导入plt模块：
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
我们还需要从 NumPy 导入默认的随机数生成器。我们可以这样导入：
from numpy.random import default_rng
rng = default_rng(12345)
最后，我们需要从 Shapely 导入MultiPoint类：
from shapely.geometry import MultiPoint
操作方法…
按照以下步骤找到一组随机生成点的凸包：
- 首先，我们生成一个二维数组的随机数：
raw_points = rng.uniform(-1.0, 1.0, size=(50, 2))
- 接下来，我们创建一个新图形，并在这个图形上绘制这些原始样本点：
fig, ax = plt.subplots()
ax.plot(raw_points[:, 0], raw_points[:, 1], "k.")
ax.set_axis_off()
这些随机生成的点可以在下图中看到。这些点大致分布在一个正方形区域内：
图 8.8：平面上的一组点
- 接下来，我们构建一个
MultiPoint对象，收集所有这些点并将它们放入一个单一对象中：
points = MultiPoint(raw_points)
- 现在，我们使用
convex_hull属性获取这个MultiPoint对象的凸包：
convex_hull = points.convex_hull
- 然后，我们创建一个 Matplotlib
Polygon补丁，可以在我们的图中绘制，以显示找到的凸包的结果：
patch = mpl.patches.Polygon(convex_hull.exterior, alpha=0.5,
ec="k", lw=1.2)
- 最后，我们将
Polygon补丁添加到图中，以显示凸包：
ax.add_patch(patch)
随机生成的点的凸包可以在下图中看到：
图 8.9：平面上一组点的凸包
工作原理…
Shapely 包是围绕 GEOS 库的 Python 包装器，用于几何分析。Shapely 几何对象的convex_hull属性调用 GEOS 库中的凸包计算例程，从而产生一个新的 Shapely 对象。从这个教程中，我们可以看到一组点的凸包是一个多边形，其顶点是离“中心”最远的点。
构造贝塞尔曲线
贝塞尔曲线，或B 样条，是一族曲线，在矢量图形中非常有用-例如，它们通常用于高质量的字体包中。这是因为它们由少量点定义，然后可以用来廉价地计算沿曲线的大量点。这允许根据用户的需求来缩放细节。
在本教程中，我们将学习如何创建一个表示贝塞尔曲线的简单类，并计算沿其路径的若干点。
准备工作
在本教程中，我们将使用导入为np的 NumPy 包，导入为plt的 Matplotlib pyplot模块，以及 Python 标准库math模块中导入为binom的comb例程：
from math import comb as binom
import matplotlib.pyplot as plt
import numpy as np
如何做…
按照以下步骤定义一个表示贝塞尔曲线的类，该类可用于计算沿曲线的点：
- 第一步是设置基本类。我们需要为实例属性提供控制点（节点）和一些相关的数字：
class Bezier:
def __init__(self, *points):
self.points = points
self.nodes = n = len(points) - 1
self.degree = l = points[0].size
- 仍然在
__init__方法中，我们生成贝塞尔曲线的系数，并将它们存储在实例属性的列表中：
self.coeffs = [binom(n, i)*p.reshape((l, 1)) for i,
p in enumerate(points)]
- 接下来，我们定义一个
__call__方法，使类可调用。我们将实例中的节点数加载到本地变量中，以便清晰明了：
def __call__(self, t):
n = self.nodes
- 接下来，我们重新整理输入数组，使其包含单行：
t = t.reshape((1, t.size))
- 现在，我们使用实例的
coeffs属性中的每个系数生成值数组的列表：
vals = [c @ (t**i)*(1-t)**(n-i) for i,
c in enumerate(self.coeffs)]
- 最后，我们对步骤 5中构造的所有数组进行求和，并返回结果数组：
return np.sum(vals, axis=0)
- 现在，我们将通过一个示例来测试我们的类。我们将为此示例定义四个控制点：
p1 = np.array([0.0, 0.0])
p2 = np.array([0.0, 1.0])
p3 = np.array([1.0, 1.0])
p4 = np.array([1.0, 3.0])
- 接下来，我们为绘图设置一个新的图形，并用虚线连接线绘制控制点：
fig, ax = plt.subplots()
ax.plot([0.0, 0.0, 1.0, 1.0], [0.0, 1.0, 1.0, 3.0], "*--k")
ax.set(xlabel="x", ylabel="y", title="Bezier curve with
4 nodes, degree 3")
- 然后，我们使用步骤 7中定义的四个点创建我们的
Bezier类的新实例：
b_curve = Bezier(p1, p2, p3, p4)
- 现在，我们可以使用
linspace创建 0 到 1 之间等间距点的数组，并计算沿着贝塞尔曲线的点：
t = np.linspace(0, 1)
v = b_curve(t)
- 最后，我们在之前绘制的控制点上绘制这条曲线：
ax.plot(v[0,:], v[1, :])
我们绘制的贝塞尔曲线可以在下图中看到。正如你所看到的，曲线从第一个点（0, 0）开始，结束于最终点（1, 3）：
图 8.10：使用四个节点构造的三次贝塞尔曲线
工作原理…
贝塞尔曲线由一系列控制点描述，我们以递归方式构造曲线。一个点的贝塞尔曲线是一个保持在该点的常数曲线。具有两个控制点的贝塞尔曲线是这两个点之间的线段：
当我们添加第三个控制点时，我们取对应点之间的线段，这些点是由一个较少点构成的贝塞尔曲线的曲线。这意味着我们使用以下公式构造具有三个控制点的贝塞尔曲线：
这种构造可以在下图中看到：
图 8.11：使用递归定义构造二次贝塞尔曲线。黑色虚线显示了两条线性贝塞尔曲线。
这种构造方式继续定义了任意数量控制点上的贝塞尔曲线。幸运的是，在实践中我们不需要使用这种递归定义，因为我们可以将公式展开成曲线的单一公式，即以下公式：
在这里，p[i]元素是控制点，t是一个参数，而
是二项式系数。请记住，t参数是生成曲线点的变化量。我们可以分离前述求和中涉及t的项和不涉及t的项。这定义了我们在步骤 2中定义的系数，每个系数由以下代码片段给出：
binom(n, i)*p.reshape((l, 1))
我们在这一步中对每个点p进行了 reshape，以确保它被排列为列向量。这意味着每个系数都是一个列向量（作为 NumPy 数组），由二项式系数缩放的控制点组成。
现在，我们需要指定如何在不同的t值上评估贝塞尔曲线。这就是我们利用 NumPy 包中的高性能数组操作的地方。在形成系数时，我们将控制点 reshape 为列向量。在步骤 4中，我们将输入t值 reshape 为行向量。这意味着我们可以使用矩阵乘法运算符将每个系数乘以相应的（标量）值，具体取决于输入的t。这就是步骤 5中列表推导式中发生的情况。在下一行中，我们将l×1数组乘以1×N数组，得到一个l×N数组：
c @ (t**i)*(1-t)**(n-i)
我们为每个系数都得到一个这样的数组。然后，我们可以使用np.sum例程来对这些l×N数组中的每一个进行求和，以得到贝塞尔曲线上的值。在本示例中，输出数组的顶行包含曲线的x值，底行包含曲线的y值。在指定axis=0关键字参数时，我们必须小心确保sum例程对我们创建的列表进行求和，而不是对该列表包含的数组进行求和。
我们定义的类是使用贝塞尔曲线的控制点进行初始化的，然后用于生成系数。曲线值的实际计算是使用 NumPy 完成的，因此这种实现应该具有相对良好的性能。一旦创建了这个类的特定实例，它的功能就非常像一个函数，正如你所期望的那样。但是，这里没有进行类型检查，所以我们只能用 NumPy 数组作为参数来调用这个“函数”。
还有更多…
贝塞尔曲线是使用迭代构造定义的，其中具有n个点的曲线是使用连接由第一个和最后一个n-1点定义的曲线来定义的。使用这种构造跟踪每个控制点的系数将很快导致我们用来定义前述曲线的方程。这种构造还导致贝塞尔曲线的有趣和有用的几何特性。
正如我们在这个配方的介绍中提到的，贝塞尔曲线出现在许多涉及矢量图形的应用程序中，比如字体。它们也出现在许多常见的矢量图形软件包中。在这些软件包中，通常会看到二次贝塞尔曲线，它们由三个点的集合定义。然而，你也可以通过提供两个端点以及这些点上的梯度线来定义一个二次贝塞尔曲线。这在图形软件包中更常见。生成的贝塞尔曲线将沿着梯度线离开每个端点，并在这些点之间平滑地连接曲线。
我们在这里构建的实现对于小型应用程序来说性能相对较好，但对于涉及在大量t值上渲染具有大量控制点的曲线的应用程序来说是不够的。对于这一点，最好使用一个用编译语言编写的低级软件包。例如，bezier Python 软件包使用编译的 Fortran 后端进行计算，并提供比我们在这里定义的类更丰富的接口。
当然，贝塞尔曲线可以自然地扩展到更高的维度。结果是一个贝塞尔曲面，使它们成为非常有用的通用工具，用于高质量、可伸缩的图形。
进一步阅读
计算几何中一些常见算法的描述可以在以下书籍中找到：Press, W.H., Teukolsky, S.A., Vetterling, W.T., and Flannery, B.P., 2007. Numerical recipes: the art of scientific computing*. 3rd ed. Cambridge: Cambridge University Press*。
有关计算几何中一些问题和技术的更详细描述，请查阅以下书籍：O’Rourke, J., 1994. Computational geometry in C*. Cambridge: Cambridge University Press*。
第十章：寻找最优解
在本章中，我们将讨论寻找给定情况下最佳结果的各种方法。这被称为优化，通常涉及最小化或最大化目标函数。目标函数是一个接受多个参数作为参数并返回代表给定参数选择的成本或回报的单个标量值的函数。关于最小化和最大化函数的问题实际上是相互等价的，因此我们只会在本章讨论最小化目标函数。最小化函数f(x)等同于最大化函数*-f*(x)。在我们讨论第一个配方时将提供更多细节。
我们可以利用的算法来最小化给定函数取决于函数的性质。例如，包含一个或多个变量的简单线性函数与具有许多变量的非线性函数相比，可用的算法不同。线性函数的最小化属于线性规划范畴，这是一个发展完善的理论。对于非线性函数，我们通常利用函数的梯度（导数）来寻找最小点。我们将讨论几种不同类型函数的最小化方法。
寻找单变量函数的极小值和极大值特别简单，如果函数的导数已知，可以轻松完成。如果不知道导数，则适用于适当配方的方法。最小化非线性函数配方中的注释提供了一些额外细节。
我们还将提供一个非常简短的介绍博弈论。广义上讲，这是一个围绕决策制定的理论，并在经济学等学科中具有广泛的影响。特别是，我们将讨论如何在 Python 中将简单的双人游戏表示为对象，计算与某些选择相关的回报，并计算这些游戏的纳什均衡。
我们将首先看如何最小化包含一个或多个变量的线性和非线性函数。然后，我们将继续研究梯度下降方法和使用最小二乘法进行曲线拟合。最后，我们将通过分析双人游戏和纳什均衡来结束本章。
在本章中，我们将涵盖以下配方：
最小化简单线性函数
最小化非线性函数
使用梯度下降法进行优化
使用最小二乘法拟合数据的曲线
分析简单的双人游戏
计算纳什均衡
让我们开始吧！
技术要求
在本章中，我们将像往常一样需要 NumPy 包、SciPy 包和 Matplotlib 包。我们还将需要 Nashpy 包用于最后两个配方。这些包可以使用您喜欢的包管理器（如pip）进行安装：
python3.8 -m pip install numpy scipy matplotlib nashpy
本章的代码可以在 GitHub 存储库的Chapter 09文件夹中找到，网址为github.com/PacktPublishing/Applying-Math-with-Python/tree/master/Chapter%2009。
查看以下视频以查看代码的实际操作：bit.ly/2BjzwGo。
最小化简单线性函数
在优化中我们面临的最基本问题是找到函数取得最小值的参数。通常，这个问题受到参数可能值的一些限制的约束，这增加了问题的复杂性。显然，如果我们要最小化的函数也很复杂，那么这个问题的复杂性会进一步增加。因此，我们必须首先考虑线性函数，它们的形式如下：
为了解决这些问题，我们需要将约束转化为计算机可用的形式。在这种情况下，我们通常将它们转化为线性代数问题（矩阵和向量）。一旦完成了这一步，我们就可以使用 NumPy 和 SciPy 中的线性代数包中的工具来找到我们所寻求的参数。幸运的是，由于这类问题经常发生，SciPy 有处理这种转化和随后求解的例程。
在这个配方中，我们将使用 SciPy optimize模块的例程来解决以下受限线性最小化问题：
这将受到以下条件的约束：
准备工作
对于这个配方，我们需要在别名np下导入 NumPy 包，以plt的名称导入 Matplotlib pyplot模块，以及 SciPy optimize模块。我们还需要从mpl_toolkits.mplot3d导入Axes3D类，以使 3D 绘图可用：
import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
如何做…
按照以下步骤使用 SciPy 解决受限线性最小化问题：
- 将系统设置为 SciPy 可以识别的形式：
A = np.array([
[2, 1], # 2*x0 + x1 <= 6
[-1, -1] # -x0 - x1 <= -4
])
b = np.array([6, -4])
x0_bounds = (-3, 14) # -3 <= x0 <= 14
x1_bounds = (2, 12) # 2 <= x1 <= 12
c = np.array([1, 5])
- 接下来，我们需要定义一个评估线性函数在x值处的例程，这是一个向量（NumPy 数组）：
def func(x):
return np.tensordot(c, x, axes=1)
- 然后，我们创建一个新的图并添加一组
3d轴，我们可以在上面绘制函数：
fig = plt.figure()
ax = fig.add_subplot(projection="3d")
ax.set(xlabel="x0", ylabel="x1", zlabel="func")
ax.set_title("Values in Feasible region")
- 接下来，我们创建一个覆盖问题区域的值网格，并在该区域上绘制函数的值：
X0 = np.linspace(*x0_bounds)
X1 = np.linspace(*x1_bounds)
x0, x1 = np.meshgrid(X0, X1)
z = func([x0, x1])
ax.plot_surface(x0, x1, z, alpha=0.3)
- 现在，我们在函数值平面上绘制与临界线
2*x0 + x1 == 6对应的线，并在我们的图上绘制落在范围内的值：
Y = (b[0] - A[0, 0]*X0) / A[0, 1]
I = np.logical_and(Y >= x1_bounds[0], Y <= x1_bounds[1])
ax.plot(X0[I], Y[I], func([X0[I], Y[I]]), "r", lw=1.5)
- 我们重复这个绘图步骤，对第二条临界线
x0 + x1 == -4：
Y = (b[1] - A[1, 0]*X0) / A[1, 1]
I = np.logical_and(Y >= x1_bounds[0], Y <= x1_bounds[1])
ax.plot(X0[I], Y[I], func([X0[I], Y[I]]), "r", lw=1.5)
- 接下来，我们着色位于两条临界线之间的区域，这对应于最小化问题的可行区域：
B = np.tensordot(A, np.array([x0, x1]), axes=1)
II = np.logical_and(B[0, ...] <= b[0], B[1, ...] <= b[1])
ax.plot_trisurf(x0[II], x1[II], z[II], color="b", alpha=0.5)
函数值在可行区域上的图可以在以下图片中看到：
图 9.1：突出显示了可行区域的线性函数值
正如我们所看到的，位于这个着色区域内的最小值发生在两条临界线的交点处。
- 接下来，我们使用
linprog来解决带有我们在步骤 1中创建的边界的受限最小化问题。我们在终端中打印出结果对象：
res = optimize.linprog(c, A_ub=A, b_ub=b, bounds=
(x0_bounds, x1_bounds))
print(res)
- 最后，我们在可行区域上绘制最小函数值：
ax.plot([res.x[0]], [res.x[1]], [res.fun], "k*")
更新后的图可以在以下图片中看到：
图 9.2：在可行区域上绘制的最小值
在这里，我们可以看到linprog例程确实发现了最小值在两条临界线的交点处。
工作原理…
受限线性最小化问题在经济情况中很常见，您尝试在保持参数的其他方面的同时最小化成本。实际上，优化理论中的许多术语都反映了这一事实。解决这些问题的一个非常简单的算法称为单纯形法，它使用一系列数组操作来找到最小解。从几何上讲，这些操作代表着改变单纯形的不同顶点（我们在这里不定义），正是这一点赋予了算法其名称。
在我们继续之前，我们将简要概述单纯形法用于解决受限线性优化问题的过程。我们所面临的问题不是一个矩阵方程问题，而是一个矩阵不等式问题。我们可以通过引入松弛变量来解决这个问题，将不等式转化为等式。例如，通过引入松弛变量s[1]，可以将第一个约束不等式重写如下：
只要*s[1]不是负数，这就满足了所需的不等式。第二个约束不等式是大于或等于类型的不等式，我们必须首先将其更改为小于或等于类型。我们通过将所有项乘以-1 来实现这一点。这给了我们在配方中定义的矩阵A的第二行。引入第二个松弛变量s[2]*后，我们得到了第二个方程：
从中，我们可以构建一个矩阵，其列包含两个参数变量*x[1]和x[2]以及两个松弛变量s[1]和s[2]的系数。该矩阵的行代表两个边界方程和目标函数。现在可以使用对该矩阵进行初等行操作来解决这个方程组，以获得最小化目标函数的x[1]和x[2]*的值。由于解决矩阵方程很容易且快速，这意味着我们可以快速高效地最小化线性函数。
幸运的是，我们不需要记住如何将不等式系统化简为线性方程组，因为诸如linprog之类的例程会为我们完成这一工作。我们只需将边界不等式提供为矩阵和向量对，包括每个系数，以及定义目标函数的单独向量。linprog例程负责制定和解决最小化问题。
实际上，单纯形法并不是linprog例程用于最小化函数的算法。相反，linprog使用内点算法，这更有效率。（可以通过提供method关键字参数并使用适当的方法名称将方法设置为simplex或revised-simplex。在打印的输出结果中，我们可以看到只用了五次迭代就达到了解决方案。）该例程返回的结果对象包含最小值发生的参数值存储在x属性中，该最小值存储在fun属性中，以及有关解决过程的各种其他信息。如果方法失败，那么status属性将包含一个数值代码，描述方法失败的原因。
在本文的步骤 2中，我们创建了一个代表此问题的目标函数的函数。该函数以单个数组作为输入，该数组包含应在其上评估函数的参数空间值。在这里，我们使用了 NumPy 的tensordot例程（带有axes=1）来评估系数向量c与每个输入x的点积。在这里我们必须非常小心，因为我们传递给函数的值在后续步骤中将是一个 2×50×50 的数组。普通的矩阵乘法（np.dot）在这种情况下不会给出我们所期望的 50×50 数组输出。
在步骤 5和6中，我们计算了临界线上的点，这些点具有以下方程：
然后我们计算了相应的z值，以便绘制在由目标函数定义的平面上的线。我们还需要“修剪”这些值，以便只包括在问题中指定范围内的值。
还有更多…
本文介绍了受约束的最小化问题以及如何使用 SciPy 解决它。然而，相同的方法也可以用于解决受约束的最大化问题。这是因为最大化和最小化在某种意义上是对偶的，即最大化函数f(x)等同于最小化函数*-f*(x)，然后取其负值。事实上，我们在本文中使用了这一事实，将第二个约束不等式从≥改为≤。
在这个示例中，我们解决了一个只有两个参数变量的问题，但是相同的方法将适用于涉及两个以上这样的变量的问题（除了绘图步骤）。我们只需要向每个数组添加更多的行和列，以考虑这增加的变量数量 - 这包括提供给例程的边界元组。在处理非常大量的变量时，例程也可以与稀疏矩阵一起使用，以获得额外的效率。
linprog例程的名称来自线性规划，用于描述这种类型的问题 - 找到满足一些矩阵不等式的x的值，受其他条件的限制。由于与矩阵理论和线性代数的理论有非常紧密的联系，因此在线性规划问题中有许多非常快速和高效的技术，这些技术在非线性情境中是不可用的。
最小化非线性函数
在上一个示例中，我们看到了如何最小化一个非常简单的线性函数。不幸的是，大多数函数都不是线性的，通常也没有我们希望的良好性质。对于这些非线性函数，我们不能使用为线性问题开发的快速算法，因此我们需要设计可以在这些更一般情况下使用的新方法。我们将使用的算法称为 Nelder-Mead 算法，这是一种健壮且通用的方法，用于找到函数的最小值，并且不依赖于函数的梯度。
在这个示例中，我们将学习如何使用 Nelder-Mead 单纯形法来最小化包含两个变量的非线性函数。
准备工作
在这个示例中，我们将使用导入为np的 NumPy 包，导入为plt的 Matplotlib pyplot模块，从mpl_toolkits.mplot3d导入的Axes3D类以启用 3D 绘图，以及 SciPy optimize模块：
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy import optimize
如何做…
以下步骤向您展示如何使用 Nelder-Mead 单纯形法找到一般非线性目标函数的最小值：
- 定义我们将最小化的目标函数：
def func(x):
return ((x[0] - 0.5)**2 + (x[1] + 0.5)**2)*
np.cos(0.5*x[0]*x[1])
- 接下来，创建一个值网格，我们可以在上面绘制我们的目标函数：
x_r = np.linspace(-1, 1)
y_r = np.linspace(-2, 2)
x, y = np.meshgrid(x_r, y_r)
- 现在，我们在这个点网格上评估函数：
z = func([x, y])
- 接下来，我们创建一个带有
3d轴对象的新图，并设置轴标签和标题：
fig = plt.figure(tight_layout=True)
ax = fig.add_subplot(projection="3d")
ax.tick_params(axis="both", which="major", labelsize=9)
ax.set(xlabel="x", ylabel="y", zlabel="z")
ax.set_title("Objective function")
- 现在，我们可以在我们刚刚创建的轴上将目标函数绘制为表面：
ax.plot_surface(x, y, z, alpha=0.7)
- 我们选择一个初始点，我们的最小化例程将从该点开始迭代，并在表面上绘制这个点：
x0 = np.array([-0.5, 1.0])
ax.plot([x0[0]], [x0[1]], func(x0), "r*")
可以在以下图像中看到目标函数表面的绘图，以及初始点。在这里，我们可以看到最小值似乎出现在 x 轴上约 0.5，y 轴上约-0.5 的位置：
图 9.3：具有起始值的非线性目标函数
- 现在，我们使用
optimize包中的minimize例程来找到最小值，并打印它产生的result对象：
result = optimize.minimize(func, x0, tol=1e-6, method=
"Nelder-Mead")
print(result)
- 最后，我们在目标函数表面上绘制
minimize例程找到的最小值：
ax.plot([result.x[0]], [result.x[1]], [result.fun], "r*")
包括minimize例程找到的最小点的目标函数的更新绘图可以在以下图像中看到：
图 9.4：具有起始点和最小点的目标函数
它是如何工作的…
Nelder-Mead 单纯形法-不要与线性优化问题的单纯形法混淆-是一种简单的算法，用于找到非线性函数的最小值，即使目标函数没有已知的导数也可以工作。（这不适用于此示例中的函数；使用基于梯度的方法的唯一收益是收敛速度。）该方法通过比较单纯形的顶点处的目标函数值来工作，在二维空间中是一个三角形。具有最大函数值的顶点通过相反的边被“反射”，并执行适当的扩展或收缩，实际上将单纯形“下坡”移动。
SciPyoptimize模块中的minimize例程是许多非线性函数最小化算法的入口点。在这个示例中，我们使用了 Nelder-Mead 单纯形算法，但还有许多其他可用的算法。其中许多算法需要对函数的梯度有所了解，该梯度可能会被算法自动计算。可以通过向method关键字参数提供适当的名称来使用该算法。
“最小化”例程返回的result对象包含有关求解器找到的解决方案的大量信息-或者如果发生错误，则未找到-。特别是，计算出的最小值发生的期望参数存储在结果的x属性中，而函数值存储在fun属性中。
“最小化”例程需要函数和x0的起始值。在这个示例中，我们还提供了一个容差值，最小值应该使用tol关键字参数计算。更改此值将修改计算解的准确度。
还有更多…
Nelder-Mead 算法是“无梯度”最小化算法的一个例子，因为它不需要任何关于目标函数的梯度（导数）的信息。有几种这样的算法，通常涉及在指定点评估目标函数，然后使用这些信息朝向最小值移动。一般来说，无梯度方法的收敛速度比梯度下降模型慢。但是，它们可以用于几乎任何目标函数，即使很难精确计算梯度或通过近似手段计算。
通常，优化单变量函数比多维情况更容易，并且在 SciPyoptimize库中有其专用函数。minimize_scalar例程对单变量函数执行最小化，并且在这种情况下应该使用而不是minimize。
在优化中使用梯度下降方法
在上一个示例中，我们使用 Nelder-Mead 单纯形算法最小化包含两个变量的非线性函数。这是一种相当健壮的方法，即使对目标函数了解甚少也可以工作。然而，在许多情况下，我们对目标函数了解更多，这一事实使我们能够设计更快和更有效的最小化函数的算法。我们可以通过利用函数的梯度等属性来做到这一点。
多于一个变量的函数的梯度描述了函数在各个分量方向上的变化率。这是函数对每个变量的偏导数的向量。从这个梯度向量中，我们可以推断出函数在哪个方向上增长最快，反之亦然，从任意给定位置开始，函数在哪个方向上下降最快。这为梯度下降方法最小化函数提供了基础。算法非常简单：给定一个起始位置x，我们计算在这个x处的梯度以及梯度最快下降的相应方向，然后沿着那个方向迈出一小步。经过几次迭代，这将从起始位置移动到函数的最小值。
在这个示例中，我们将学习如何实现基于最陡下降算法的算法，以在有界区域内最小化目标函数。
准备工作
对于这个示例，我们需要导入 NumPy 包作为np，导入 Matplotlib 的pyplot模块作为plt，并从mpl_toolkits.mplot3d导入Axes3D对象：
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
如何做…
在接下来的步骤中，我们将实现一个简单的梯度下降方法，以最小化具有已知梯度函数的目标函数（实际上我们将使用一个生成器函数，以便我们可以看到方法的工作方式）：
- 我们将首先定义一个
descend例程，它将执行我们的算法。函数声明如下：
def descend(func, x0, grad, bounds, tol=1e-8, max_iter=100):
- 接下来，我们需要实现这个例程。我们首先定义将在方法运行时保存迭代值的变量：
xn = x0
xnm1 = np.inf
grad_xn = grad(x0)
- 然后，我们开始我们的循环，这将运行迭代。我们立即检查是否在继续之前正在取得有意义的进展：
for i in range(max_iter):
if np.linalg.norm(xn - xnm1) < tol:
break
- 方向是梯度向量的负。我们计算一次并将其存储在
direction变量中：
direction = -grad_xn
- 现在，我们更新先前和当前的值，分别为
xnm1和xn，准备进行下一次迭代。这结束了descend例程的代码：
xnm1 = xn
xn = xn + 0.2*direction
- 现在，我们可以计算当前值的梯度并产生所有适当的值：
grad_xn = grad(xn)
yield i, xn, func(xn), grad_xn
这结束了descend例程的定义。
- 我们现在可以定义一个要最小化的样本目标函数：
def func(x):
return ((x[0] - 0.5)**2 + (x[1] + 0.5)**2)*np.cos(0.5*x[0]*x[1])
- 接下来，我们创建一个网格，我们将评估并绘制目标函数：
x_r = np.linspace(-1, 1)
y_r = np.linspace(-2, 2)
x, y = np.meshgrid(x_r, y_r)
- 一旦网格创建完成，我们可以评估我们的函数并将结果存储在
z变量中：
z = func([x, y])
- 接下来，我们创建目标函数的三维表面图：
surf_fig = plt.figure(tight_layout=True)
surf_ax = surf_fig.add_subplot(projection="3d")
surf_ax.tick_params(axis="both", which="major", labelsize=9)
surf_ax.set(xlabel="x", ylabel="y", zlabel="z")
surf_ax.set_title("Objective function")
surf_ax.plot_surface(x, y, z, alpha=0.7)
- 在我们开始最小化过程之前，我们需要定义一个初始点
x0。我们在前一步中创建的目标函数图上绘制这一点：
x0 = np.array([-0.8, 1.3])
surf_ax.plot([x0[0]], [x0[1]], func(x0), "r*")
目标函数的表面图，以及初始值，可以在以下图像中看到：
图 9.5：具有初始位置的目标函数表面
- 我们的
descend例程需要一个评估目标函数梯度的函数，因此我们将定义一个：
def grad(x):
c1 = x[0]**2 - x[0] + x[1]**2 + x[1] + 0.5
cos_t = np.cos(0.5*x[0]*x[1])
sin_t = np.sin(0.5*x[0]*x[1])
return np.array([
(2*x[0]-1)*cos_t - 0.5*x[1]*c1*sin_t,
(2*x[1]+1)*cos_t - 0.5*x[0]*c1*sin_t
])
- 我们将在轮廓图上绘制迭代，因此我们设置如下：
cont_fig, cont_ax = plt.subplots()
cont_ax.set(xlabel="x", ylabel="y")
cont_ax.set_title("Contour plot with iterates")
cont_ax.contour(x, y, z, levels=30)
- 现在，我们创建一个变量，将x和y方向的边界作为元组的元组保存。这些是步骤 10中
linspace调用中的相同边界：
bounds = ((-1, 1), (-2, 2))
- 我们现在可以使用
for循环驱动descend生成器来产生每个迭代，并将步骤添加到轮廓图中：
xnm1 = x0
for i, xn, fxn, grad_xn in descend(func, x0, grad, bounds):
cont_ax.plot([xnm1[0], xn[0]], [xnm1[1], xn[1]], "k*--")
xnm1, grad_xnm1 = xn, grad_xn
- 循环完成后，我们将最终值打印到终端：
print(f"iterations={i}")
print(f"min val at {xn}")
print(f"min func value = {fxn}")
前面打印语句的输出如下：
iterations=37
min val at [ 0.49999999 -0.49999999]
min func value = 2.1287163880894953e-16
在这里，我们可以看到我们的例程使用了 37 次迭代来找到大约在(0.5，-0.5)处的最小值，这是正确的。
可以在以下图像中看到带有其迭代的轮廓图：
图 9.6：梯度下降迭代到最小值的目标函数轮廓图
在这里，我们可以看到每次迭代的方向 - 由虚线表示 - 是目标函数下降最快的方向。最终迭代位于目标函数“碗”的中心，这是最小值出现的地方。
它是如何工作的…
这个配方的核心是descend例程。在这个例程中定义的过程是梯度下降法的一个非常简单的实现。在给定点计算梯度由grad参数处理，然后用direction = -grad推断迭代的旅行方向。我们将这个方向乘以一个固定的比例因子（有时称为学习率）0.2 的值来获得缩放步骤，然后通过添加0.2*direction到当前位置来进行这一步。
这个配方的解决方案需要 37 次迭代才能收敛，这比最小化非线性函数配方中的 Nelder-Mead 单纯形算法需要 58 次迭代要好一点。（这并不是一个完美的比较，因为我们改变了这个配方的起始位置。）这种性能在很大程度上取决于我们选择的步长。在这种情况下，我们将最大步长固定为方向向量的 0.2 倍。这使得算法简单，但并不特别高效。
在这个配方中，我们选择将算法实现为生成器函数，以便我们可以看到每一步的输出，并在迭代中绘制在等高线图上。在实践中，我们可能不想这样做，而是在迭代完成后返回计算出的最小值。为此，我们可以简单地删除yield语句，并在函数的最后，即主函数的缩进位置（即不在循环内），用return xn替换它。如果你想防止不收敛，你可以使用for循环的else特性来捕捉循环完成的情况，因为它已经到达了迭代器的末尾而没有触发break关键字。这个else块可以引发一个异常，以表明算法未能稳定到一个解决方案。在这个配方中，我们用来结束迭代的条件并不能保证该方法已经达到了最小值，但这通常是这种情况。
还有更多…
在实践中，你通常不会为自己实现梯度下降算法，而是使用来自库（如 SciPy optimize模块）的通用例程。我们可以使用与前一个配方中使用的相同的minimize例程来执行多种不同算法的最小化，包括几种梯度下降算法。这些实现可能比这样的自定义实现具有更高的性能和更强大的鲁棒性。
我们在这个配方中使用的梯度下降方法是一个非常天真的实现，可以通过允许例程在每一步选择步长来大大改进。（允许选择自己步长的方法有时被称为自适应方法。）这种改进的困难部分是选择在这个方向上采取的步长大小。为此，我们需要考虑单变量函数，它由以下方程给出：
在这里，x*[n]表示当前点，d[n]表示当前方向，t是一个参数。为了简单起见，我们可以使用 SciPy optimize模块中的minimize_scalar最小化例程来处理标量值函数。不幸的是，这并不像简单地传入这个辅助函数并找到最小值那样简单。我们必须限定t的可能值，以便计算出的最小化点x[n]+ td*[n]位于我们感兴趣的区域内。
要理解我们如何限制t的值，我们必须首先从几何上看这个构造。我们引入的辅助函数在给定方向上沿着一条直线评估目标函数。我们可以将其想象为通过当前x*[n]点在d[n]方向上穿过的表面的单个横截面。算法的下一步是找到最小化沿着这条直线的目标函数值的步长t*，这是一个标量函数，要最小化它要容易得多。然后边界应该是t值的范围，在这个范围内，这条直线位于由x和y边界值定义的矩形内。我们确定这条直线穿过这些x和y边界线的四个值，其中两个将是负值，另外两个将是正值。（这是因为当前点必须位于矩形内。）我们取两个正值中的最小值和两个负值中的最大值，并将这些边界传递给标量最小化例程。这是通过以下代码实现的：
alphas = np.array([
(bounds[0][0] - xn[0]) / direction[0], # x lower
(bounds[1][0] - xn[1]) / direction[1], # y lower
(bounds[0][1] - xn[0]) / direction[0], # x upper
(bounds[1][1] - xn[1]) / direction[1] # y upper
])
alpha_max = alphas[alphas >= 0].min()
alpha_min = alphas[alphas < 0].max()
result = minimize_scalar(lambda t: func(xn + t*direction),
method="bounded", bounds=(alpha_min, alpha_max))
amount = result.x
一旦选择了步长，唯一剩下的步骤就是更新当前的xn值，如下所示：
xn = xn + amount * direction
使用这种自适应步长会增加例程的复杂性，但性能大大提高。使用这个修改后的例程，该方法在只进行了三次迭代后就收敛了，远少于本食谱中使用的朴素代码（37 次迭代）或上一个食谱中使用的 Nelder-Mead 单纯形算法（58 次迭代）。通过在梯度函数中提供更多信息，减少了迭代次数，这正是我们预期的。
我们创建了一个函数，返回给定点处函数的梯度。我们在开始之前手动计算了这个梯度，这并不总是容易或甚至可能的。相反，更常见的是用数值计算的梯度替换这里使用的“解析”梯度，这个数值梯度是使用有限差分或类似算法估计的。这对性能和准确性有影响，就像所有近似一样，但鉴于梯度下降方法提供的收敛速度的提高，这些问题通常是次要的。
梯度下降类型的算法在机器学习应用中特别受欢迎。大多数流行的 Python 机器学习库，包括 PyTorch、TensorFlow 和 Theano，都提供了用于自动计算数据数组梯度的实用工具。这使得梯度下降方法可以在后台使用以提高性能。
梯度下降方法的一种流行变体是随机梯度下降，其中梯度是通过随机抽样来估计的，而不是使用整个数据集。这可以显著减少方法的计算负担，但会以更慢的收敛速度为代价，特别是对于高维问题，比如在机器学习应用中常见的问题。随机梯度下降方法通常与反向传播结合，构成了机器学习应用中训练人工神经网络的基础。
基本随机梯度下降算法有几种扩展。例如，动量算法将先前的增量纳入到下一个增量的计算中。另一个例子是自适应梯度算法，它纳入每个参数的学习率，以改善涉及大量稀疏参数的问题的收敛速度。
使用最小二乘法拟合数据曲线
最小二乘法是一种强大的技术，用于从相对较小的潜在函数族中找到最能描述特定数据集的函数。这种技术在统计学中特别常见。例如，最小二乘法用于线性回归问题-在这里，潜在函数族是所有线性函数的集合。通常，我们尝试拟合的这个函数族具有相对较少的参数，可以调整以解决问题。
最小二乘法的思想相对简单。对于每个数据点，我们计算残差的平方-点的值与给定函数的期望值之间的差异-并尝试使这些残差的平方和尽可能小（因此最小二乘法）。
在这个配方中，我们将学习如何使用最小二乘法来拟合一组样本数据的曲线。
做好准备
对于这个配方，我们将需要导入 NumPy 包，通常作为np，并导入 Matplotlib pyplot模块作为plt：
import numpy as np
import matplotlib.pyplot as plt
我们还需要从 NumPy random模块中导入默认随机数生成器的实例，如下所示：
from numy.random import default_rng
rng = default_rng(12345)
最后，我们需要从 SciPy optimize模块中的curve_fit例程：
from scipy.optimize import curve_fit
如何做…
以下步骤向您展示如何使用curve_fit例程来拟合一组数据的曲线：
- 第一步是创建样本数据：
SIZE = 100
x_data = rng.uniform(-3.0, 3.0, size=SIZE)
noise = rng.normal(0.0, 0.8, size=SIZE)
y_data = 2.0*x_data**2 - 4*x_data + noise
- 接下来，我们生成数据的散点图，以查看是否可以识别数据中的潜在趋势：
fig, ax = plt.subplots()
ax.scatter(x_data, y_data)
ax.set(xlabel="x", ylabel="y", title="Scatter plot of sample data")
我们生成的散点图可以在下面的图像中看到。在这里，我们可以看到数据显然不遇线性趋势（直线）。由于我们知道趋势是多项式，我们的下一个猜测将是二次趋势。这就是我们在这里使用的：
图 9.7：样本数据的散点图。我们可以看到数据不遵循线性趋势
- 接下来，我们创建一个代表我们希望拟合的模型的函数：
def func(x, a, b, c):
return a*x**2 + b*x + c
- 现在，我们可以使用
curve_fit例程将模型函数拟合到样本数据中：
coeffs, _ = curve_fit(func, x_data, y_data)
print(coeffs)
# [ 1.99611157 -3.97522213 0.04546998]
- 最后，我们在散点图上绘制最佳拟合曲线，以评估拟合曲线描述数据的效果如何：
x = np.linspace(-3.0, 3.0, SIZE)
y = func(x, coeffs[0], coeffs[1], coeffs[2])
ax.plot(x, y, "k--")
更新后的散点图可以在下面的图像中看到：
图 9.8：散点图，使用最小二乘法找到的最佳拟合曲线
在这里，我们可以看到我们找到的曲线相当合理地拟合了数据。
它是如何工作的…
curve_fit例程执行最小二乘拟合，将模型的曲线拟合到样本数据中。在实践中，这相当于最小化以下目标函数：
这里，配对（x[i]，y[i]）是样本数据中的点。在这种情况下，我们正在优化一个三维参数空间，每个参数都有一个维度。该例程返回估计的系数-参数空间中的点，其中目标函数被最小化-和包含拟合的协方差矩阵的估计的第二个变量。在这个配方中，我们忽略了这一点。
从curve_fit例程返回的估计协方差矩阵可以用于为估计的参数提供置信区间。这是通过取对角线元素的平方根除以样本大小（在这个配方中为 100）来完成的。这给出了估计的标准误差，当乘以与置信度对应的适当值时，给出了置信区间的大小。（我们在第六章中讨论了置信区间，使用数据和统计。）
您可能已经注意到，curve_fit例程估计的参数接近，但并非完全等于我们用于定义步骤 1中的样本数据的参数。这些参数并非完全相等的原因是我们向数据中添加了正态分布的噪声。在这个示例中，我们知道数据的基本结构是二次的 - 也就是二次多项式 - 而不是其他更神秘的函数。实际上，我们不太可能知道数据的基本结构，这就是我们向样本添加噪声的原因。
还有更多…
SciPy optimize模块中还有另一个用于执行最小二乘拟合的例程，名为least_squares。这个例程的签名略微不太直观，但确实返回了一个包含有关优化过程的更多信息的结果对象。然而，这个例程的设置方式可能更类似于我们在*它是如何工作的…*部分中构造基础数学问题的方式。要使用这个例程，我们定义目标函数如下：
def func(params, x, y):
return y - (params[0]*x**2 + params[1]*x + params[2])
我们将这个函数与参数空间中的起始估计x0一起传递，例如(1, 0, 0)。目标函数的额外参数func可以使用args关键字参数传递；例如，我们可以使用args=(x_data, y_data)。这些参数被传递到目标函数的x和y参数中。总之，我们可以使用以下调用least_squares来估计参数：
results = least_squares(func, [1, 0, 0], args=(x_data, y_data))
从least_squares例程返回的results对象实际上与本章中描述的其他优化例程返回的对象相同。它包含诸如使用的迭代次数、过程是否成功、详细的错误消息、参数值以及最小值处的目标函数值等细节。
分析简单的两人游戏
博弈论是数学的一个分支，涉及决策和策略分析。它在经济学、生物学和行为科学中有应用。许多看似复杂的情况可以简化为一个相对简单的数学游戏，可以以系统的方式进行分析，找到“最优”解决方案。
博弈论中的一个经典问题是囚徒困境，其原始形式如下：两个同谋被抓住，必须决定是保持沉默还是对对方作证。如果两人都保持沉默，他们都要服刑 1 年；如果一个作证而另一个不作证，作证者将获释，而另一个将服刑 3 年；如果两人都相互作证，他们都将服刑 2 年。每个同谋应该怎么做？事实证明，在对对方有任何合理的不信任的情况下，每个同谋可以做出的最佳选择是作证。采用这种策略，他们将不会服刑或最多服刑 2 年。
由于这本书是关于 Python 的，我们将使用这个经典问题的变体来说明这个问题的普遍性。考虑以下问题：你和你的同事必须为客户编写一些代码。你认为你可以用 Python 更快地编写代码，但你的同事认为他们可以用 C 更快地编写代码。问题是，你应该选择哪种语言来进行项目？
你认为你可以用 Python 代码比 C 代码快 4 倍，所以你用速度 1 写 C，用速度 4 写 Python。你的同事说他们可以比 Python 稍微更快地写 C，所以他们用速度 3 写 C，用速度 2 写 Python。如果你们两个都同意一种语言，那么你们按照你预测的速度编写代码，但如果你们意见不一致，那么更快的程序员的生产力会降低 1。我们可以总结如下：
| 同事/你 | C | Python |
|---|---|---|
| C | 3 / 1 | 3 / 2 |
| Python | 2 / 1 | 2 / 4 |
在这个配方中，我们将学习如何在 Python 中构建一个对象来表示这个简单的双人游戏，然后对这个游戏的结果进行一些基本分析。
准备工作
对于这个配方，我们需要导入 NumPy 包为np，导入 Nashpy 包为nash：
import numpy as np
import nashpy as nash
如何做…
以下步骤向您展示了如何使用 Nashpy 创建和执行一些简单的双人游戏分析：
- 首先，我们需要创建矩阵，用于保存每个玩家（在这个例子中是您和您的同事）的支付信息。
you = np.array([[1, 3], [1, 4]])
colleague = np.array([[3, 2], [2, 2]])
- 接下来，我们创建一个
Game对象，它保存了由这些支付矩阵表示的游戏：
dilemma = nash.Game(you, colleague)
- 我们使用索引表示法计算给定选择的效用：
print(dilemma[[1, 0], [1, 0]]) # [1 3]
print(dilemma[[1, 0], [0, 1]]) # [3 2]
print(dilemma[[0, 1], [1, 0]]) # [1 2]
print(dilemma[[0, 1], [0, 1]]) # [4 2]
- 我们还可以根据做出特定选择的概率计算预期效用：
print(dilemma[[0.1, 0.9], [0.5, 0.5]]) # [2.45 2.05]
它是如何工作的…
在这个配方中，我们构建了一个代表非常简单的双人战略游戏的 Python 对象。 这里的想法是有两个“玩家”需要做决定，每个玩家的选择组合都会给出一个特定的支付值。 我们的目标是找到每个玩家可以做出的最佳选择。 假设玩家同时进行一次移动，即没有人知道对方的选择。 每个玩家都有一个确定他们所做选择的策略。
在步骤 1中，我们创建两个矩阵 - 每个玩家一个 - 分配给每个选择组合的支付值。 这两个矩阵由 Nashpy 的Game类包装，它提供了一个方便且直观（从博弈论的角度来看）的接口来处理游戏。 通过使用索引表示法传递选择，我们可以快速计算给定选择组合的效用。
我们还可以根据一种策略提供预期效用的计算，其中选择是根据某种概率分布随机选择的。 语法与先前描述的确定性情况相同，只是我们为每个选择提供了一个概率向量。 我们根据您选择 Python 的概率计算预期效用为 90％，而您的同事选择 Python 的概率为 50％。 预期速度分别为 2.45 和 2.05。
还有更多…
在 Python 中，还有一种计算博弈论的替代方法。 Gambit 项目是一组用于博弈论计算的工具，具有 Python 接口（www.gambit-project.org/）。 这是一个成熟的项目，建立在 C 库周围，并提供比 Nashpy 更高的性能。
计算纳什均衡
纳什均衡是一个双人战略游戏 - 类似于我们在分析简单的双人游戏配方中看到的游戏 - 代表每个玩家都看到“最佳可能”结果的“稳态”。 但是，这并不意味着与纳什均衡相关联的结果是最好的。 纳什均衡比这更微妙。 纳什均衡的非正式定义如下：在其中没有个别玩家可以改善他们的结果的行动配置，假设所有其他玩家都遵守该配置。
我们将探讨纳什均衡的概念，使用经典的猜拳游戏。 规则如下。 每个玩家可以选择以下选项之一：石头，纸或剪刀。 石头打败剪刀，但输给纸； 纸打败石头，但输给剪刀； 剪刀打败纸，但输给石头。 如果两名玩家做出相同选择的游戏是平局。 在数值上，我们用+1 表示赢，-1 表示输，0 表示平局。 从中，我们可以构建一个双人游戏，并计算该游戏的纳什均衡。
在这个配方中，我们将计算猜拳游戏的纳什均衡。
准备工作
对于这个配方，我们需要导入 NumPy 包为np，导入 Nashpy 包为nash：
import numpy as np
import nashpy as nash
如何做…
以下步骤向您展示了如何计算简单的两人游戏的纳什均衡：
- 首先，我们需要为每个玩家创建一个收益矩阵。我们将从第一个玩家开始：
rps_p1 = np.array([
[ 0, -1, 1], # rock payoff
[ 1, 0, -1], # paper payoff
[-1, 1, 0] # scissors payoff
])
- 第二个玩家的收益矩阵是
rps_p1的转置：
rps_p2 = rps_p1.transpose()
- 接下来，我们创建
Game对象来表示游戏：
rock_paper_scissors = nash.Game(rps_p1, rps_p2)
- 我们使用支持枚举算法计算游戏的纳什均衡：
equilibria = rock_paper_scissors.support_enumeration()
- 我们遍历均衡，并打印每个玩家的策略：
for p1, p2 in equilibria:
print("Player 1", p1)
print("Player 2", p2)
这些打印语句的输出如下：
Player 1 [0.33333333 0.33333333 0.33333333]
Player 2 [0.33333333 0.33333333 0.33333333]
它是如何工作的…
纳什均衡在博弈论中非常重要，因为它们允许我们分析战略游戏的结果，并确定有利的位置。它们最早由约翰·F·纳什在 1950 年描述，并在现代博弈论中发挥了重要作用。一个两人游戏可能有许多纳什均衡，但任何有限的两人游戏必须至少有一个。问题在于找到给定游戏的所有可能的纳什均衡。
在这个示例中，我们使用了支持枚举，它有效地枚举了所有可能的策略，并筛选出了那些是纳什均衡的策略。在这个示例中，支持枚举算法只找到了一个纳什均衡，这是一个混合策略。这意味着唯一没有改进的策略是随机选择其中一个选项，每个选项的概率为 1/3。这对于任何玩过石头剪刀布的人来说都不足为奇，因为无论我们做出什么选择，我们的对手都有 1/3 的机会（随机）选择打败我们选择的动作。同样，我们有 1/3 的机会平局或赢得比赛，因此我们在所有这些可能性中的期望值如下：
在不知道我们的对手会选择哪一个选项的情况下，没有办法改进这个预期结果。
还有更多…
Nashpy 软件包还提供了其他计算纳什均衡的算法。具体来说，vertex_enumeration方法在Game对象上使用顶点枚举算法，而lemke_howson_enumeration方法使用Lemke Howson算法。这些替代算法具有不同的特性，对于某些问题可能更有效。
另请参阅
Nashpy 软件包的文档包含了有关所涉及的算法和博弈论的更详细信息。其中包括了一些关于博弈论的文本参考。这些文档可以在nashpy.readthedocs.io/en/latest/找到。
进一步阅读
通常情况下，Numerical Recipes书是数值算法的良好来源。第十章，其他主题，涉及函数的最大化和最小化：
- Press, W.H., Teukolsky, S.A., Vetterling, W.T., and Flannery, B.P., 2017. Numerical recipes: the art of scientific computing. 3rd ed. Cambridge: Cambridge University Press。
有关优化的更具体信息可以在以下书籍中找到：
Boyd, S.P. and Vandenberghe, L., 2018. Convex optimization*. Cambridge: Cambridge University Press*。
*Griva, I., Nash, S., and Sofer, A., 2009. Linear and nonlinear optimization.*2nd ed. Philadelphia: Society for Industrial and Applied Mathematics。
最后，以下书籍是博弈论的良好入门：
- Osborne, M.J., 2017. An introduction to game theory*. Oxford: Oxford University Press*。