Python 数学应用（三）（1）

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

python3.8 -m pip install statsmodels sklearn fbprophet

conda install fbprophet

import statsmodels.api as sm
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import default_rng
rng = default_rng(12345)

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")

pred_x = sm.add_constant(x)

model1 = sm.OLS(y1, pred_x).fit()
print(model1.summary())

model2 = sm.OLS(y2, pred_x).fit()
print(model2.summary())

model_x = sm.add_constant(np.linspace(0, 5))

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')

from numpy.random import default_rng
rng = default_rng(12345)

import statsmodels.api as sm

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:
          ![](https://gitee.com/OpenDocCN/freelearn-python-zh/raw/master/docs/app-math-py/img/151a3a34-831e-40fc-9d29-ec1db98ab665.png)
        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`:

The output of the `print` statement is as follows:

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:

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:
          ![](https://gitee.com/OpenDocCN/freelearn-python-zh/raw/master/docs/app-math-py/img/0839eae4-836a-497f-9f90-8dd8f2cd4a17.png)
        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:

We also need several components from the `scikit-learn` package to perform logistic regression. These can be imported as follows:

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:

2.  Now, we use two of our three predictor variables to create our response variable as a series of Boolean values:

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:

The resulting plot can be seen in the following figure:
          ![](https://gitee.com/OpenDocCN/freelearn-python-zh/raw/master/docs/app-math-py/img/e972cdff-197a-4940-8f33-dbe849ddb774.png)
        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:

5.  Next, we prepare some extra data, different from what we used to fit the model, to test the accuracy of our model:

6.  Then, we generate predicted results based on our logistic regression model:

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:

The report that's generated by this routine looks as follows:

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:
          ![](https://gitee.com/OpenDocCN/freelearn-python-zh/raw/master/docs/app-math-py/img/ef09118a-1d49-472f-9b70-c9925e410036.png)
        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:
          ![](https://gitee.com/OpenDocCN/freelearn-python-zh/raw/master/docs/app-math-py/img/4edf13c3-10f1-42a0-a475-91264d1ea732.png)
        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:

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:

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:

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:
          ![](https://gitee.com/OpenDocCN/freelearn-python-zh/raw/master/docs/app-math-py/img/518372fa-7c61-409c-90aa-f0a13459785f.png)
        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:

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:

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:
          ![](https://gitee.com/OpenDocCN/freelearn-python-zh/raw/master/docs/app-math-py/img/093fb316-c237-48f3-b47b-a62afcde7afa.png)
        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:

6.  Now, we fit the model to the data and get the resulting model. We print a summary of these results to the Terminal:

The summary data given for the fitted model is as follows:

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:

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:
          ![](https://gitee.com/OpenDocCN/freelearn-python-zh/raw/master/docs/app-math-py/img/6dc54aaa-0b81-403f-b4d8-2a67cb360f1e.png)
        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*:

The updated plot can be seen in the following figure:
          ![](https://gitee.com/OpenDocCN/freelearn-python-zh/raw/master/docs/app-math-py/img/67b312d7-3986-4bf8-aba6-a534a0cf9b89.png)
        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:
          ![](https://gitee.com/OpenDocCN/freelearn-python-zh/raw/master/docs/app-math-py/img/7a4a2e01-301a-4448-9ad2-32ab71c37a97.png)
        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:
          ![](https://gitee.com/OpenDocCN/freelearn-python-zh/raw/master/docs/app-math-py/img/c71efac3-1651-4fc9-932c-e3930e7d80e7.png)
        Here, the variants of *θ[i]* are parameters. Putting these two models together gives us an ARMA(1, 1) model, which has the following form:
          ![](https://gitee.com/OpenDocCN/freelearn-python-zh/raw/master/docs/app-math-py/img/fcb7eb3b-9d02-4a2e-b954-b260ba87f8e0.png)
        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:

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:

2.  As usual, the next step is to plot the time series so that we can visually identify the trend of the data:

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:
          ![](https://gitee.com/OpenDocCN/freelearn-python-zh/raw/master/docs/app-math-py/img/a5081c24-e36c-48ac-86e5-a5ec2cfccee0.png)
        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:

4.  Now, we plot the ACF and PACF for the differenced time series: 

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:
          ![](https://gitee.com/OpenDocCN/freelearn-python-zh/raw/master/docs/app-math-py/img/dc3a6b97-bdd4-4c6c-822a-d40b00a0a2e8.png)
        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:

The summary information that's printed looks as follows: