TowardsDataScience 博客中文翻译 2020（八百九十）

原文：TowardsDataScience Blog

协议：CC BY-NC-SA 4.0

用 Python 中的 SARIMA 进行时间序列预测

原文：https://towardsdatascience.com/time-series-forecasting-with-sarima-in-python-cda5b793977b?source=collection_archive---------2-----------------------

关于使用 Python 使用 SARIMA 进行时间序列建模的实践教程

摩根·豪斯尔在 Unsplash 上的照片

在之前的文章中，我们介绍了移动平均过程 MA(q) ，和自回归过程 AR§ 。我们将它们结合起来，形成了 ARMA(p，q)和 ARIMA(p，d，q)模型来模拟更复杂的时间序列。

现在，给模型添加最后一个组件:季节性。

本文将涵盖:

• 季节性 ARIMA 模型
• 使用真实数据的完整建模和预测项目

Github 上有笔记本和数据集。

我们开始吧！

关于 Python 中时间序列分析的完整课程，涵盖统计和深度学习模型，请查看我新发布的 课程 ！

萨里玛模型

到目前为止，我们还没有考虑时间序列中季节性的影响。然而，这种行为肯定存在于许多情况下，如礼品店销售或飞机乘客总数。

季节性 ARIMA 模型或萨里玛的写法如下:

萨里玛符号

你可以看到我们为时间序列的季节性部分添加了 P、D 和 Q。它们与非季节性成分的术语相同，因为它们涉及季节性周期的后移。

在上面的公式中， m 是每年或一段时间内的观察次数。如果我们分析季度数据， m 等于 4。

ACF 和 PACF 图

AR 和 MA 模型的季节性部分可以从 PACF 和 ACF 图中推断出来。

在 SARIMA 模型中，只有一个阶为 1、周期为 12 的季节性移动平均过程，表示为:

• 在滞后 12 时观察到峰值
• PACF 季节滞后的指数衰减(滞后 12，24，36，…)

类似地，对于仅具有阶为 1、周期为 12 的季节性自回归过程的模型:

• ACF 季节性滞后的指数衰减(滞后 12，24，36，…)
• 在 PACF 中，在滞后 12 时观察到峰值

系统模型化

建模过程与非季节性 ARIMA 模型相同。在这种情况下，我们只需要考虑额外的参数。

使时间序列平稳和根据最低 AIC 选择模型所需的步骤仍留在建模过程中。

让我们用真实世界的数据集来涵盖一个完整的示例。

项目—为强生公司的季度每股收益建模

我们将重温强生公司的季度每股收益(EPS)数据集。这是一个非常有趣的数据集，因为有一个移动平均过程在起作用，而且我们有季节性，这是用 SARIMA 进行一些练习的完美时机。

和往常一样，我们首先导入所有必要的库来进行分析

from statsmodels.graphics.tsaplots import plot_pacf
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from statsmodels.tsa.stattools import adfuller
import matplotlib.pyplot as plt
from tqdm import tqdm_notebook
import numpy as np
import pandas as pdfrom itertools import productimport warnings
warnings.filterwarnings('ignore')%matplotlib inline

现在，让我们读入数据帧中的数据:

data = pd.read_csv('jj.csv')

然后，我们可以显示一个时间序列图:

plt.figure(figsize=[15, 7.5]); # Set dimensions for figure
plt.plot(data['date'], data['data'])
plt.title('Quarterly EPS for Johnson & Johnson')
plt.ylabel('EPS per share ($)')
plt.xlabel('Date')
plt.xticks(rotation=90)
plt.grid(True)
plt.show()

强生公司季度每股收益

很明显，时间序列不是静态的，因为它的平均值在整个时间内不是恒定的，我们看到数据中的方差在增加，这是异方差的标志。

为了确保这一点，让我们绘制 PACF 和 ACF:

plot_pacf(data['data']);
plot_acf(data['data']);

PACF 和 ACF

同样，不能从这些图中推断出任何信息。您可以使用扩展的 Dickey-Fuller 测试进一步测试平稳性:

ad_fuller_result = adfuller(data['data'])
print(f'ADF Statistic: {ad_fuller_result[0]}')
print(f'p-value: {ad_fuller_result[1]}')

ADF 测试结果

由于 p 值很大，我们不能拒绝零假设，必须假设时间序列是非平稳的。

现在，让我们取对数差，努力使它保持稳定:

data['data'] = np.log(data['data'])
data['data'] = data['data'].diff()
data = data.drop(data.index[0])

绘制新数据应给出:

plt.figure(figsize=[15, 7.5]); # Set dimensions for figure
plt.plot(data['data'])
plt.title("Log Difference of Quarterly EPS for Johnson & Johnson")
plt.show()

强生公司季度每股收益记录表

厉害！现在，我们仍然在上面的图中看到季节性。因为我们处理的是季度数据，所以我们的周期是 4。因此，我们将在 4:

# Seasonal differencingdata['data'] = data['data'].diff(4)
data = data.drop([1, 2, 3, 4], axis=0).reset_index(drop=True)

绘制新数据:

plt.figure(figsize=[15, 7.5]); # Set dimensions for figure
plt.plot(data['data'])
plt.title("Log Difference of Quarterly EPS for Johnson & Johnson")
plt.show()

完美！请记住，虽然我们在 4 个月的时间内取了差值，但季节差值(D)的顺序是 1，因为我们只取了一次差值。

现在，让我们再次运行扩展的 Dickey-Fuller 测试，看看我们是否有一个平稳的时间序列:

ad_fuller_result = adfuller(data['data'])
print(f'ADF Statistic: {ad_fuller_result[0]}')
print(f'p-value: {ad_fuller_result[1]}')

事实上，p 值足够小，我们可以拒绝零假设，我们可以认为时间序列是平稳的。

看看 ACF 和 PACF:

plot_pacf(data['data']);
plot_acf(data['data']);

我们可以从 PACF 中看到，我们在滞后 1 处有一个显著的峰值，这表明 AR(1)过程。此外，我们在滞后 4 处有另一个峰值，表明 1 阶的季节性自回归过程(P = 1)。

查看 ACF 图，我们仅在滞后 1 处看到一个显著的峰值，表明存在非季节性 MA(1)过程。

尽管这些图可以让我们大致了解正在进行的过程，但最好测试多个场景，并选择产生最低 AIC 的模型。

因此，让我们编写一个函数来测试 SARIMA 模型的一系列参数，并输出一个性能最佳的模型位于顶部的表:

def optimize_SARIMA(parameters_list, d, D, s, exog):
    """
        Return dataframe with parameters, corresponding AIC and SSE

        parameters_list - list with (p, q, P, Q) tuples
        d - integration order
        D - seasonal integration order
        s - length of season
        exog - the exogenous variable
    """

    results = []

    for param in tqdm_notebook(parameters_list):
        try: 
            model = SARIMAX(exog, order=(param[0], d, param[1]), seasonal_order=(param[2], D, param[3], s)).fit(disp=-1)
        except:
            continue

        aic = model.aic
        results.append([param, aic])

    result_df = pd.DataFrame(results)
    result_df.columns = ['(p,q)x(P,Q)', 'AIC']
    #Sort in ascending order, lower AIC is better
    result_df = result_df.sort_values(by='AIC', ascending=True).reset_index(drop=True)

    return result_df

请注意，我们将只测试参数 P、P、q 和 q 的不同值。我们知道季节和非季节积分参数都应该是 1，季节长度是 4。

因此，我们生成所有可能的参数组合:

p = range(0, 4, 1)
d = 1
q = range(0, 4, 1)
P = range(0, 4, 1)
D = 1
Q = range(0, 4, 1)
s = 4parameters = product(p, q, P, Q)
parameters_list = list(parameters)
print(len(parameters_list))

你应该看到我们有 256 种不同的组合！现在，我们的函数将根据我们的数据拟合 256 种不同的 SARIMA 模型，以找到 AIC 最低的模型:

result_df = optimize_SARIMA(parameters_list, 1, 1, 4, data['data'])
result_df

结果表

从表中可以看出，最佳模型是:SARIMA(0，1，2)(0，1，2，4)。

我们现在可以拟合模型并输出其摘要:

best_model = SARIMAX(data['data'], order=(0, 1, 2), seasonal_order=(0, 1, 2, 4)).fit(dis=-1)
print(best_model.summary())

最佳模式总结

在这里，你可以看到表现最好的模型同时具有季节性和非季节性移动平均过程。

根据上面的总结，您可以找到系数的值及其 p 值。请注意，从 p 值来看，所有系数都是显著的。

现在，我们可以研究残差:

best_model.plot_diagnostics(figsize=(15,12));

模型诊断

从正常的 Q-Q 图，我们可以看到，我们几乎有一条直线，这表明没有系统偏离正常。此外，右下角的相关图表明残差中没有自相关，因此它们实际上是白噪声。

我们已经准备好绘制模型的预测图，并预测未来:

data['arima_model'] = best_model.fittedvalues
data['arima_model'][:4+1] = np.NaNforecast = best_model.predict(start=data.shape[0], end=data.shape[0] + 8)
forecast = data['arima_model'].append(forecast)plt.figure(figsize=(15, 7.5))
plt.plot(forecast, color='r', label='model')
plt.axvspan(data.index[-1], forecast.index[-1], alpha=0.5, color='lightgrey')
plt.plot(data['data'], label='actual')
plt.legend()plt.show()

模型预测

瞧啊。

结论

恭喜你！您现在了解了什么是季节性 ARIMA(或 SARIMA)模型，以及如何使用它进行建模和预测。

通过以下课程了解有关时间序列的更多信息:

• 在 Python 中应用时间序列分析

干杯🍺

使用 SQL 进行时间序列预测——比您想象的要简单

原文：https://towardsdatascience.com/time-series-forecasting-with-sql-its-easier-than-you-think-1f5b362d0c81?source=collection_archive---------9-----------------------

是的，SQL 现在可以做到这一点。

时间序列预测是我通常用 Python 做的一项任务。你可能已经习惯了其他语言，比如 R 或者 Julia ，但是我敢打赌你从来没有想到过这种类型的任务。如果是这样的话——继续读下去——你会惊讶地发现只用 SQL 就能完成这么多事情。

晨酿在 Unsplash 上拍摄的照片

我以前写过关于使用 SQL 执行分类任务的文章，所以如果你对此感兴趣的话，一定要看看:

[## 使用 SQL 进行机器学习——比你想象的要简单

当 Python 不是一个选项时该做什么。包括代码。

towardsdatascience.com](/machine-learning-with-sql-its-easier-than-you-think-c6aae9064d5a)

时间序列不同于一般的机器学习任务。你不可能训练一次模型，在生产中用几个月。时间序列模型必须使用完整的历史数据进行训练，并且新的数据点可能每小时、每天、每周或每月都会出现—因项目而异。

这就是为什么在硬件资源有限的情况下，在数据库中进行训练是有益的。 Python 几乎总是会比数据库消耗更多的资源。

我们将再次使用甲骨文云。这是免费的，所以请注册并创建一个 OLTP 数据库实例(版本 19c，有 0.2TB 的存储)。完成后，下载云钱包并通过SQL Developer——或任何其他工具建立连接。

这将花费你至少 10 分钟，但这是一件相当简单的事情，所以我不会在这上面浪费时间。

厉害！让我们继续数据加载。

数据加载

在进行任何类型的预测之前，我们需要一些数据。任何时间序列教程的事实上的标准数据集是航空乘客数据集。下载它，并在安全的地方保存一分钟。

我们需要创建一个保存数据集的表，接下来我们就这么做。下面是 SQL 语句:

CREATE TABLE airline_passengers(
    air_period DATE,
    air_passengers INTEGER
);

我们现在可以通过导入数据功能加载数据集:

当弹出一个模态窗口时，只需提供下载的 CSV 文件的路径，然后点击几次 Next 。使用您的最佳判断选择列，并选择日期格式为YYYY-MM。

完成后，我们的数据集就可以使用了:

厉害！我们现在可以继续进行模型训练和预测。

模特培训

我们的数据集有 144 行。我们不会对其整体进行模型训练。我们将保留最后 12 行用于评估。

从训练开始，我们需要创建一个指向训练数据的VIEW。方法如下:

CREATE OR REPLACE VIEW src_passengers AS
    SELECT * FROM airline_passengers 
    WHERE air_period < TO_DATE('1960-01-01', 'YYYY-MM-DD');

src_passengers视图现在保存了前 132 行——这正是我们想要的。

接下来，我们将声明一个处理模型训练的简短的PL/SQL片段:

DECLARE
    v_setlst DBMS_DATA_MINING.SETTING_LIST;
BEGIN
    v_setlst(DBMS_DATA_MINING.ALGO_NAME) := DBMS_DATA_MINING.ALGO_EXPONENTIAL_SMOOTHING;
    v_setlst(DBMS_DATA_MINING.EXSM_INTERVAL) := DBMS_DATA_MINING.EXSM_INTERVAL_MONTH;
    v_setlst(DBMS_DATA_MINING.EXSM_PREDICTION_STEP) := '12';
    v_setlst(DBMS_DATA_MINING.EXSM_MODEL) := DBMS_DATA_MINING.EXSM_HW;
    v_setlst(DBMS_DATA_MINING.EXSM_SEASONALITY) := '12';

    DBMS_DATA_MINING.CREATE_MODEL2(
        model_name => 'AIRLINE_PSG_FORECAST',
        mining_function => 'TIME_SERIES',
        data_query => 'SELECT * FROM src_passengers',
        set_list => v_setlst,
        case_id_column_name => 'air_period',
        target_column_name => 'air_passengers'
    );
END;
/

让我们把这个片段分解一下，让它更容易理解:

• DBMS_DATA_MINING.ALGO_NAME -时间序列预测算法类型，指数平滑是目前唯一可用的算法
• DBMS_DATA_MINING.EXSM_INTERVAL -表示数据集的间隔。我们的数据按月间隔存储，因此有了EXSM_INTERVAL_MONTH值
• DBMS_DATA_MINING.PREDICTION_STEP -做出多少预测。12(一年)不错
• DBMS_DATA_MINING.EXSM_MODEL -本质上是指数平滑模型的超参数组合。我选择使用三重指数平滑法或霍尔特-温特斯法。这里是可用算法的完整列表。
• DBMS_DATA_MINING.EXSM_SEASONALITY -表示一个季节持续多长时间

声明之后，我们可以在DBMS_DATA_MINING.CREATE_MODEL2过程的帮助下创建一个时间序列模型(顺便说一下，这是一个很好的命名约定)。以下是解释:

• model_name -任意，随心所欲命名型号
• mining_function -设置为TIME_SERIES，原因很清楚
• data_query -模型如何获得训练数据
• set_list -之前声明的设置列表，告诉 Oracle 如何实际训练模型
• case_id_column_name -包含日期值的列的名称
• target_column_name -包含数值的列的名称(我们试图预测的内容)

就是这样！如果你能理解这一点，你就知道如何用 SQL 训练时间序列模型。

您现在可以运行PL/SQL代码片段了。需要几秒钟才能完成。完成后，您可以继续下一部分。

模型评估

让我们看看我们的模型表现有多好。为此，我准备了以下声明:

SELECT
    a.case_id AS time_period,
    b.air_passengers AS actual,
    ROUND(a.prediction, 2) AS predicted,
    ROUND((b.air_passengers - a.prediction), 2) AS difference,
    ROUND((ABS(b.air_passengers - a.prediction) / b.air_passengers) * 100, 2) AS pct_error
FROM 
    dm$p0airline_psg_forecast a, 
    airline_passengers b
WHERE 
    a.case_id = b.air_period
AND a.case_id >= TO_DATE('1960-01-01', 'YYYY-MM-DD');

它将实际数据与 Holt-Winters 算法所做的预测进行比较，并比较绝对误差和百分比误差。以下是上述 SQL 语句的输出:

厉害！就我们投入的工作而言，我们的模型并没有那么糟糕。让我们在下一部分总结一下。

在你走之前

这工作量很大——不要争论了。尽管如此，我们仍有改进的方法。一个想法立即浮现在脑海中——创建一个将返回最佳算法的函数。

您可以将所有可能的算法存储在一个数组中，然后在遍历数组的同时训练一个模型，跟踪每个模型的性能。但这是另一个话题了。

感谢阅读。

加入我的私人邮件列表，获取更多有用的见解。

喜欢这篇文章吗？成为 中等会员 继续无限制学习。如果你使用下面的链接，我会收到你的一部分会员费，不需要你额外付费。

[## 通过我的推荐链接加入 Medium-Dario rade ci

作为一个媒体会员，你的会员费的一部分会给你阅读的作家，你可以完全接触到每一个故事…

medium.com](https://medium.com/@radecicdario/membership)

原载于 2020 年 9 月 22 日 https://betterdatascience.com。

Python 代码中基于统计模型的时间序列预测

原文：https://towardsdatascience.com/time-series-forecasting-with-statistical-models-in-python-code-da457a46d68a?source=collection_archive---------9-----------------------

在我之前的文章中，我展示了使用脸书先知 Python API(可用的统计模型之一)预测数字广告支出是多么容易。在第二部分中，我们将更多地讨论如何应用统计方法进行时间序列分析。

问题:根据 2 年多的历史每日广告支出，预测未来 2 个月的数字广告支出。

下面快速浏览一下数据文件(从 2017 年 1 月 1 日到 2019 年 9 月 23 日的每日广告支出)

创建培训和测试数据集:

目录

1. 介绍

• 时间序列分析
• 什么是流行的统计模型
• 单变量统计时间序列建模的关键要素

2.ARIMA

• ACF 和 PACF 图
• 萨里玛
• 案例研究:用 SARIMA 预测广告支出

3.美国教育考试服务中心

• 美国教育考试服务中心
• 霍尔特-温特季节方法
• 案例研究:用霍尔特-温特季节平滑法预测广告支出

4.对比模型

5.结案摘要

1.1 时间序列分析

细节在我之前的帖子中解释过，这里。

1.2 我们可以使用什么时间序列统计模型？

a.单变量统计时间序列建模(例如，平均和平滑模型、有/无季节项的 ARIMA 模型)

b.多变量统计时间序列建模(如外部回归变量、风险值)

c.附加或组件模型(如第一部分中涉及的脸书预言家)

d.结构时序建模(例如，贝叶斯结构时序建模、分层时序建模)。

在本文中，我将重点解释萨里玛和霍尔特-温特斯。

1.3 关键组件

我们需要我们的数据平稳性，这意味着(1)常数均值(2)常数方差(3)自协方差不依赖于时间，如果我们想使用时间序列数据的统计建模。

那如果我的数据不符合这个统计假设呢？换句话说，非平稳数据是什么样的？

非平稳数据背后的两个共同特征是:( 1)趋势~均值不随时间变化;( 2)季节性~方差不随时间变化;( 3)自协方差依赖于时间。

去趋势化和消除季节性的方法有哪些？

• 转换(例如，日志、sqrt 等)
• 平滑(例如，滚动平均等)
• 差异
• 分解
• 多项式拟合(例如，拟合回归模型)

2.1 ARIMA 模型中的关键术语

格式:ARIMA(p，d，q) — (1) p: AR 项(2) d: (3) q: MA 项

图片来自https://www.slideshare.net/21_venkat?

• AR(自回归)项将序列中的下一步描述为前一时间步观察值的线性函数。
• I(积分)项，序列的差分预处理步骤，使序列稳定。
• MA(移动平均)项将序列中的下一步描述为前一时间步平均过程的残差的线性函数。

模型的符号包括将 AR §、I (d)和 MA (q)模型的阶指定为 ARIMA 函数的参数，例如 ARIMA(p，d，q)。ARIMA 模型也可以用来开发 AR、MA 和 ARMA 模型。

那我们怎么确定 p，d，q 呢？

2.2 ACF 和 PACF 图

ACF(自相关函数)描述了时间序列与其自身滞后版本之间的相关性(例如，Y(t)与 Y(t-1)的相关性)

PACF(部分自相关函数)提供了由每个连续滞后项解释的附加相关性。

1. 用 PACF 图确定 p
2. 使用 ACF 图确定 q

[## 识别 ARIMA 模型中 AR 和 MA 项的顺序

非季节性 ARIMA 模型说明(pdf 文件)季节性和非季节性 ARIMA 模型幻灯片(pdf 文件)介绍…

people.duke.edu](https://people.duke.edu/~rnau/411arim3.htm)

2.3 萨里玛

SARIMA(季节性自回归综合移动平均)在季节性水平上结合了 ARIMA 模型。

格式:SARIMA(p，D，q) (P，D，Q)s 其中 s 代表季节周期。

2.3 案例研究:用 SARIMA 预测广告支出

我创建了 test _ stationarity 来检查平稳性。

该函数绘制了滚动平均值和标准差，并通过进行增强迪基-富勒测试输出 p 值。

检查原始数据:test _ stationary(df1[’ Spend ']。dropna())

与临界值相比，时间序列是非平稳的。

让我们做一个对数转换:df1[’ log _ Spend ‘]= NP . log(df1[’ Spend '])

test _ stationary(df1[’ log _ Spend ']。diff(1)。dropna())

太好了。利用对数变换，时间序列在 5%的阈值处是平稳的。

让我们尝试第一个差异:test _ stationary(df1[’ Spend ']。diff(1)。dropna())

令人印象深刻。结果是事件更好:时间序列在 1%的阈值处是稳定的，具有一阶差分。

建立 SARIMA 模型，预测 2019 年 7 月 23 日至 2019 年 9 月 23 日两个月的广告支出

现在，让我们通过从sk learn . metrics import mean _ squared _ error，mean_absolute_error 计算 mse 和 mae 来检查该模型的性能

让我们想象一下预测:

SARIMA 模型能够捕捉趋势和每周季节性，但在 MAE 为 5813 的波峰和波谷上有点偏离。

3.指数平滑(ETS)

因为时间序列数据在一段时间内自然是随机的，我们通常希望平滑数据，为此我们将使用 ETS，使用指数方法为过去的观察值分配较少的权重。在 wiki 上有一个像样的解释。

类似地，我们也希望通过应用 ETS 将时间序列分解成趋势(T)、季节(S)和不规则或误差(E)分量。

三种常见的 ETS 类型:

• **线性:**双指数平滑:
• **加法:**三重指数平滑
• **乘法:**也称三重指数平滑

来自:https://www . stat . Berkeley . edu/~ arturof/Teaching/stat 248/lab 10 _ part 1 . html

3.2 霍尔特-温特斯季节性方法

热冬季节法是三重指数平滑法的一种。它有三个关键部分:(1)一个水平(2)一个趋势(3)季节成分霍尔特-温特斯加法和霍尔特-温特斯乘法。

3.3 案例研究:用霍尔特-温特季节平滑法预测广告支出

建立热冬附加模型，预测 2019 年 7 月 23 日至 2019 年 9 月 23 日两个月的广告支出

通过从 sklearn.metrics 导入 mean_squared_error，mean_absolute_error 计算 mse 和 mae

H-W 模型能够捕捉趋势和每周季节性，但在 MAE 为 6055 的波峰和波谷上有点偏离。

4.比较这两个模型

让我们在同一个图中可视化两个模型的预测。

从图中，我们可以看出 SARIMA 模型比 Holt_Winter 模型表现稍好，MSE 和 MAE 都较低。尽管这两个模型都不能完美地捕捉到波峰和波谷，但在这种情况下，它们仍然对业务有用。根据数据，平均每月广告支出为 200 万美元以上。然而，基于这两种模型的平均寿命约为 6000。换句话说，对于一个平均每月广告支出为 200 万美元的企业来说，两个月的广告支出预测将下降约 6000 美元。这是一个相当不错的预测！

5.结案摘要

在本文中，单变量预测方法在广告支出数据上表现良好。但是用这些方法组合/合并新的信号，例如事件、天气，并不容易。这些统计方法对缺失数据也非常敏感，通常不能很好地预测长期数据。

在下一篇文章中，我将展示如何使用深度学习来预测同一数据集上的时间序列！

有用的参考资料:

[## 实验 10，第 1 部分-指数平滑

首先，我们区分方法和模型。预测方法是一种算法，它提供一个点…

www.stat.berkeley.edu](https://www.stat.berkeley.edu/~arturof/Teaching/STAT248/lab10_part1.html) [## 指数平滑法

指数平滑是一种使用指数窗口平滑时间序列数据的经验法则…

en.wikipedia.org](https://en.wikipedia.org/wiki/Exponential_smoothing) [## ARIMA 模型介绍

非季节性 ARIMA 模型说明(pdf 文件)季节性和非季节性 ARIMA 模型幻灯片(pdf 文件)介绍…

people.duke.edu](https://people.duke.edu/~rnau/411arim.htm)

时间序列土地覆盖挑战:深度学习视角

原文：https://towardsdatascience.com/time-series-land-cover-challenge-a-deep-learning-perspective-6a953368a2bd?source=collection_archive---------15-----------------------

针对 TiSeLaC 挑战的不同 DL 架构快速回顾

蒙彼利埃大学在 2017 年提出的 TiSeLaC 挑战[1】(TiSeLaC 针对时间序列土地覆盖)在于预测卫星图像时间序列中的土地覆盖类像素。

照片由 Kelly Lacy 从 Pexels

目录

1. 什么是时间序列卫星图像
2. TiSeLaC 数据集呢？
3. TiSeLaC 分类任务
4. 结论

1.什么是时间序列卫星图像？

时间序列卫星图像是对卫星图像时间维度的补充。换句话说，它是一系列定期拍摄的卫星图片，以便不仅使用图片中嵌入的空间信息，而且使用时间维度来进行预测，无论这些预测是分类、检测还是分割。

X. Yang 和 C. P. Lo [2]在 2002 年对卫星图像时间序列进行了历史上最引人注目的使用，他们揭示并量化了在亚特兰大加速城市发展的背景下森林的损失和城市扩张。

2.TiSeLaC 数据集呢？

这些卫星图片是 2014 年以 30 米分辨率拍摄的 23 张 2866 * 2633 像素的留尼汪岛图片。这些像素中的每一个都由 10 个特征组成:7 个表面反射率，代表每个独立多光谱带(OLI)的测量值:超蓝、蓝、绿、红、NIR、SWIR1 和 SWIR2。我们还发现了 3 个互补的辐射指数，分别是归一化差异植被指数、归一化差异水指数以及最后的亮度指数。

带有 10 维细节数据的图形像素(来源: 托马斯·迪·马蒂诺 )

有这么多波段，我们只能想象波段之间的潜在相关性。我们现在将深入探索性数据分析。

2.1 探索性数据分析

这里要注意的第一件重要的事情是，在每幅图像包含的 2866*2633 个像素中，只有 81714 个像素被保留在训练集中，测试集中有 17973 个像素。这意味着我们不处理整个图像，因为只有像素的子集被处理和分析。

以下动画显示了每个像素的蓝色、绿色和红色分量随时间的变化:

使用 RGB 波段的卫星数据的 23 天 gif 动画(来源: 托马斯·迪·马蒂诺 )

此外，为了对每个类栅格有一个总体了解，我通过计算 23 天中每一天的 10 个波段各自的平均值和标准差，为每个类创建了一个图。

“每一类”的时间序列意味着(来源: 托马斯·迪·马蒂诺 )

正如我们在那里看到的，一些不同的类别看到它们的特征随着时间的推移被高度扰乱(例如水类别或岩石和裸露土壤是最明显的)。相反，森林类和农作物类的结果相当稳定。

每个类别的特征变化中的这种明显多样性是一个“好迹象”，对于自动分类器来说，能够将它们彼此区分开应该不会太难。然而，这些结果充其量是一般化的，可能会产生误导，它们只是提供了关于数据多样性的见解，但并不能证明分类任务(或至少是高度精确的分类)会很容易。

我们现在将探索类的分布。TiSeLaC 组织者对像素进行的二次采样的第一个目的是平衡等级分布。然后，我们希望看到一个稍微平衡的数据集:

训练集的类的柱状图(来源: 托马斯·迪·马蒂诺 )

我们看到，数据集实际上是倾斜的，因为“其他作物”类的训练样本少于 2000 个，而城市地区、森林和稀疏植被都是向上采样的。

对于测试集，我们有:

训练集的类的条形图(来源: 托马斯·迪·马蒂诺 )

我们在这里有一个有点类似的图，每个类都比训练集中的少 4 到 6 倍，除了“*其他作物”*在测试集中更少。

2.2 分类问题的初步提示

由于这个数据集是由二次抽样数据组成的，在某种程度上，只有一些像素被保留，我们不能像图像处理一样处理这个问题，因为我们只能访问图像的一部分。

然而，我们可以从信号处理的角度来看问题，更具体地说，可以从时间序列分类的角度来看问题，由于像素坐标数据，时间序列分类具有关于所述时间序列在空间中的定位的额外信息。

3.TiSeLaC 分类任务

一旦我们采取了时间序列分类的方向，我们就可以比较最近流行的针对时间序列文学的深度学习的不同模型。

3.1 多重单峰网络:多通道深度卷积神经网络

最受欢迎的模型之一是在[3]中开发和研究的多通道深度卷积神经网络(即 MCDCNN)。这种架构希望通过在输入的每个维度上独立(即并行)应用卷积来利用多模态时间序列数据的不同特征之间的假定独立性。

我自己用 Python 和 Keras 实现了一个。

首先，我描述了 10 个通道各自的架构:

然后我用一个连接层将它们混合成一个模型。

对于训练过程，我使用了一个 SGD ，其学习率为 0.01 ，一个衰减为 0.0005，一个批量为 64 个实例，以及 120 个时期。

考虑到 TiSeLaC 组织者使用的 F1 评分标准,我得到了 0.867 的分

然而，对这个数据集的评估是，我们正在处理多模态时间序列。

3.2 多模式网络:时代-CNN

事实上，人们可能会认为指数的值与反射率测量值相关。这一点在“*岩石和裸土”*类地块中尤为明显，在第 11 天左右，多条带中的峰值非常明显。

这个想法在[4]中被赵 B、陆 H、陈 S、刘 J 和吴 D 讨论过，他们引入了时间的概念。该模型在多个方面不同于 MCDCNN:

• 它使用 sigmoid 输出层的 MSE 损失，而不是通常的分类交叉熵和 softmax 输出层；
• 它使用平均池代替通常的最大池；
• 在最后一个卷积层之后不存在池层

我在 Keras 和 Python 中实现的 Time-CNN 如下:

这一次，我使用了 1e-3 的学习率，没有衰减的 Adam 优化器，128 的批量和 100 个时期的(使用提前停止，通常在 50 个时期左右停止训练过程)。****

通过这种技术，我取得了 0.878 的 F1 分数，这是正确的，但仍然不令人满意。

正如我们所看到的，两个模型在考虑如何完成这项任务时采取了完全不同的角度:一个模型假设时间序列是不相关的，而另一个模型认为 10 个时间序列是一个独特的实体。这就是我认为使用这两种想法并将它们合并到一个模型中会有所帮助的地方。

3.3 多模式和多单模式架构的组合

这个想法在[5]中被 TiSeLaC 竞赛的一组研究人员深入探讨，他们的解决方案获得了冠军。我的实现深受他们工作的启发，虽然我没有达到与他们相同的性能(原因是他们使用像素的预处理空间特征表示以及通常的时间序列，并使用不同初始化的多个模型的打包技术完成)，但我仍然设法获得了令人满意的分数。

我的架构如下:

我提议的最终网络架构，由 3 种不同模型组合而成(来源: 托马斯·迪·马蒂诺 )

所提出的架构使用 3 种不同的模型，从左边的多变量模型**、中间的 10 个单变量模型和右边的用于位置信息的聚集模型开始。**

第一个模型是单变量模型，仅包含 3 个卷积层，中间没有池层。

第二个模型使用 10 个单变量模型，具有不同的连接级别，以处理在网络的不同级别设计的功能，有点类似于 UNet 或 ResNet 会做的事情。

然后，第三个模型仅通过预处理和缩放的像素坐标到达最终的完全连接的层。

然后，将这些特征提取模型输出中的每一个连接起来，用通常的全连接层进行分类:

有了这个架构，我获得了 0.930 的 F1 分数，0.001 的学习率，256** 的批量， 50 个历元和一个默认参数为的 Adam 优化器。**

TiSeLaC 排行榜(来源: TiSeLaC 网站 )

有了这个分数，我的解决方案将位于 GIT 团队之上。然而，它仍然需要通过在高性能计算机上运行可能的 GridSearch 算法以及为像素位置找到一个好的预处理思想来彻底调整超参数，而不是简单的缩放过程。

4.结论

通过这个项目，我们看到了如今的 DL 架构在对时间序列数据进行分类时是多么有效，以及如何将单峰和多峰分析结合起来以获得更好的性能。

作为一个更普遍的结论，我们已经看到深度学习在处理序列方面有多好。

最后一点，我想感谢 TiSeLaC 竞赛的组织者公开了这个数据集。通过这种方式，我已经能够处理它，并了解更多关于构建不需要数百万个参数的高效模型的信息。

整个代码和执行程序都保存在我的个人 GitHub 上的一个 Jupyter 笔记本里，在 这个链接 。去看看！

文献学

**[1] R. Gaetano，D. Ienco。 TiSeLaC:时间序列土地覆盖分类挑战数据集。法国蒙彼利埃的 UMR·泰蒂斯。2017

[2] X .杨，C. P. Lo。使用时间序列的卫星图像检测佐治亚州亚特兰大大都市地区的土地利用和土地覆盖变化。国际遥感杂志。第 1775-1798 页。2002.

**[3]郑，刘，陈，葛，赵。利用多通道深度卷积神经网络进行多元时间序列分类。计算机科学前沿。2014.

[4]赵 B，陆 H，陈 S，刘 J，吴 D. 卷积神经网络用于时间序列分类。系统工程与电子学杂志。第 162-169 页。2017 年 2 月。

[5]迪毛罗、韦尔加里、巴西尔、文托拉、埃斯波里图。用于卫星图像时序分类的深度时空表示的端到端学习。2017。

大气 CO2 浓度(ppm)的时间序列建模，1958–2019

原文：https://towardsdatascience.com/time-series-modeling-for-atmospheric-co2-concentration-ppm-1958-2019-8425fadcb927?source=collection_archive---------14-----------------------

R +中季节性 ARIMA 建模的分步方法(源代码)

在本帖中，我将介绍历史二氧化碳浓度(1958 年至 2019 年)的时间序列建模。季节性 ARIMA 建模将用于建立模型。Rmarkdown 格式的源代码可以在作者的 GitHub 链接这里找到。知道人们常说二氧化碳浓度应该在本世纪下半叶开始之前达到 460 ppm 以下，具体来说，这篇文章试图回答这个问题:

• 如果二氧化碳浓度的当前趋势继续下去(一切照旧)，那么在本世纪下半叶达到 460 ppm 的可能性有多大？

好，我们开始吧！

初始化

首先要做的是导入这个分析中需要的包(如果您还没有安装它们，请预先在 R 中安装这些包):

library(tidyverse)
library(forecast)
library(lubridate)
library(car)
library(scales)
library(patchwork)
library(kableExtra)

数据导入:

本分析将使用在夏威夷莫纳罗亚天文台测量的“1960/3–2019/12”期间的月 CO2 浓度(ppm)数据。此数据的链接可从以下网址获得:

【https://www.esrl.noaa.gov/gmd/ccgg/trends/data.html

现在，要读取 R 中的数据，必须考虑几点:

• 数据有注释，这不是我们分析的兴趣。因此，我们用*comment.char = ‘#’ *，
• 数据有 7 列，但其中一些列(如年和月标签)尚未写入数据，因此在导入时，我们分配以下列名(年、月、时间、Co2_con、插值、趋势、天)。
• 使用空格分隔列，因此 sep =“”将被添加到代码中。
• 因为我们包含了列名，所以可以包含 header = F。

data <- read.delim('ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_mm_mlo.txt', comment.char = '#', header = F, sep = '', col.names = c('Year','Month','Time','Co2_Concentration','Interpolated','Trend','Days'))

使数据整洁:

查看任何 NA 值:

which(is.na(data))## integer(0)

好，我们有完整的测量数据！但是，在读取数据时我们看到一些-99.99 的值！—注意，如注释中所述，这些值是在测量不可用时的值，因此对于这些点(741 次测量中的 7 次)，我们使用插值列:

data_cc <- data %>% 
    mutate(
        Co2_Con = case_when(
            Co2_Concentration == -99.99 ~ Interpolated,
            TRUE ~ Co2_Concentration
        )
    )

让我们看看列类型:

sapply(data_cc, class)##              Year             Month              Time Co2_Concentration 
##         "integer"         "integer"         "numeric"         "numeric" 
##      Interpolated             Trend              Days           Co2_Con 
##         "numeric"         "numeric"         "integer"         "numeric"

我们可以看到列类型采用了适当的格式，但是我们可以添加名为 Date 的新列，它以标准的时间序列格式给出测量日期。

数据转换

这里 Lubridate 包提供了一个简单的方法来转换我们的年和月列到日期:

data_cc$Date <- ymd(paste0(data$Year, " ", data$Month, " ", "15"))

此外，我们可以在我们想要做的分析中看到，我们不需要所有的列，因此我们可以选择我们分析中需要的列:

data_cc_sel <- data_cc %>% 
    select(Year, Month, Date, Co2_Con )

此外，我们需要有一部分数据，以测试我们基于训练数据开发的模型-因此，在这里，我们认为 2017 年、 2018 年和 2019 年的数据是测试数据，其余的是训练数据。

data_cc_sel_test <- data_cc_sel %>% 
    filter(Year > 2016)
data_cc_sel_train <- data_cc_sel %>% 
    filter(Year <= 2016)

数据可视化

现在，让我们先把数据可视化，

ggplot(data_cc_sel,aes(Date, Co2_Con)) +
    geom_line(color='blue') +
    xlab("Year, Month") +
    scale_x_date(date_labels = "%Y-%m", date_breaks = "5 year") +
    theme(axis.text.x = element_text(face = "bold", color = "#993333", 
                           size = 12, angle = 45, hjust = 1)) +
    ylab("CO2 Concentration (ppm)") +
    #scale_x_continuous(breaks = trans_breaks(identity, identity, n = 10))
    scale_y_continuous() +
    theme(axis.text.y = element_text(face = "bold", color = "#993333", 
                           size = 10, hjust = 1),axis.title.y = element_text(size = 10))

现在，有时使用“拼凑”包将所有总数据集、训练和测试一个接一个地绘制出来也很好:

p1 <- ggplot(data_cc_sel,aes(Date, Co2_Con)) +
    geom_line(color='blue') +
    xlab("Year, Month") +
    scale_x_date(date_labels = "%Y-%m", date_breaks = "5 year") +
    theme(axis.text.x = element_text(face = "bold", color = "#993333", 
                           size = 12, angle = 45, hjust = 1)) +
    ylab("CO2 Concentration (ppm)") +
    #scale_x_continuous(breaks = trans_breaks(identity, identity, n = 10))
    scale_y_continuous() +
    theme(axis.text.y = element_text(face = "bold", color = "#993333", 
                           size = 10, hjust = 1),axis.title.y = element_text(size = 8)) p2 <- ggplot(data_cc_sel_train,aes(Date, Co2_Con)) +
    geom_line(color='blue') +
    xlab("Year, Month") +
    scale_x_date(date_labels = "%Y-%m", date_breaks = "5 year") +
    theme(axis.text.x = element_text(face = "bold", color = "#993333", 
                           size = 12, angle = 45, hjust = 1)) +
    ylab("CO2 Concentration (ppm)") +
    #scale_x_continuous(breaks = trans_breaks(identity, identity, n = 10))
    scale_y_continuous() +
    theme(axis.text.y = element_text(face = "bold", color = "#993333", 
                           size = 10, hjust = 1), axis.title.y = element_text(size = 8)) p3 <- ggplot(data_cc_sel_test,aes(Date, Co2_Con)) +
    geom_line(color='blue') +
    xlab("Year, Month") +
    scale_x_date(date_labels = "%Y-%m", date_breaks = "1 year") +
    theme(axis.text.x = element_text(face = "bold", color = "#993333", 
                           size = 12, angle = 45, hjust = 1)) +
    ylab("CO2 Concentration (ppm)") +
    #scale_x_continuous(breaks = trans_breaks(identity, identity, n = 10))
    scale_y_continuous() +
    theme(axis.text.y = element_text(face = "bold", color = "#993333", 
                           size = 10, hjust = 1), axis.title.y = element_text(size = 8)) (p2 | p3 ) /
      p1

建模:

在时间序列分析中，关于趋势，我们需要知道的前三件事是:

• 数据是静态的吗？
• 回答:不是，我们在图中看到了清晰的趋势，因此 CO2 浓度的特性取决于时间(例如，数据的平均值取决于时间，是非平稳的)。
• 数据有季节性吗？
• 回答:是的，我们在数据中肯定能看到季节性。现在，知道了数据的非平稳性和季节性，它建议使用季节性差异来建模数据。要回答，
• 自相关函数和偏相关是怎样的？

这是来自预测包的 ACF 和 PACF 的图:

Co2_train <- ts(data_cc_sel_train$Co2_Con, start = c(1958,3), frequency = 12)
Co2_train %>% ggtsdisplay()

很明显，数据显示需要差分(因为 ACF 不会下降到零)，现在我们对进行普通差分，滞后时间为 12:

Co2_train %>% diff(lag=12) %>% diff() %>% ggtsdisplay()

现在，情况有所好转，我们基本上消除了这种趋势。我们以 ARIMA(p，D，q)(P，D，Q)[12]中的 d=D = 1 开始模型(原因是我们在上面的图中区分了季节性和非季节性部分。)

现在，我们必须有一些 P，Q，P，Q 的起始参数。所以，让我们看看上面的 ACF 和 PACF:

• 在季节性滞后中，ACF 中有一个显著的尖峰，表明可能存在 MA(1)项。所以，起点是 Q = 1
• 在非季节性差异数据的图中，ACF 图有三个尖峰，这可能暗示季节性 MA(3)项，q=3。

因此，我们开始使用 ARIMA(0，1，3)(3，1，1)[12]并在 AR 和 MA 项中进行变化。这里，在保持阶常数(D，D)的同时，我们使用 AICs 值来判断模型的质量。(尽量减少 AIC)

aicsvalue <- function(p,q,P,Q) {
  fit <- Arima(Co2_train, order=c(p,1,q),seasonal=list(order=c(P,1,Q),period=12),
                  lambda = "auto"
              )
  return(fit$aicc)
}model_eva <- data.frame(Model_name=c("ARIMA(0,1,3)(3,1,1)[12]","ARIMA(0,1,1)(3,1,1)[12]","ARIMA(1,1,0)(1,1,0)[12]",
                                     "ARIMA(1,1,2)(1,1,0)[12]","ARIMA(1,1,3)(0,1,1)[12]","ARIMA(1,1,1)(1,1,0)[12]",
                                     "ARIMA(1,1,1)(1,1,0)[12]","ARIMA(1,1,0)(1,1,1)[12]","ARIMA(1,1,1)(0,1,1)[12]" ), AICc=c(aicsvalue(0,3,3,1),aicsvalue(0,1,3,1),aicsvalue(1,0,1,0),                                                            aicsvalue(1,2,1,0),aicsvalue(1,3,0,1),aicsvalue(1,1,1,0),                                                           aicsvalue(1,1,1,0),aicsvalue(1,0,1,1), aicsvalue(1,1,0,1)))

基于以上分析，将选择ARIMA(1，1，1)(0，1，1)[12]，但是我们需要检查残差以避免任何过拟合和欠拟合，以查看 Ljung-Box 测试残差是否类似白噪声。

(fit_minaicc <- Arima(Co2_train, order=c(1,1,1),seasonal=list(order=c(0,1,1),period=12),
                  lambda = "auto"
              ))
checkresiduals(fit_minaicc, lag=36)
fit_minaicc$aicc

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,1)(1,1,1)[12] with drift
## Q* = 32.406, df = 30, p-value = 0.3489
## 
## Model df: 6\.   Total lags used: 36

现在，我们可以看到残差非常类似于白噪声，p 值也很高，模型通过了 Ljong-Box 测试。(然而，必须提到的是，一些 ACF 刚刚达到蓝线的边界，然而，我不认为它会对预测产生实质性影响——有时很难让模型通过所有测试。)

然而，这并不是模型选择的终点。这里比较一下模型在测试数据上的表现。我们寻求最小化 RMSE 的模型。

Co2_test <- ts(data_cc_sel_test$Co2_Con, start = c(2017,1), frequency = 12)
mm <- accuracy(forecast(fit_minaicc,h=35)$mean, Co2_test )

本节比较了上一节中提供的 9 种车型的 RMSE 值。

rmse_eva <- function(p,d,q,P,D,Q) {
  fit <- Arima(Co2_train, order=c(p,d,q),seasonal=list(order=c(P,D,Q),period=12),
                  lambda = "auto", include.drift = T
              )
  mm <- accuracy(forecast(fit,h=35)$mean, Co2_test)
  return(mm[2])}rmse_eva <- data.frame(Model_name=c(
"ARIMA(0,1,3)(3,1,1)[12]","ARIMA(0,1,1)(3,1,1)[12]","ARIMA(1,1,0)(1,1,0)[12]",
"ARIMA(1,1,2)(1,1,0)[12]","ARIMA(1,1,3)(0,1,1)[12]","ARIMA(1,1,1)(1,1,0)[12]",
"ARIMA(1,1,1)(1,1,0)[12]","ARIMA(1,1,0)(1,1,1)[12]","ARIMA(1,1,1)(0,1,1)[12]"
 ), RMSE=c(                        
rmse_eva(0,1,3,3,1,1),rmse_eva(0,1,1,3,1,1),rmse_eva(1,1,0,1,1,0),                                                  rmse_eva(1,1,2,1,1,0),rmse_eva(1,1,3,0,1,1),rmse_eva(1,1,1,1,1,0),                                                  rmse_eva(1,1,1,1,1,0),rmse_eva(1,1,0,1,1,1),rmse_eva(1,1,1,0,1,1)))

结果表明，模型 ARIMA(1，1，1)(0，1，1)[12]没有最小 RMSE 值，但它非常接近最小值，但它在 AICc 值中是最小的。最后，已知模型残差遵循白噪声，选择模型 ARIMA(1，1，1)(0，1，1)[12]来预测包装，因为它的参数较少且简单，同时保持最小的 AIC。(注意:如果标准被设置为 RMSE，这可能会导致我们的模型非常复杂，有许多参数)

Co2_train %>%
  Arima(order=c(1,1,1),seasonal=list(order=c(0,1,1),period=12),
                lambda = "auto"
              ) %>%
  forecast(h=400) %>%
  autoplot() +
    ylab("H02 sales (million scripts)") + xlab("Year") +
  autolayer(Co2_test)

让我们放大模型预测和测试数据，以直观地查看模型性能:

prediction <- forecast(fit_minaicc,h=35) 
data_cc_sel_test$prediction <- prediction$mean
data_test_pre_tidy <- gather(data_cc_sel_test, "type", "Co2", -Year,-Month,-Date)## Warning: attributes are not identical across measure variables;
## they will be droppedggplot(data_test_pre_tidy,aes(Date, Co2,color=type)) +
    geom_line() +
    xlab("Year, Month") +
    scale_x_date(date_labels = "%Y-%m", date_breaks = "1 year") +
    theme(axis.text.x = element_text(face = "bold", color = "#993333", 
                           size = 12, angle = 45, hjust = 1)) +
    ylab("CO2 Concentration (ppm)") +
    #scale_x_continuous(breaks = trans_breaks(identity, identity, n = 10))
    scale_y_continuous() +
    theme(axis.text.y = element_text(face = "bold", color = "#993333", 
                           size = 10, hjust = 1), axis.title.y = element_text(size = 8))

洞察力:

现在，给定开发的模型，我们想要回答的问题是:

给定开发的模型，2050 年达到 460 ppm 的几率有多大？要回答这个问题，我们首先需要建立 2050 年二氧化碳浓度的累积分布:

prediction1 <- forecast(fit_minaicc,h=396, level = c(80,90)) 
p10 <- prediction1$upper[396,2]
p50 <- prediction1$mean[396]
sd_calc <- (p10-p50)/1.28Co2_con_2050 <- rnorm(10^6,p50,sd_calc)
cdf_co2_con_2050 <- ecdf(Co2_con_2050)
cdf_co2_con_2050_data <- data.frame(Co2_con_2050)
ggplot(cdf_co2_con_2050_data, aes(Co2_con_2050)) + stat_ecdf(geom = "step", color='blue') +
  geom_vline(xintercept = 460, color='red') +
  geom_hline(yintercept = cdf_co2_con_2050(460), color='red') +
  theme(axis.text.x = element_text(face = "bold", color = "#993333", 
                           size = 12, angle = 0, hjust = 1)) +
  scale_x_continuous(breaks=c(400,425,450, 460,475,500,525, 550), limits = c(425,525)) +
  scale_y_continuous(breaks=c(seq(0,1,0.1)), limits = c(0,1)) +
  ylab('Cumulative Distribution') +
  xlab("Co2 Concentraion(ppm) at 2050")## Warning: Removed 238 rows containing non-finite values (stat_ecdf).

现在，有了累积分布，我们可以问这个问题:

• 到 2050 年，二氧化碳浓度(ppm)保持在 460 以下的概率有多大？

cdf_co2_con_2050(460)
## [1] 0.089823

可以看到，答案是 9%左右。

使用 Scikit、Pandas 和 Numpy 进行时间序列建模

原文：https://towardsdatascience.com/time-series-modeling-using-scikit-pandas-and-numpy-682e3b8db8d1?source=collection_archive---------0-----------------------

直观地利用季节性来提高模型准确性

图片由作者

欢迎学习时间序列分析的第 2 部分！在本帖中，我们将通过建模时间序列数据来学习我们的方法。这是我上一篇关于时间序列数据的文章的延续。

在我们之前的博客文章中，我们讨论了什么是时间序列数据，如何格式化这样的数据以最大化其效用，以及如何处理缺失数据。我们还学习了如何按月、周、年等对时间序列数据进行重采样，并计算滚动平均值。我们深入研究了趋势、季节性、一阶差分和自相关等概念。如果你熟悉大部分的东西，你就可以开始了。如果你需要复习，你可以在谷歌上快速搜索这些话题，或者在这里阅读我之前的文章。

在我们开始之前说几句话:

毫无疑问，还有其他更好的用于时间序列预测的软件包，比如 ARIMA 的或者脸书的专有软件 Prophet 。然而，这篇文章的灵感来自一个朋友的带回家的作业，该作业要求她只能使用 Scikit、Numpy 和 Pandas(否则将面临立即取消资格！).

让我们深入我们的数据集

我们将使用公开的数据集开放电力系统数据。你可以在这里下载数据。它包含 2006-2017 年的电力消耗、风力发电和太阳能发电。

url='[https://raw.githubusercontent.com/jenfly/opsd/master/opsd_germany_daily.csv'](https://raw.githubusercontent.com/jenfly/opsd/master/opsd_germany_daily.csv')
data = pd.read_csv(url,sep=",")

将Date列设置为索引后，我们的数据集看起来是这样的:

*# to explicitly convert the date column to type DATETIME*
data['Date'] = pd.to_datetime(data['Date'])
data = data.set_index('Date')

定义建模任务

预测的目标

我们的目标是从这个时间序列数据集中预测Consumption(理想情况下是未来未知的日期)。

训练和测试设备

我们将使用 10 年的数据进行培训，即 2006-2016 年，使用去年的数据进行测试，即 2017 年。

工作指标

为了评估我们的模型有多好，我们将使用 R-squared 和均方根误差(但将打印所有相关指标，以便您进行最终通话)。

助手功能

为了打印与回归任务相关的所有性能指标(比如 MAE 和 R-square)，我们将定义regression_results函数。

import sklearn.metrics as metrics
def regression_results(y_true, y_pred): *# Regression metrics*
    explained_variance=metrics.explained_variance_score(y_true, y_pred)
    mean_absolute_error=metrics.mean_absolute_error(y_true, y_pred) 
    mse=metrics.mean_squared_error(y_true, y_pred) 
    mean_squared_log_error=metrics.mean_squared_log_error(y_true, y_pred)
    median_absolute_error=metrics.median_absolute_error(y_true, y_pred)
    r2=metrics.r2_score(y_true, y_pred) print('explained_variance: ', round(explained_variance,4))    
    print('mean_squared_log_error: ', round(mean_squared_log_error,4))
    print('r2: ', round(r2,4))
    print('MAE: ', round(mean_absolute_error,4))
    print('MSE: ', round(mse,4))
    print('RMSE: ', round(np.sqrt(mse),4))

特征工程

作为基线，我们选择了一个简单的模型，该模型基于以下因素预测今天的消费价值

• 昨天的消费值和；
• 昨天和前天的消费值之差。

*# creating new dataframe from consumption column*
data_consumption = data[['Consumption']]*# inserting new column with yesterday's consumption values*
data_consumption.loc[:,'Yesterday'] = 
data_consumption.loc[:,'Consumption'].shift()*# inserting another column with difference between yesterday and day before yesterday's consumption values.* data_consumption.loc[:,'Yesterday_Diff'] = data_consumption.loc[:,'Yesterday'].diff()*# dropping NAs*
data_consumption = data_consumption.dropna()

定义训练集和测试集

X_train = data_consumption[:'2016'].drop(['Consumption'], axis = 1)
y_train = data_consumption.loc[:'2016', 'Consumption']X_test = data_consumption['2017'].drop(['Consumption'], axis = 1)
y_test = data_consumption.loc['2017', 'Consumption']

时间序列数据的交叉验证

在数据科学面试中经常出现的一个问题是:你会对时间序列数据使用哪种交叉验证技术？

你可能会被有史以来最受欢迎的 K 倍交叉验证所吸引(相信我，直到最近——不要问最近！—我不知道除了 K-fold 之外还有 CV 技术。不幸的是，这不是正确的答案。原因是，它没有考虑到时间序列数据有一些自然的顺序，标准 k 倍交叉验证中的随机化没有保持这种顺序。

对时间序列数据进行交叉验证的一个更好的替代方法(比 K-fold CV)是正向链接策略。

在正向链接中，假设有 3 个折叠，训练集和验证集看起来像:

• 折叠 1:培训[1]，验证[2]
• 折叠 2:培训[1 2]，验证[3]
• 折叠 3:培训[1 2 3]，验证[4]

其中 1，2，3，4 代表年份。这样，连续的训练集是它们之前的训练集的超集。

幸运的是，sklearn 有一个使用TimeSeriesSplit实现这种列车测试分割的规定。

from sklearn.model_selection import TimeSeriesSplit

TimeSerieSplit函数将分割数作为输入。由于我们的培训数据有 11 个独特的年份(2006 -2016)，我们将设置n_splits = 10。这样我们就有了整洁的训练和验证集:

• 折叠 1:培训[2006]，验证[2007]
• 折叠 2:培训[2006 年 2007 年]，验证[2008 年]
• 折叠 3:培训[2006 年 2007 年 2008 年]，验证[2009 年]
• 折叠 4:培训[2006 年 2007 年 2008 年 2009 年]，验证[2010 年]
• 折叠 5:培训[2006 年 2007 年 2008 年 2009 年 2010 年]，验证[2011 年]
• 折叠 6:培训[2006 年 2007 年 2008 年 2009 年 2010 年 2011 年]，验证[2012 年]
• 折叠 7:培训[2006 年 2007 年 2008 年 2009 年 2010 年 2011 年 2012 年]，验证[2013 年]
• 折叠 8:培训[2006 年 2007 年 2008 年 2009 年 2010 年 2011 年 2012 年 2013 年]，验证[2014 年]
• 折叠 9:培训[2006 年 2007 年 2008 年 2009 年 2010 年 2011 年 2012 年 2013 年 2014 年]，验证[2015 年]
• 折叠 10:培训[2006 年 2007 年 2008 年 2009 年 2010 年 2011 年 2012 年 2013 年 2014 年 2015 年]，验证[2016 年]

抽查算法

*# Spot Check Algorithms*models = []
models.append(('LR', LinearRegression()))
models.append(('NN', MLPRegressor(solver = 'lbfgs')))  #neural network
models.append(('KNN', KNeighborsRegressor())) 
models.append(('RF', RandomForestRegressor(n_estimators = 10))) # Ensemble method - collection of many decision trees
models.append(('SVR', SVR(gamma='auto'))) # kernel = linear*# Evaluate each model in turn*
results = []
names = []
for name, model in models:
    *# TimeSeries Cross validation*
 tscv = TimeSeriesSplit(n_splits=10)

 cv_results = cross_val_score(model, X_train, y_train, cv=tscv, scoring='r2')
 results.append(cv_results)
 names.append(name)
 print('%s: %f (%f)' % (name, cv_results.mean(), cv_results.std()))

*# Compare Algorithms*
plt.boxplot(results, labels=names)
plt.title('Algorithm Comparison')
plt.show()

KNN 和 RF 表现同样出色。但我个人更喜欢 RF，因为这种集合模型(将多个‘个体’(不同的)模型结合在一起，并提供卓越的预测能力。)几乎可以开箱即用，这也是它们非常受欢迎的原因之一。

网格搜索超参数

我在之前的文章中讨论了网格搜索超参数的必要性。

超参数的最佳组合使模型的性能最大化，而不会导致高方差问题(过度拟合)。

执行网格搜索的 Python 代码如下:

from sklearn.model_selection import GridSearchCVmodel = RandomForestRegressor()
param_search = { 
    'n_estimators': [20, 50, 100],
    'max_features': ['auto', 'sqrt', 'log2'],
    'max_depth' : [i for i in range(5,15)]
}tscv = TimeSeriesSplit(n_splits=10)
gsearch = GridSearchCV(estimator=model, cv=tscv, param_grid=param_search, scoring = rmse_score)gsearch.fit(X_train, y_train)
best_score = gsearch.best_score_
best_model = gsearch.best_estimator_

如果你注意到上面的代码，我们已经通过设置scoring = rmse_score定义了一个自定义 计分器，而不是使用 sklearn 中定义的通用计分指标之一。我们将自定义计分器定义如下:

from sklearn.metrics import make_scorerdef rmse(actual, predict):predict = np.array(predict)
    actual = np.array(actual)distance = predict - actualsquare_distance = distance ** 2mean_square_distance = square_distance.mean()score = np.sqrt(mean_square_distance)return scorermse_score = make_scorer(rmse, greater_is_better = False)

根据测试数据检查最佳模型性能

y_true = y_test.values
y_pred = best_model.predict(X_test)regression_results(y_true, y_pred)

这对初学者来说并不坏。让我们看看是否可以进一步改进我们的模型。

特征工程回报

到目前为止，我们一直使用(t-1)th日的值来预测t日的值。现在，让我们也使用(t-2)天的值来预测消耗量:

*# creating copy of original dataframe*
data_consumption_2o = data_consumption.copy()*# inserting column with yesterday-1 values*
data_consumption_2o['Yesterday-1'] = data_consumption_2o['Yesterday'].shift()*# inserting column with difference in yesterday-1 and yesterday-2 values.*
data_consumption_2o['Yesterday-1_Diff'] = data_consumption_2o['Yesterday-1'].diff()*# dropping NAs*
data_consumption_2o = data_consumption_2o.dropna()

重置列车和测试装置

X_train_2o = data_consumption_2o[:'2016'].drop(['Consumption'], axis = 1)
y_train_2o = data_consumption_2o.loc[:'2016', 'Consumption']X_test = data_consumption_2o['2017'].drop(['Consumption'], axis = 1)
y_test = data_consumption_2o.loc['2017', 'Consumption']

检查使用“新”预测器的“最佳”随机森林是否表现更好

model = RandomForestRegressor()
param_search = { 
    'n_estimators': [20, 50, 100],
    'max_features': ['auto', 'sqrt', 'log2'],
    'max_depth' : [i for i in range(5,15)]
}tscv = TimeSeriesSplit(n_splits=10)
gsearch = GridSearchCV(estimator=model, cv=tscv, param_grid=param_search, scoring = rmse_score)gsearch.fit(X_train_2o, y_train_2o)
best_score = gsearch.best_score_
best_model = gsearch.best_estimator_y_true = y_test.values
y_pred = best_model.predict(X_test)regression_results(y_true, y_pred)

好消息！！我们已经显著降低了 RMSE 和 MAE 值，而 R 平方值也上升了。

特征工程反击

让我们看看增加太阳能生产的价值是否在某种程度上有利于预测电力消耗。

data_consumption_2o_solar = data_consumption_2o.join(data[['Solar']])data_consumption_2o_solar = data_consumption_2o_solar.dropna()

重置训练/测试+网格搜索+检查性能

X_train_2o_solar = data_consumption_2o_solar[:'2016'].drop(['Consumption'], axis = 1)
y_train_2o_solar = data_consumption_2o_solar.loc[:'2016', 'Consumption']X_test = data_consumption_2o_solar['2017'].drop(['Consumption'], axis = 1)
y_test = data_consumption_2o_solar.loc['2017', 'Consumption']model = RandomForestRegressor()
param_search = { 
    'n_estimators': [20, 50, 100],
    'max_features': ['auto', 'sqrt', 'log2'],
    'max_depth' : [i for i in range(5,15)]
}tscv = TimeSeriesSplit(n_splits=5)
gsearch = GridSearchCV(estimator=model, cv=tscv, param_grid=param_search, scoring = rmse_score)gsearch.fit(X_train_2o_solar, y_train_2o_solar)
best_score = gsearch.best_score_
best_model = gsearch.best_estimator_y_true = y_test.values
y_pred = best_model.predict(X_test)regression_results(y_true, y_pred)

瞧，现在模型的性能更好了。

可变重要性图

imp = best_model.feature_importances_
features = X_train_2o_solar.columns
indices = np.argsort(imp)plt.title('Feature Importances')
plt.barh(range(len(indices)), imp[indices], color='b', align='center')
plt.yticks(range(len(indices)), [features[i] for i in indices])
plt.xlabel('Relative Importance')
plt.show()

正如我们所见，太阳能发电量并不像其他基于时间的预测指标那样是一个强大的电力消耗预测指标。

特征工程的最后阶段

如果你阅读了我前一篇博文第一部分的叙述，你会记得我们的数据集有一些季节性因素，更准确地说是每周季节性。因此，将给定日期前一周的消费值作为模型的输入更有意义。这意味着，如果模型试图预测 1 月 8 日的消费值，它必须获得 1 月 1 日的消费信息。

data_consumption_2o_solar_weeklyShift = data_consumption_2o_solar.copy()data_consumption_2o_solar_weeklyShift['Last_Week'] = data_consumption_2o_solar['Consumption'].shift(7)data_consumption_2o_solar_weeklyShift = data_consumption_2o_solar_weeklyShift.dropna()

重置训练/测试+网格搜索+检查性能

X_train_2o_solar_weeklyShift = data_consumption_2o_solar_weeklyShift[:'2016'].drop(['Consumption'], axis = 1)
y_train_2o_solar_weeklyShift = data_consumption_2o_solar_weeklyShift.loc[:'2016', 'Consumption']X_test = data_consumption_2o_solar_weeklyShift['2017'].drop(['Consumption'], axis = 1)
y_test = data_consumption_2o_solar_weeklyShift.loc['2017', 'Consumption']model = RandomForestRegressor()
param_search = { 
    'n_estimators': [20, 50, 100],
    'max_features': ['auto', 'sqrt', 'log2'],
    'max_depth' : [i for i in range(5,15)]
}tscv = TimeSeriesSplit(n_splits=10)
gsearch = GridSearchCV(estimator=model, cv=tscv, param_grid=param_search, scoring = rmse_score)gsearch.fit(X_train_2o_solar_weeklyShift, y_train_2o_solar_weeklyShift)
best_score = gsearch.best_score_
best_model = gsearch.best_estimator_y_true = y_test.values
y_pred = best_model.predict(X_test)regression_results(y_true, y_pred)

我们又做了一次…误差进一步减小，r 平方增加。我们可以继续添加更多相关的功能，但我想你现在已经明白了！

特征重要性图

正如我们正确假设的那样，第(t-7)天的价值比第(t-1)天的价值具有更强的预测力。

结论

在本文中，我们学习了如何对时间序列数据建模，对时间序列数据进行交叉验证，以及微调我们的模型超参数。我们还成功地将预测功耗的 RMSE 从 85.61 降低到 54.57。

在本系列的第 3 部分中，我们将进行一个案例研究，分析呼叫中心生成的时间序列数据，主要是分析放弃率的(可怕)增量。敬请关注…

直到下次:)

我喜欢写循序渐进的初学者指南、操作指南、面试问题、解码 ML/AI 中使用的术语等。如果你想完全访问我的所有文章(以及其他媒体上的文章)，那么你可以使用 我的链接这里 注册。

* [## 检测语音数据中的情感:使用 Huggingface 微调 HuBERT

构建自定义数据加载器、实验日志、改进指标的技巧和 GitHub repo，如果您想了解…

towardsdatascience.com](/fine-tuning-hubert-for-emotion-recognition-in-custom-audio-data-using-huggingface-c2d516b41cd8) [## 了解 Python 导入，init。py 和 pythonpath —一劳永逸

了解如何导入包和模块(以及两者之间的区别)

towardsdatascience.com](/understanding-python-imports-init-py-and-pythonpath-once-and-for-all-4c5249ab6355) [## 处理时序带回家的作业:Python 中的一个案例研究

利用呼叫中心分析改善客户支持

medium.com](https://medium.com/@vishi2020/tackling-the-time-series-take-home-assignment-a-case-study-in-python-b2a3bd78d956) [## 数据科学家的 Python 高效编码指南

我每天用来编写干净代码的技巧和窍门

towardsdatascience.com](/data-scientists-guide-to-efficient-coding-in-python-670c78a7bf79)*

用 ARIMA 时间序列建模预测未来房价

原文：https://towardsdatascience.com/time-series-modeling-with-arima-to-predict-future-house-price-9b180c3bbd2f?source=collection_archive---------12-----------------------

从 0 到 1 的时间序列分析和建模

艾萨克·史密斯在 Unsplash 上拍摄的照片

时间序列曾经是我在研究生院时试图避免的话题，因为我的同龄人分享的课程非常理论化，课堂上讨论的所有案例都与金融有关，不幸的是，我当时对金融不感兴趣。多亏了我的数据科学训练营，我又有了一次接触时间序列的机会，我发现它在许多情况下都非常实用和有用。这一次，我使用时间序列分析和模型来预测投资布鲁克林的 5 个最佳邮政编码，我和我丈夫打算在那里购买一套公寓。在这篇博文中，我将分享你需要知道的关于时间序列的基础知识，以及我是如何一步步用 ARIMA 模型预测房价的。

首先我要打下一个时间序列建模项目的架构。我将在后面的章节中详细解释。

第一步:数据处理

第二步:数据探索和可视化

第三步:确定模型方法和 KPI

步骤 4:在训练集上开发模型，并使用测试集进行验证

第五步:微调模型并进行预测

请记住，时间序列只是在一段时间内以一致的时间间隔测量的一系列定义明确的数据点。时间序列分析帮助我们理解隐藏的模式，有意义的特征和数据的统计。

第一步:数据处理——将原始数据转换成以时间为索引的时间序列

在处理来自 Zillow.com 的数据集时，我必须首先选择该城市，并使用 pd.melt 将数据框从宽格式重塑为长格式，然后转换为时间序列数据。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlinedfm.set_index('Month', inplace = True)

第二步:数据探索和可视化——季节分解、与 ACF 和 PCAF 的自相关等。

通过可视化，我们可以识别数据中潜在的趋势和故事。我们来看看布鲁克林的房价随着时间的变化。

从 1996 年到 2007 年可以看到总体上升趋势，随后从 2008 年到 2010 年年中出现波动。2011 年开始，房价更加稳定，再次持续上涨。我们可以非常自信地说，2008 年的房地产市场崩溃是这次波动的原因，我们确实希望跳过这段时间，以便对未来有一个更准确的预测。

2011 年后趋势更加稳定

我们关心的有三个重要特征:平稳性、季节性和自相关性。

平稳性

大多数时间序列模型都假设时间序列是平稳的，这意味着时间序列的统计特性如均值、方差等是平稳的。保持不变。理想情况下，我们希望有一个平稳的时间序列来建模。 Dickey-Fuller 检验可以用来检验一个时间序列是否平稳。注意零假设是:时间序列不是平稳的。

from statsmodels.tsa.stattools import adfullerdftest = adfuller(ts)

季节性

季节性是指在一段固定时间内重复的周期性变化和模式。有时，季节性可以与上升或下降趋势相结合。季节性分解总是有助于检测数据集中的季节性、趋势和任何噪声。以布鲁克林住房数据为例:

from statsmodels.tsa.seasonal import seasonal_decomposedecomposition = sm.tsa.seasonal_decompose(month_avg, model='additive')trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid*# Plot gathered statistics*
plt.figure(figsize=(12,8))
plt.subplot(411)
plt.plot(month_avg, label='Original', color='blue')
plt.legend(loc='best')
plt.subplot(412)
plt.plot(trend, label='Trend', color='blue')
plt.legend(loc='best')
plt.subplot(413)
plt.plot(seasonal,label='Seasonality', color='blue')
plt.legend(loc='best')
plt.subplot(414)
plt.plot(residual, label='Residuals', color='blue')
plt.legend(loc='best')
plt.tight_layout()

布鲁克林月平均房价的季节分解

从上面的图表中我们可以清楚地看到，随着每年的季节性，有一个上升的趋势。下一步是用 Dickey-Fuller 测试来测试我们的残差。

自相关

自相关有助于我们研究每个时间序列观测值如何与其最近(或不太最近)的过去相关联。很容易得出，明天的房价极有可能和今天的房价有关。 **ACF(自相关函数)**和 **PACF(偏自相关函数)**是两个强有力的工具。ACF 将时间序列的自相关表示为时滞的函数。PACF 试图移除观测值和先前时间步长的观测值的自相关中存在的间接相关性

检查所有邮政编码的月平均价格的 ACF 和 PACF。

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from matplotlib.pylab import rcParamsrcParams['figure.figsize']=7,5
plot_acf(month_avg); plt.xlim(0,24); plt.show()
plot_pacf(month_avg); plt.xlim(0,24); plt.ylim(-1,1);plt.show()

ACF 显示该时间序列与之前的时间周期具有自相关性，然而，PACF 没有显示出显著的偏相关。

如果我们从当前月份的值中减去 3 个月前的值，换句话说，这是一个滞后 3 的差值。我们可以在 PACF 图中看到，当滞后=2 时，存在负的偏自相关，这意味着滞后-1 差异在时间序列数据中是显著的。

plot_acf(month_avg.diff(periods=3).bfill()); plt.xlim(0,24); plt.show()
plot_pacf(month_avg.diff(periods=3).bfill()); plt.xlim(0,24); plt.ylim(-1,1);plt.show()

步骤 3:决定模型方法和 KPI

由于我们的数据集不是静态的，并且存在季节性成分，因此使用 SARIMA 模型——季节性 ARIMA(带有外生回归量的季节性自回归综合移动平均)是合理的。在不深入研究方法论的情况下，我现在将重点关注重要的参数。

根据公式 SARIMA( p ， d ， q )x( P ， D ， Q，s )，这些型号的参数如下:

• p 和季节性 P :表示自回归项(平稳序列的滞后)的个数
• d 和季节性 D :表示为了平稳化系列必须进行的差分
• q 和季节性 Q :表示移动平均项的个数(预测误差的滞后)
• s :表示时间序列的周期性(4 表示季度，12 表示年度)

KPI:使用 AIC 选择最佳参数集

第四步:开发模型并用测试集验证模型

由于布鲁克林的 Zillow 数据集中有 29 个邮政编码，我决定首先在 3 个样本邮政编码上构建 SARIMA 模型，然后遍历所有其他邮政编码。

我首先创建一个数据帧列表，每个数据帧都有一个邮政编码的信息。

zip_dfs = []
zip_list = dfm.RegionName.unique()
for x in zip_list:
    zip_dfs.append(pd.DataFrame(dfm[dfm['RegionName']==x][['MeanValue']].copy()))

然后定义 P，D，Q 和 P，D，Q，s 取 0 到 2 之间的任意值

p = d = q = range(0,2)
# Generate all different combinations of p, d and q triplets
pdq = list(itertools.product(p,d,q))
# Generate all different combinations of seasonal p, d and q triplets
pdqs = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]

萨里玛模型

ans = []for df, name in zip(zip_dfs, zip_list):
    for para1 in pdq:
        for para2 in pdqs:
            try:
                mod = sm.tsa.statespace.SARIMAX(df,
                                                order = para1,
                                                seasonal_order = para2,
                                                enforce_stationarity = False,
                                                enforce_invertibility = False)
                output = mod.fit()
                ans.append([name, para1, para2, output.aic])
                print('Result for {}'.format(name) + ' ARIMA {} x {}12 : AIC Calculated = {}'.format(para1, para2, output.aic))
            except:
                continue

然后将所有结果存储到数据框中

result = pd.DataFrame(ans, columns = ['name','pdq','pdqs','AIC'])

按最低 AIC 排序，以找到每个邮政编码的最佳参数

best_para = result.loc[result.groupby("name")["AIC"].idxmin()]

对部分数据进行动态预测，并与真实值进行比较

#Make Prediction and compare with real values
summary_table = pd.DataFrame()
Zipcode = []
MSE_Value = []
models = []
for name, pdq, pdqs, df in zip(best_para['name'], best_para['pdq'], best_para['pdqs'], zip_dfs):ARIMA_MODEL = sm.tsa.SARIMAX(df,
                                 order = pdq,
                                 seasonal_order = pdqs,
                                 enforce_stationarity = False,
                                 enforce_invertibility = False,
                                 )

    output = ARIMA_MODEL.fit()
    models.append(output)

    #get dynamic predictions starting 2017-06-01
    pred_dynamic = output.get_prediction(start=pd.to_datetime('2017-06-01'), dynamic = True, full_results = True)
    pred_dynamic_conf = pred_dynamic.conf_int()
    zip_forecasted = pred_dynamic.predicted_mean
    zip_truth = df['2017-06-01':]['MeanValue']

    sqrt_mse = np.sqrt(((zip_forecasted - zip_truth)**2).mean())

    Zipcode.append(name)
    MSE_Value.append(sqrt_mse)

summary_table['Zipcode'] = Zipcode
summary_table['Sqrt_MSE'] = MSE_Value

第五步:预测未来

下一步将是使用完整的数据集来预测未来的值。我用 3 年作为例子。

#Final Model 
forecast_table = pd.DataFrame()
current = []
forecast_3Yr = []

for zipcode, output, df in zip(Zipcode, models, zip_dfs):
    pred_3 = output.get_forecast(steps = 36)
    pred_conf_3 = pred_3.conf_int()
    forecast_3 = pred_3.predicted_mean.to_numpy()[-1]
    current.append(df['2018-04']['MeanValue'][0])
    forecast_3Yr.append(forecast_3) forecast_table['Zipcode'] = Zipcode
forecast_table['Current Value'] = current
forecast_table['3 Years Value'] = forecast_3Yr forecast_table['3Yr-ROI']=(forecast_table['3 Years Value'] - forecast_table['Current Value'])/forecast_table['Current Value']

这是我的最终结果:

• 11220 年:南日落公园(3 年投资回报率:17%-87%)
• 第 11205 名:克林顿·希尔(3 年投资回报率:16%-78%)
• 11203 年:东弗拉特布什(3 年投资回报率:8%-78%)
• 11224 年:科尼岛(3 年投资回报率:-0.5%-76%)
• 11217 年:博勒姆山(3 年投资回报率:6%-61%)

下一步

该模型纯粹基于时间序列，因此预测可能不太准确，因为显然还有许多其他因素影响房价，如经济、利率、房屋市场安全得分等。

如果要考虑其他因素，线性模型会更理想。

感谢阅读！让我知道你的想法。

有用的资源:

利用 Python 和自动 ARIMA 预测季节性时间序列

关于时间序列你需要知道的一切

用脸书预言家进行时间序列建模

原文：https://towardsdatascience.com/time-series-modeling-with-facebook-prophet-57f146a11d0d?source=collection_archive---------10-----------------------

分析和可视化时间序列数据的快速入门方法

当试图理解时间序列时，有太多的东西需要考虑。数据的总体趋势是什么？受季节性影响吗？我应该使用什么样的模型，它的性能如何？

所有这些问题都会让时间序列建模变得有点吓人，但也没那么糟糕。最近在为我的数据科学训练营做一个项目时，我尝试了脸书预言家，一个由……你知道，脸书开发的用于时间序列建模的开源包。我发现使用我的数据运行它非常快速和容易，所以在这篇文章中，我将向您展示如何做到这一点，并分享一个我编写的使用 Prophet all-in-one 进行建模和验证的函数。请继续阅读！

时间序列建模:基础

简而言之，时间序列建模就是寻找因变量和时间之间的关系。许多现实世界的数据可以以时间序列的形式出现:股票价格、温度、二氧化碳水平、我在万圣节每小时吃的零食数量等等。如果您想要预测时间序列中接下来会发生什么，那么您需要识别并消除时间对数据的影响，以便您可以将预测精力集中在任何潜在的、非基于时间的模式上。

实际上，这意味着您应该:

• 识别数据随时间变化的总体趋势。价值观总体上是上升还是下降？还是他们徘徊在一个中间点？
• 识别数据中任何重复的、季节性的模式。价值观是否每年夏天都会飙升，每年冬天都会跌至低谷？值在工作日还是周末更高？“季节性”趋势可能发生在任何时间尺度上(年、月、日、小时等)。)，只要它在你的总时间段内重复。
• 测试去除趋势和季节性后剩余的数据是否是平稳的。如果你的时间序列是平稳的，它的均值、方差和协方差从一项到下一项将是常数，而不是随时间变化。大多数时间序列模型都假设您提供给它们的数据是稳定的，因此如果您希望模型按预期工作，这是一个重要的步骤。

接下来往往是某种形式的 ARIMA (自回归综合移动平均线)建模。ARIMA 类似于线性回归，但是我们的解释性特征是我们比较每个数据点的过去时间点。

要在 Statsmodels 中运行一个 ARIMA 模型，你需要计算出有多少以及哪些过去的时间点用于比较，这可能会变得棘手且耗时。要找到最佳组合，您可能需要多次运行您的模型，每次都使用不同的术语组合(使用网格搜索或随机搜索)，然后选择产生最佳结果的一个。当您考虑更多选项时，这样做所需的时间会成倍增加！

脸书先知

脸书预言家在幕后自动化了其中的一些过程，所以你需要做的预处理和参数选择更少(你出错的机会也更少)。你可以想象，脸书对他们用户活动中基于时间的模式非常感兴趣，所以他们建立了 Prophet，使其在处理多种规模的季节性方面特别强大。Prophet 还被设计为可供开发人员和非开发人员使用，这意味着它易于编写，并且可以为您提供有吸引力的图表和方便的预测数据表，而您只需付出很少的努力。Prophet 真正希望从您那里得到的是特定格式的数据，我将在下一节介绍这一点。

为先知做准备

请记住，我在 Python 中使用 Prophet，所以如果您在 r 中使用它，情况可能会有所不同。当 Prophet 读取您的数据帧时，它会查找两个特定的列。一个必须称为“ds”并包含您的 datetime 对象，另一个必须称为“y”并包含您的因变量。确保您的数据是“长”格式，即每行代表一个时间点。

为了更容易理解，我将使用 Statsmodels 的数据集，记录 1958 年至 2001 年莫纳罗亚天文台的二氧化碳水平。正如你将看到的，这个数据既有上升趋势，又有很强的季节性模式。

# Load needed packages
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from fbprophet import Prophet as proph

# Load the "co2" dataset from sm.datasets
data_set = sm.datasets.co2.load()

# load in the data_set into pandas data_frame
CO2 = pd.DataFrame(data=data_set["data"])
CO2.rename(columns={"index": "date"}, inplace=True)

# set index to date column
CO2.set_index('date', inplace=True)

CO2.head()

现在，日期是这个数据框架的索引。

请注意，我并没有马上将我的数据帧放入 Prophet 的正确格式中；我想先用它做些事情。

二氧化碳水平每周记录一次，但假设我对月间隔更感兴趣。我将重新取样为每月一次:

# Resample from weekly to monthly
co2_month = CO2.resample('MS').mean()

# Backfill any missing values
co2_month.fillna(method='bfill', inplace=True)

co2_month.head()

现在数据代表月平均二氧化碳水平。

看看趋势和季节性:

# Plot co2_month
plt.plot(co2_month.index, co2_month.co2)

一路曲折到达顶端！

好吧，这是一个可怕的二氧化碳水平上升，但让我们继续。在使用 Prophet 之前，我们需要创建“ds”和“y”列。

# Reset the index to make the dates an ordinary column
co2_month.reset_index(inplace=True)

# Rename the columns
co2_month.columns = ['ds', 'y']
co2_month.tail()

我们有 526 行月度二氧化碳数据，为 Prophet 做了适当的标记。

现在我们准备好迎接一些预言了！

包装该函数

我想预测二氧化碳将如何随时间增长，但我也想知道我的模型表现如何。我可以用我所拥有的数据来做这件事，方法是分离出数据的一部分，预测这部分的值，然后将预测值与实际值进行比较。在这种情况下，我将去掉最后 10%的数据(53 个月)用于验证:

# Split the data 90:10
co2_train = co2_month[:473]
co2_test = co2_month[473:]

我为自己编写了一个函数，它可以接受我的训练数据、验证数据和整个数据集，并返回以下内容:

• 我的训练数据和预测的图表；
• 将我的预测与来自验证数据集的实际值进行比较的图；
• 显示基于整个数据集的长期预测的图；
• 显示整个数据集中总体趋势和季节性的图。

我的函数打印这些图，并返回包含要与验证数据和长期预测进行比较的预测的数据帧。如果您想要分析模型预测的精确值，这些会很有用。

这是我的函数:

# Define a function to do all the great stuff described above
def proph_it(train, test, whole, interval=0.95, forecast_periods1, 
             forecast_periods2):
    '''Uses Facebook Prophet to fit model to train set, evaluate  
    performance with test set, and forecast with whole dataset. The 
    model has a 95% confidence interval by default.

       Remember: datasets need to have two columns, `ds` and `y`.
       Dependencies: fbprophet, matplotlib.pyplot as plt
       Parameters:
          train: training data
          test: testing/validation data
          whole: all available data for forecasting
          interval: confidence interval (percent)
          forecast_periods1: number of months for forecast on 
              training data
          forecast_periods2: number of months for forecast on whole 
              dataset'''

    # Fit model to training data and forecast
    model = proph(interval_width=interval)
    model.fit(train)
    future = model.make_future_dataframe(periods=forecast_periods1, freq='MS')
    forecast = model.predict(future)

    # Plot the model and forecast
    model.plot(forecast, uncertainty=True)
    plt.title('Training data with forecast')
    plt.legend();

    # Make predictions and compare to test data
    y_pred = model.predict(test)

    # Plot the model, forecast, and actual (test) data
    model.plot(y_pred, uncertainty=True)
    plt.plot(test['ds'], test['y'], color='r', label='actual')
    plt.title('Validation data v. forecast')
    plt.legend();

    # Fit a new model to the whole dataset and forecast
    model2 = proph(interval_width=0.95)
    model2.fit(whole)
    future2 = model2.make_future_dataframe(periods=forecast_periods2, 
                                          freq='MS')
    forecast2 = model2.predict(future2)

    # Plot the model and forecast
    model2.plot(forecast2, uncertainty=True)
    plt.title('{}-month forecast'.format(forecast_periods2))
    plt.legend();

    # Plot the model components
    model2.plot_components(forecast);

    return y_pred, forecast2

请注意，我的函数被设置为按月处理数据，图中将有反映这一点的标题。您可以重构函数来处理不同规模的数据，但这符合我的目的。

快跑，先知，快跑！

让我们在三个数据集(训练、验证和整体)上运行这个函数。因为我任意保留了最后 10% (53 个月)的数据进行验证，所以我也任意预测超过 10 倍的数据，即 530 个月，也就是 44 年多一点。

# Run the wrapper function, supplying numbers of months for forecasts
short_term_pred, long_term_pred = proph_it(co2_train, co2_test, co2_month, 53, 530)

该功能现在将开始打印出图。这是前 90%的数据集以及对后 10%的预测:

现在我们可以将预测值与过去 53 个月的实际值进行比较:

现在利用所有可用数据进行长期预测。自然，95%的置信区间会随着时间的推移而急剧扩大。二氧化碳水平可能会降低……但也可能会升高。

这些是 Prophet 的“组件图”，它分解了数据中的趋势和季节性。看起来二氧化碳水平往往在春季达到年度最高，在秋季达到最低。

我现在可以直观地分析这些漂亮的图，如果我想看预测值本身，我手头也有。

这就是全部了！我在一个关于预测未来几年俄勒冈州房价上涨的项目中使用了这个函数。你可以在 GitHub 上查看代码。

感谢阅读！

跨贴自【jrkreiger.net】。

时间序列模型

原文：https://towardsdatascience.com/time-series-models-d9266f8ac7b0?source=collection_archive---------1-----------------------

AR，MA，ARMA，ARIMA

图片由皮克斯拜的 Gerd Altmann 提供

AR、MA、ARMA 和 ARIMA 模型用于根据为同一观测记录的先前时间点的历史数据来预测(t+1)的观测。但是，有必要确保时间序列在观测超时期间的历史数据上是平稳的。如果时间序列不是平稳的，那么我们可以在记录上应用差异因子，看看时间序列的图形是否是平稳的超时周期。

ACF(自相关函数)

自相关函数考虑了所有过去的观察值，而不考虑其对未来或当前时间段的影响。它计算 t 和(t-k)时间段之间的相关性。它包括 t 和(t-k)时间段之间的所有滞后或间隔。相关性总是使用皮尔逊相关公式来计算。

PACF(部分相关函数)

PACF 确定时间段 t 和 t-k 之间的部分相关性。它没有考虑 t 和 t-k 之间的所有时滞。例如，让我们假设今天的股价可能取决于前 3 天的股价，但它可能没有考虑昨天的股价收盘。因此，我们通过忽略两个时隙 t 和 t-k 之间的不重要的时滞，只考虑对未来时间周期有直接影响的时滞

如何区分何时使用 ACF 和 PACF？

让我们以一个村庄一年的糖果销售和收入为例。假设村里每 2 个月有一个节日，我们取出 12 个月的糖果销售和收入的历史数据。如果我们把时间画成月份，那么我们可以观察到，在计算糖果销售额时，我们只对隔月感兴趣，因为糖果销售额每两个月增加一次。但如果我们要考虑下个月的收入，那么我们就必须考虑去年全年的收入。

因此，在上述情况下，我们将使用 ACF 来找出未来产生的收入，但我们将使用 PACF 来找出下个月销售的糖果。

AR(自回归)模型

作者图片

在 t 的时间周期受到在不同时隙 t-1，t-2，t-3，……先前时间点的影响由该特定时间段的系数因子决定。任何特定公司 X 的股票价格可能取决于时间序列中所有以前的股票价格。这种模型计算过去时间序列的回归，并计算序列中的当前或未来值，称为自回归(AR)模型。

yt = β₁y-₁+β₂yₜ-₂+β₃yₜ-₃+…………+βₖyₜ-ₖ

考虑一个牛奶分销公司的例子，该公司每月在该国生产牛奶。考虑到去年生产的牛奶，我们想计算本月生产的牛奶量。我们首先计算当月所有 12 个滞后期的 PACF 值。如果任何特定月份的 PACF 值大于某个有效值，则只有这些值会被考虑用于模型分析。

例如，在上图中，值 1，2，3 到 12 显示了在给定滞后 t 的情况下，当月牛奶产量的直接影响(PACF)。如果我们考虑两个高于阈值的显著值，则该模型将被称为 AR(2)。

移动平均线模型

作者图片

t 处的时间周期受到在各个时隙 t-1、t-2、t-3、…的意外外部因素的影响…这些意想不到的影响被称为误差或残差。先前时间点的影响由该特定时间段的系数因子α决定。任何特定 X 公司的股票价格可能取决于一夜之间发生的公司合并，或者该公司因破产而关闭。这种模型计算过去时间序列的残差或误差，并计算序列的当前或未来值，称为移动平均(MA)模型。

yt = α₁ɛₜ-₁+α₂ɛₜ-₂+α₃ɛₜ-₃+…………+αₖɛₜ-ₖ

考虑一个在我生日时分发蛋糕的例子。假设你妈妈让你带糕点去派对。每年你都会错过判断参加聚会的邀请数量，并根据要求结束或多或少的蛋糕数量。实际和预期的差异导致了误差。因此，您希望避免今年的错误，因此我们对时间序列应用移动平均模型，并根据过去的总体错误计算今年所需的糕点数量。接下来，计算时间序列中所有滞后的 ACF 值。如果任何特定月份的 ACF 值大于一个重要值，则只有这些值将被考虑用于模型分析。

例如，在上图中，值 1、2、3 到 12 显示了当前月份糕点计数的总误差(ACF ),考虑到时间 t 和当前月份之间的所有滞后，给出了滞后 t。如果我们考虑高于阈值的两个重要值，那么该模型将被称为 MA(2)。

ARMA(自回归移动平均)模型

作者图片

这是一个由 AR 和 MA 模型组合而成的模型。在该模型中，预测时间序列的未来值时考虑了先前滞后和残差的影响。这里，β表示 AR 模型的系数，α表示 MA 模型的系数。

yt = β₁yₜ-₁+α₁ɛₜ-₁+β₂yₜ-₂+α₂ɛₜ-₂+β₃yₜ-₃+α₃﹍-̿+……+β﹍ y﹍+α﹍*﹍-̿*

考虑上面的图表，其中 MA 和 AR 值标有各自的有效值。让我们假设我们只考虑来自 AR 模型的 1 个重要值和来自 MA 模型的 1 个重要值。因此，ARMA 模型将从其他两个模型的组合值中获得，其数量级为 ARMA(1，1)。

ARIMA(自回归综合移动平均)模型

作者图片

我们知道，为了应用各种模型，我们必须首先将序列转换成平稳时间序列。为了达到同样的目的，我们应用差分或积分方法，从时间序列的 t 值中减去 t-1 值。在应用第一差分之后，如果我们仍然不能获得平稳的时间序列，那么我们再次应用二阶差分。

ARIMA 模型与 ARMA 模型非常相似，除了它包括一个称为积分(I)的因子，即在 ARIMA 模型中代表 I 的差分。因此，简而言之，ARIMA 模型是已经应用于模型的许多差异的组合，以使其稳定，先前的滞后数以及残差，以预测未来值。

考虑上面的图表，其中 MA 和 AR 值标有各自的有效值。让我们假设我们只考虑来自 AR 模型的 1 个重要值和来自 MA 模型的 1 个重要值。此外，该图最初是不稳定的，我们必须执行一次差分操作，以便转换成稳定的集合。因此，将从其他两个模型的组合值以及积分算子获得的 ARIMA 模型可以显示为 ARIMA(1，1，1)。

结论:

所有这些模型都给了我们对任何特定时间序列的洞察力或至少足够接近的预测。此外，哪种型号完全满足用户的需求也取决于用户。如果与其他模型相比，任何一个模型的错误率都较小，那么我们最好选择给出最接近估计值的模型。

希望这篇文章能帮助你更好地理解事情！！

基于 LSTM 的时间序列价格异常检测

原文：https://towardsdatascience.com/time-series-of-price-anomaly-detection-with-lstm-11a12ba4f6d9?source=collection_archive---------1-----------------------

强生、JNJ、Keras、Autoencoder、Tensorflow

自动编码器是一种无监督学习技术，尽管它们是使用监督学习方法训练的。目标是最小化基于损失函数的重建误差，例如均方误差。

在本帖中，我们将尝试使用 LSTM 自动编码器来检测约翰逊&约翰逊的历史股价时间序列数据中的异常。

数据可以从雅虎财经下载。我选择的时间段是从 1985 年 9 月 4 日到 2020 年 9 月 3 日。

使用 LSTM 自动编码器检测强生股票价格数据异常时，我们将遵循的步骤:

1. 用强生公司从 1985 年 9 月 4 日到 2013 年 9 月 3 日的股票价格数据训练一个 LSTM 自动编码器。我们假设没有异常，它们是正常的。
2. 使用 LSTM 自动编码器重建 2013 年 9 月 4 日至 2020 年 9 月 3 日测试数据的误差。
3. 如果测试数据的重建误差高于阈值，我们将该数据点标记为异常。

我们将分解一个 LSTM 自动编码器网络来逐层理解它们。

数据

LSTM _ 自动编码器 _ 异常. py

可视化时间序列

viz_timeseries.py

图 1

预处理

• 列车测试分离

火车 _ 测试. py

• 使数据标准化

标准化. py

• 创建序列

结合TIME_STEPS将输入数据转换成三维数组。根据 LSTM 网络的要求，阵列的形状应为[samples, TIME_STEPS, features]。

我们希望我们的网络有 30 天的记忆，所以我们设置了TIME_STEPS=30。

create _ sequences.py

建立模型

• 我们定义了重建 LSTM 自动编码器架构，其期望具有 30 个时间步长和一个特征的输入序列，并且输出具有 30 个时间步长和一个特征的序列。
• RepeatVector()重复输入 30 次。
• 设置return_sequences=True，那么输出仍然是一个序列。
• TimeDistributed(Dense(X_train.shape[2]))加在最后得到输出，其中X_train.shape[2]是输入数据中的特征数。

LSTM _ 自动编码器. py

图 2

训练模型

LSTM _ 自动编码器 _ 火车 _ 模型. py

图 3

plt.plot(history.history['loss'], label='Training loss')
plt.plot(history.history['val_loss'], label='Validation loss')
plt.legend();

图 4

model.evaluate(X_test, y_test)

确定异常

• 在训练数据上找到 MAE loss。
• 将训练数据中的最大 MAE 损失值作为reconstruction error threshold。
• 如果测试集中一个数据点的重建损失大于这个reconstruction error threshold值，那么我们将把这个数据点标记为异常。

LSTM _ 火车 _ 损失. py

图 5

test _ 测试 _ 损失. py

图 6

测试 _ 损耗 _ vs _ 阈值. py

图 7

anomalies = test_score_df.loc[test_score_df['anomaly'] == True]
anomalies.shape

如您所见，测试集中有 22 个数据点超过了reconstruction error threshold。

可视化异常

plot_anomalies.py

图 8

模型发现 3 月份出现一些低价异常，4 月份出现高价异常。据详细记载， JNJ 股票在 3 月份触及 2020 年低点，但由于对其冠状病毒疫苗的乐观预期，不到一个月之后迅速回升至高点。

Jupyter 笔记本可以在 Github 上找到。祝你一周愉快！

参考资料:

[## LSTM 自动编码器简介-机器学习掌握

LSTM 自动编码器是一个使用编码器-解码器 LSTM 的序列数据自动编码器的实现

machinelearningmastery.com](https://machinelearningmastery.com/lstm-autoencoders/) [## 基于 Keras 的时间序列数据异常检测

由 Coursera 项目网提供。在这个动手介绍异常检测的时间序列数据与 Keras…

www.coursera.org](https://www.coursera.org/projects/anomaly-detection-time-series-keras) [## Keras 文档:使用自动编码器的时间序列异常检测

作者:pavithrasv 创建日期:2020/05/31 最近修改时间:2020/05/31 描述:检测时间序列中的异常…

keras.io](https://keras.io/examples/timeseries/timeseries_anomaly_detection/)

基于空气质量传感器数据的时间序列模式识别

原文：https://towardsdatascience.com/time-series-pattern-recognition-with-air-quality-sensor-data-4b94710bb290?source=collection_archive---------5-----------------------

一个具有真实传感器数据的面向真实客户的项目

Marcin Jozwiak 在 Unsplash 上拍摄的照片

1。简介
2。探索性数据分析
∘ 2.1 模式变化
∘ 2.2 特征之间的相关性
3 .异常检测与模式识别
∘ 3.1 点异常检测(系统故障)
∘ 3.2 集体异常检测(外部事件)
∘ 3.3 聚类与模式识别(外部事件)
4 .结论
参考文献
附录—所选特征的散布矩阵
关于我

注:详细的项目报告和本帖使用的数据集可以在我的 GitHub 页面找到。

1.介绍

这个项目是我为客户做的自由数据科学工作的一部分。不需要保密协议，该项目不包含任何敏感信息。因此，我决定展示该项目的数据分析和建模部分，作为我个人数据科学作品集的一部分。客户信息已被匿名化。

在该项目中，提供了两个数据集，每个数据集由一周的空气质量传感器读数组成。它们用于完成以下四项任务:

1.发现数据集中的异常以自动标记事件

2.将异常归类为“系统故障”或“外部事件”

3.从数据集中的模式中提供任何其他有用的结论

4.可视化数据集中要素的相互依赖关系

在本报告中，我将简要介绍我用于数据分析、特征关联可视化、自动标记“系统故障”和“外部事件”的机器学习技术以及我从数据中得出的发现的步骤。

2.探索性数据分析

我在这部分的代码和结果可以在这里找到。

数据集带有两个 CSV 文件，都可以从我的 GitHub 页面访问。我首先用 Python 将它们导入并连接成一个熊猫数据帧。除了我们感兴趣的 11 个特征之外，进行了一些重新排列以移除列:

• 臭氧
• 硫化氢
• 总挥发性有机化合物
• 二氧化碳
• PM 1
• PM 2.5
• PM 10
• 温度(内部和外部)
• 湿度(内部和外部)。

在美国东部时间(GMT-4)时区，时间戳从 2020 年 5 月 26 日到 6 月 9 日(总共 14 整天)。通过减法，发现每个读数之间有不同的间隔，从 7 秒到 3552 秒不等。下面的表 1 中列出了前 5 个频繁的时间间隔，其中大多数都接近 59 和 60 秒，因此可以得出结论，传感器每分钟都在读取。然而，如果不涉及故意干扰，阅读间隔的不一致性可能值得研究，因为它可能会在未来的时间序列分析中造成麻烦。

表 1:传感器测量的前 5 个时间间隔

对于每个特征，时间序列数据处于不同的尺度上，因此为了更好的可视化和机器学习效率，它们被归一化。然后对它们进行绘图和视觉检查，以发现任何有趣的模式。

2.1 模式变化

一些特征似乎在特定时间点共享相似的模式变化。三个最重要的因素(外部温度、外部湿度和臭氧)如下图 1 所示。可以清楚地看到，用粉红色突出显示的区域倾向于具有平坦的信号，而未突出显示的区域是正弦曲线。

图 1:温度(外部)、湿度(外部)和臭氧的读数。它们在相同的时间点经历相似的模式变化。

根据常识，室外温度在中午达到最高点，并在晚上下降，我开始怀疑在这 14 天期间不同的测试环境的可能性。为了测试这个想法，从加拿大天气统计中查询多伦多天气数据[1]。温度和相对湿度被叠加，并与该数据集中的外部温度和湿度进行比较。该图如图 2 所示。可以看出，实际温度和湿度以正弦方式波动。温度和湿度读数的大部分与天气数据有很好的相关性，而用粉红色突出显示的区域保持相对不变。我没有获得任何关于测量环境的相关信息，但是从图中可以合理地推断出该装置在 14 天期间在室内和室外环境之间进行了重新定位。这也将在第 3.3 节的自动异常检测中进行测试。

图 2:温度(外部)和湿度(外部)与多伦多天气数据叠加。与其他区域相比，粉红色突出显示的区域保持相对不变。

2.2 特征之间的相关性

相关性是一种研究两个定量的连续变量之间关系的技术，以表示它们之间的相互依赖性。在不同的相关技术中，皮尔逊相关是最常用的一种，它测量两个变量之间的关联强度。其相关系数从-1 到 1，其中 1 代表最强的正相关，-1 代表最强的负相关，0 代表无相关。计算每对数据集之间的相关系数，并绘制成热图，如表 2 所示。所选特征的散点图也绘制在图中并附在附录部分。

首先要注意的是，PM 1、PM 2.5 和 PM 10 之间的相关性很高，这意味着它们总是以相同的方式波动。臭氧与二氧化碳负相关，与温度(内部)和温度(外部)正相关。另一方面，令人惊讶的是，在温度(内部)和温度(外部)之间没有发现任何显著的相关性，这可能是由于仪器具有优异的隔热性能。然而，由于没有提供相关知识，因此无法对这一发现的合理性做出结论。除臭氧外，温度(内部)也与二氧化碳、硫化氢和三种颗粒物指标呈负相关。相反，温度(外部)与湿度(内部)和三种颗粒物测量值正相关，而与湿度(外部)负相关，正如图 1 中的时间序列图所示。

表 2:特征间皮尔逊相关系数的热图。

3.异常检测和模式识别

在本节中，将基于数据集检查各种异常检测方法。数据没有标签，因此没有知识或分类规则来区分“系统故障”、“外部事件”和其他。也没有提供仪器和实验的任何细节。因此，我在这一部分的结果可能会偏离预期，但我会尽最大努力做出假设，定义问题，然后根据我的个人经验完成它们。该部分由三部分组成:点异常检测、集体异常检测和聚类。

3.1 点异常检测(系统故障)

点异常或全局异常值是那些完全在通常信号范围之外的数据点，没有任何近邻的支持。这通常是由人为或系统错误引起的，需要在数据清理过程中移除，以便在预测建模中获得更好的性能。在这个数据集中，通过假设“系统故障”等同于这种点异常，有几个特征值得研究，例如图 3 所示的例子。

图 3:自动标记之前(左)和之后(右)的时间序列信号。点异常(异常值)标记为红色。从上到下:湿度(内部)、总挥发性有机化合物和二氧化碳。

这里，从湿度(内部)到总挥发性有机化合物到二氧化碳，每一个都代表了点异常检测任务的独特复杂性。在第一个例子中，三个异常值位于 0 级，因此一个简单的布尔过滤器可以标记这些数据点。在第二种情况下，异常值明显偏离我们感兴趣的信号，因此可以使用线性阈值来分离出异常值。这两种情况都很容易实现，因为它们可以通过纯粹基于经验的方法来完成。当涉及到第三种情况时，不可能使用线性阈值来分离出异常值，因为即使它们偏离其邻居，这些值也可能不像在其他时间点的通常信号那样大。

对于这种情况，有很多方法可以处理。一种简单的方法是计算每个时间点的滚动平均值或中值，并测试实际值是否在通过从中心线加减一定波动范围计算的预测区间内。因为我们在处理异常值，所以滚动中位数更稳健，所以在这个数据集中使用。

图 4:放大的二氧化碳时间序列图，带有滚动中值和预测区间

从图 4 中可以更清楚地看到这种方法:基于滚动中值设置等距预测区间。这里，滚动窗口设置为 5，这意味着对于每个数据点，我们取其最近的四个邻居，并计算中位数作为预测的中心。则围绕中心填充 0.17 的预测间隔。任何外部点都被视为异常值。

使用这种方法检测点异常是简单而有效的。但是，它也有不足之处，在处理更复杂的数据时可能不够可靠。在该模型中，有两个参数:滚动窗口大小和预测区间大小，这两个参数都是通过给定数据的实验手动定义的。我们基本上解决了从情况 2 到情况 3 时遇到的问题，如图 3 所示，方法是根据时间对信号和异常值之间的分类边界启用自调整能力。然而，带宽是固定的，因此在点异常的定义随时间变化的情况下，它变得不太有用。例如，超级便宜机票的定义在正常工作日和假日季节可能完全不同。

这就是机器学习发挥作用的时候，模型可以从数据中学习，以便在时间发生变化时自行调整参数。它将知道在特定时间点何时数据点可以被分类为点异常或不被分类为点异常。然而，我不能为这个数据集在这个任务上构建这样一个有监督的机器学习模型，因为前提是有标记的数据，而我们的没有。我们仍然建议在未来沿着这条路走下去，因为有了这样的模型，无论系统或数据有多复杂，都会产生最准确的结果。

即使有监督的学习方法在这个任务中不起作用，也有无监督的学习技术是有用的，例如聚类，也在 3.3 小节中讨论。聚类可以通过使用相似性度量(如向量空间中数据点之间的距离)对未标记的数据进行分组。然后通过选择远离聚类中心的异常点来区分异常点。然而，为了在该数据集中使用聚类进行点异常检测，我们必须遵循这样的假设，即外部事件是用多个特征的范围来定义的，而不是分别对待每个单独的时间序列。关于这一点的更多细节将在以下两个小节 3.2 和 3.3 中讨论。

3.2 集体异常检测(外部事件)

如果我们将“系统故障”定义为点异常，那么“外部事件”有两个方向。其中之一是将其定义为出现在每个时间序列信号中的集体异常。集体异常的概念与点异常相反。点异常是与通常信号有很大偏差的不连续值，而集体异常通常是连续的，但值超出预期，如在某些时间点显著增加或减少。在这一小节中，所有 11 个特征被分别视为一个时间序列。接下来的任务是找出它们各自发生的突变。

对于这样的问题，一种典型的方法是调整时间序列预测的用法:我们可以将模型拟合到某个时间段之前，然后预测该时间段之后的值。然后将实际值进行比较，看它是否落入预测区间。它与前一小节中使用的滚动中值方法非常相似，但这里只使用了以前的时间点，而不是使用两个方向的邻居。

该模型也有不同的选项。像 SARIMA 这样的传统时间序列预测模型是一个很好的选择，但是该模型可能不够复杂，无法适应我在第 2.1 节和第 3.3 节中提到的“模式”。另一种选择是为时间序列训练一个监督回归模型，这种模型现在被广泛使用。

想法很简单:使用滑动窗口的概念从时间序列中提取特征，如表 3 和图 5 所示。滑动窗口大小(蓝色)设置为与所需的特征号 k 相同。然后，对于时间序列中的每个数据点(橙色)，特征是从滞后 1 到滞后 k 之前的数据点值。因此，具有 N 个样本的时间序列能够被转换成 N-k 个观察值和 k 个特征的表。接下来，通过实现“正向链接”的概念，每个点由使用从指数 0 到 k-1 的观察值训练的回归模型预测。除了主回归模型之外，还训练了两个具有不同显著性水平的分位数回归变量来预测预测区间的上限和下限，这样我们就能够判断实际值是高于还是低于区间范围。

表 3:通过滑动窗口的特征提取算法

图 5:滑动窗口和正向链接的概念[2]。使用滑动窗口从时间序列中提取特征和目标。训练过程基于“正向链接”。

这个方法应用于给定的数据集，下面的图 6 显示了一个例子。为了加快训练速度，每小时对臭氧时间序列进行采样，提取特征，并使用 Scikit-Learn 将其输入三个梯度增强回归器模型(1 个主回归器和 2 个分位数回归器)。选择显著性水平，使得预测水平代表 90%的置信区间(在图 6 顶部显示为绿色)。然后将实际值与预测区间进行比较，并在图 6 底部用红色(意外增加)和蓝色(意外减少)标记。

图 6:通过应用分位数回归方法，用臭氧时间序列数据(每小时采样)标记结果。顶部:数据和预测区间；底部:带有显示在预测间隔之上或之下的标志的数据。

结果可能还不是非常令人印象深刻，因为围绕回归模型选择和超参数微调还需要做更多的工作。然而，它已经显示出应对这些突然增加和减少的能力。使用机器学习模型的一个好处是，当输入数据时，模型会自己学习和进化。从图 6 底部可以看出，最后三个山丘(大约 240 小时后)的标记点比之前的少。这不仅是因为幅度较小，而且是因为模型正在从以前的经验中学习，并开始适应“思想”,它现在处于“山区”,应预料到周期性波动。因此，不难得出结论，如果输入更多的数据实例，模型性能会越来越好。

除了这种分位数回归模型，深度学习模型如 LSTM 可能能够实现更好的性能。长短期记忆(LSTM)是一种专门的人工递归神经网络(RNN)，由于其反馈连接的特殊设计，它是序列建模的最新选择之一。但是，设置和微调网络架构需要更长的时间和精力，这超出了本项目的时间限制，因此本报告中未将其作为可展示的内容。

另一方面，正如我在前面的章节中提到的，所提供的数据没有显示哪些数据点被认为是异常的标签。通过限制模型的选择和性能，它确实给本小节中讨论的集体异常检测任务带来了困难。未来如果提供一些标签，就变成了半监督或监督学习问题，这样就变得更容易达到更好的效果。

3.3 聚类和模式识别(外部事件)

正如我在上面 3.2 小节中提到的，识别“外部事件”可以从两个方向着手:一个是分别处理每个时间序列并监控传感器信号发生的任何意外变化，另一个是假设事件同时影响多个特征，这样我们希望通过查看不同特征中显示的独特特征来区分事件。在这种情况下，如果我们已经标记了数据，这将是一个常见的分类问题，但即使没有标签，我们仍然能够通过聚类来处理。

聚类是一种无监督的机器学习技术，它根据特征发现数据之间的相似性，并将相似的数据对象分组到聚类中。它可以作为一个独立的工具来深入了解数据分布，也可以作为其他算法的预处理步骤。有许多不同的聚类方法，这里我使用两种最常用的方法:K-Means 和 DBSCAN。

K-means 是用于聚类的划分方法之一。它将对象随机划分为非空子集，并不断添加新对象和调整质心，直到在优化每个对象和质心之间的距离平方和时达到局部最小值。另一方面，基于密度的噪声应用空间聚类(DBSCAN)是一种基于密度的方法，其中聚类被定义为密度连接点的最大集合[3]。

图 7:投影到前两个主成分的聚类结果。左:K-Means；右图:DBSCAN

图 8:投影到前 3 个主成分的聚类结果。左:K-Means；右图:DBSCAN

图 9:温度时间序列的聚类结果(外部)。左:K-Means；右图:DBSCAN

主成分分析(PCA)是一种降维技术，它创造了新的不相关变量，以增加可解释性和最小化信息损失。在本项目中，对标准化数据应用 K-means 和 DBSCAN 算法后，执行 PCA，并使用前 2 个和 3 个主成分在 2D(图 7)和 3D(图 8)中绘制聚类结果。此外，为了从另一个角度查看聚类结果，制作了带标签的时间序列图，温度(外部)图如图 9 所示。

从图中可以清楚地看出，这两种方法都能够区分我在第 2.1 节中提到的室内/室外模式变化。主要区别是基于划分的 K-means 方法对昼夜交替引起的星等变化更敏感。数据集中的许多变量都受到同时发生的这种明显的正弦变化的影响，包括温度(外部和内部)、湿度(外部和内部)和臭氧。k 均值倾向于区别对待波峰和波谷。另一方面，基于密度的 DBSCAN 不太关心星等差异，而是更关注密度分布。因此，它将整个正弦部分聚集成一个质量云，如图 7 和图 8 所示。

在现阶段不可能评论哪种聚类方法比另一种更好，因为它们有足够的独特性来满足不同的兴趣。如果我们对区别对待正弦信号的高低部分更感兴趣，我们将使用 K 均值；如果只想区分室内/室外模式，那么 DBSCAN 更好。另外，由于是无监督的学习任务，除了可视化和凭经验判断，没有办法量化模型之间的性能。在未来，如果提供一些标记数据，结果可以变成半监督学习任务，并且可以获得更多关于模型选择的直觉。

4.结论

在本文中，我将简要介绍探索性数据分析和相关性分析的方法和结果，以及用于点异常检测、集体异常检测和聚类的三种不同建模管道的构造。

在探索性数据分析部分，发现传感器读数时间间隔变化很大。即使它们中的大多数包含在一分钟左右，不一致性问题仍然值得研究，因为它可能会降低分析任务的效率和性能。我们发现测量环境受时间序列图变化的影响，后来通过与实际的多伦多天气数据以及聚类结果保持一致而得到保证。此外，还对特征之间的相关性进行了研究和举例。出现了一些困惑，如温度(内部)和温度(外部)之间的奇怪关系，这需要通过实验或设备本身来研究。

在异常检测部分，由于“系统故障”和“外部事件”没有明确定义，我将项目分成三个不同的任务。点异常被定义为严重偏离和不连续的数据点。这里使用了滚动中值方法来成功地自动化标记这种点异常的过程。另一方面，集体异常被定义为偏离的数据点集合，通常被视为突然增加或减少。这项任务是通过从时间序列数据中提取特征，然后训练回归模型来完成的。还使用 K-mean 和 DBSCAN 对数据集执行聚类，这两种方法都发挥了它们的优势，并通过利用它们相似和不相似的特征成功地对数据进行聚类。

本项目中介绍的所有异常检测模型都只是原型，没有广泛的模型部分和微调。如果付出更多的努力，并通过获得更多的数据知识，它们中的每一个都有很大的潜力进化成更好的形式。对于点异常，有更多基于机器学习的异常检测技术，如隔离森林和局部异常因子，以适应更复杂的数据形式。对于集体异常，最先进的 LSTM 值得投入精力，特别是在时间序列数据和序列建模方面。对于聚类来说，还有许多其他的方法，比如层次聚类和基于网格的聚类。他们有能力获得类似的出色表现。

当然，这些未来的方向都是在没有标注数据的前提下建议的。如果有经验的工程师或科学家能够就哪些类型的数据被视为“系统故障”或“外部事件”给出他们的见解，那么通过将任务转化为半监督或监督学习问题，肯定会取得更令人兴奋的进展，届时将有更多的工具可供选择。

感谢您的阅读！如果你喜欢这篇文章，请关注我的频道和/或 成为我今天的推荐会员 (非常感谢🙏).我会继续写下去，分享我关于数据科学的想法和项目。如果你有任何问题，请随时联系我。

[## 通过我的推荐链接加入 Medium 周(乔)徐

阅读周(Joe)徐(以及媒体上成千上万的其他作家)的每一个故事。您的会员费直接支持…

zhouxu-ds.medium.com](https://zhouxu-ds.medium.com/membership)

关于我

我是赛诺菲的数据科学家。我拥抱技术，每天都在学习新技能。欢迎您通过媒体博客、 LinkedIn 或 GitHub 联系我。我的观点是我自己的，而不是我雇主的观点。

请看我的其他文章:

• 使用 Python 和 Flask-RESTful 为机器学习模型构建 REST API
• 使用 Elasticsearch 的实时类型预测搜索(AWS OpenSearch)
• 理解分类问题中的 Sigmoid、Logistic、Softmax 函数和交叉熵损失(对数损失)
• 利润最大化的贷款违约预测
• 使用 Berka 数据集进行贷款违约预测

参考

[1]加拿大天气统计:https://www.weatherstats.ca/

[2]时间序列机器学习框架:https://towards data science . com/Time-Series-Machine-Learning-regression-Framework-9ea 33929009 a

[3] J. Han，M. Kamber 和 J. Pei，数据挖掘:概念和技术，2011 年第 3 版。

附录—所选特征的散布矩阵

超越测试数据的时间序列预测

原文：https://towardsdatascience.com/time-series-prediction-beyond-test-data-3f4625019fd9?source=collection_archive---------3-----------------------

许多(如果不是全部的话)与时间序列预测相关的例子仅用于测试数据，缺乏如何预测真实未来日期的例子。这个简单的例子展示了如何在测试数据集之外进行预测。

来源:Pixabay

我当时的任务是构建一个大规模时间序列预测解决方案。我最终在单一解决方案中使用了多种方法的组合——预言家、ARIMA 和 LSTM 神经网络(运行在 Keras/TensorFlow 之上)。使用 Prophet ( 用 Flask 提供 Prophet 模型—预测未来)和 ARIMA 可以直接计算未来日期的预测，两者都提供了返回给定未来范围的预测的函数。如果你是新手，LSTM 的情况就不一样了——这将需要大量的时间来研究如何预测真实的未来日期(大多数例子只展示了如何根据测试数据集进行预测)。

我发现了一个很好的例子，它帮助我解决了我的任务——一个使用长短期记忆(LSTM)网络进行时间序列预测的快速例子。在这篇文章中，我将展示如何预测洗发水销售月度数据，主要基于上面例子中的代码。

可以从 Kaggle 下载一个包含洗发水月销售数据的数据集。

LSTM 模型是用 TensorFlow 2.0 和 Keras 建立的。

首先，我们需要获取和解析时间序列数据:

def parser(x):
    return pd.datetime.strptime('190'+x, '%Y-%m')df = pd.read_csv('shampoo.csv', parse_dates=[0], index_col=0, date_parser=parser)
df.tail()

对于 LSTM 模型训练，必须对数据进行缩放，这是通过来自 sklearn 库的 MinMaxScaler 完成的:

train = dfscaler = MinMaxScaler()
scaler.fit(train)
train = scaler.transform(train)

时间序列预测是一项序列预测任务。日期对于 ML 算法来说并不重要，它试图识别过去和未来数据之间的依赖模式。有许多方法可以将原始序列数据转换成适合 LSTM 神经网络训练的数据。我更喜欢使用 Keras TimeseriesGenerator(在这里阅读更多信息— 如何在 Keras 中使用 TimeseriesGenerator 进行时间序列预测)。我们将预测 12 个月的洗发水销售额，因为这将使用过去 12 个月的时间窗口来确定预测:

n_input = 12
n_features = 1
generator = TimeseriesGenerator(train, train, length=n_input, batch_size=6)

LSTM 模特培训。对于输入和目标，我们将使用整个时间序列数据集，模型将在整个数据集上进行训练。这将允许为训练集中可用的最后日期之后的日期生成预测:

model = Sequential()
model.add(LSTM(200, activation='relu', input_shape=(n_input, n_features)))
model.add(Dropout(0.15))
model.add(Dense(1))optimizer = keras.optimizers.Adam(learning_rate=0.001)
model.compile(optimizer=optimizer, loss='mse')history = model.fit_generator(generator,epochs=100,verbose=1)

一旦模型训练完毕，我们就可以生成未来日期的预测。这里的想法是——使用最近的 12 个步骤，为未来的一个步骤生成一个预测。向数组中添加新的预测，从同一数组中删除第一个条目，并用 12 步的更新数组预测下一步。

pred_list = []batch = train[-n_input:].reshape((1, n_input, n_features))for i in range(n_input):   
    pred_list.append(model.predict(batch)[0]) 
    batch = np.append(batch[:,1:,:],[[pred_list[i]]],axis=1)

创建包含未来月份的数据框，我们将为这些月份分配计算的预测数据:

add_dates = [df.index[-1] + DateOffset(months=x) for x in range(0,13) ]
future_dates = pd.DataFrame(index=add_dates[1:],columns=df.columns)

创建包含日期和预测数据的数据框。最后，与训练数据框架合并，以绘制显示训练数据和真实未来预测的图表:

df_predict = pd.DataFrame(scaler.inverse_transform(pred_list),
                          index=future_dates[-n_input:].index, columns=['Prediction'])df_proj = pd.concat([df,df_predict], axis=1)

结果显示为 Plotly:

plot_data = [
    go.Scatter(
        x=df_proj.index,
        y=df_proj['Sales'],
        name='actual'
    ),
    go.Scatter(
        x=df_proj.index,
        y=df_proj['Prediction'],
        name='prediction'
    )
]plot_layout = go.Layout(
        title='Shampoo sales prediction'
    )
fig = go.Figure(data=plot_data, layout=plot_layout)
pyoff.iplot(fig)

预测的真实未来销售额，以黄色显示:

源代码可以在 GitHub 上获得。