TowardsDataScience 2023 博客中文翻译(一百三十七)

原文:TowardsDataScience

协议:CC BY-NC-SA 4.0

修复 Prophet 的预测问题

原文:towardsdatascience.com/fixing-prophets-forecasting-issue-b473afe2cc70

第一步:限制疯狂的趋势

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传 Tyler Blume

·发布于Towards Data Science ·10 分钟阅读·2023 年 1 月 24 日

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由Hunter Haley提供,来源于Unsplash

目前,Prophet 预测准确性问题已不是秘密。它在多个基准和预测竞赛中一再交付糟糕的结果。然而,它仍然是最常用的预测算法之一……

所以…是时候通过一些小的调整来解决困扰它的问题,并(希望)提高它的预测准确性了。

TSUtilities 项目

预言者的趋势问题

奇怪的是,Prophet 的一个主要吸引力也是其核心弱点之一。它确实提供了一个引人注目的趋势,包含了变更点和线性段,分析起来非常易于消化。但是,有时,这种拟合趋势的方式既过度拟合又不足拟合——在序列末尾的一个简单水平偏移可能会让趋势在预测范围内无限延伸,同时在残差中留下大量有用信号。

这篇文章旨在解决第一个问题:不受限制的疯狂趋势

为了说明这个问题,我们将使用 M4 数据集。这些数据集都是开源的,可以在 M-competitions 的github上找到。

首先,让我们看一下每周数据集,特别是这个数据集中第 52 个时间序列。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由作者提供

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由作者提供

如果你之前使用过 Prophet,这些图表应该很熟悉。

在这里,我们对 100 的时间范围进行预测——远远超出了 M4 竞赛中对该系列的处理。但这展示了 Prophet 可能会如何失控,趋势从 7000 骤降至约 5500。这意味着时间序列在几年的过程中完全崩溃。从视觉上来看,我明白模型捕捉到的内容,并且在预测范围的近期,遵循这一点可能是可以的。然而,从长期来看,这种趋势可能会导致灾难性的结果。这个想法很重要,所以我会用粗体文字重新写一下:

短期内跟随趋势,长期内控制趋势。

在我们进入详细内容之前,让我们澄清一些术语。我见过‘dampen’、‘dampened’、‘dampening’和‘damped’这些词在不同来源和时代中交替使用。我知道有一个‘正确的’词,我曾经谷歌过,但我会根据自己的感觉使用它们。只要知道它们在本文中意义相同,如果你有偏好,务必告诉我!

现在开始修复问题。

修复趋势

通常,在时间序列背景下,‘damped’趋势出现时讨论的是双指数平滑。对我们来说,我们将采用更简单的方法,更符合‘经验法则’,而不是严格的统计方法。比如——指数衰减。

使用指数衰减,我们将从原始趋势非常接近开始,但在进入时间范围的过程中会逐渐远离它。同时,我们将接近一个‘目标’,该目标将作为一个连贯的参数提供。

让我们看一个简单的例子,首先是一个漂亮的曲线:

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

sns.set_style('darkgrid')
y = np.linspace(0, 100, 100)
plt.plot(y)
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

作者提供的图片

真是漂亮的一条线!

接下来,我们将其分解为在处理时间序列时预期的不同部分。这里我们假设趋势等于实际信号,没有噪音或季节性。

y_train = y[:80]
future_y = y[80:]
future_trend = future_y

现在我们将使用‘future_trend’来抑制这个信号。为此,我们将安装一个名为‘TSUtilities’的包,它是我在其他包ThymeBoostLazyProphetTimeMurmur中使用的各种工具的集合,我希望将其集中管理。

pip install TSUtilities

仍在整理和结构化,但 TrendDampen 在这里,并可以像这样导入:

from TSUtilities.TSTrend.trend_dampen import TrendDampen

创建类时传递的两个参数是:

damp_factor:一个介于 0 和 1 之间的浮点数,其中 0 表示‘完全抑制’,1 表示‘不抑制’。

damp_style:如何‘抑制’趋势,通过平滑来‘平滑’地抑制趋势,使用指数衰减。

dampener = TrendDampen(damp_factor=.7,
                       damp_style='smooth')
dampened_trend = dampener.dampen(future_trend)

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

作者提供的图片

在这里我们使用 0.7 的 damp_factor,这意味着新趋势实现了旧趋势大约 70%的效果。我们在这里唯一需要传递给dampen方法的是预测的趋势组件。

正如前面提到的,我们从接近原始趋势开始,但会越来越远。

现在这对我们帮助不大,因为我们传入了一个硬编码的 damp_factor 值,但我们可以利用这个方法中内置的一些进一步逻辑。如果我们传入拟合的组件:

  1. 训练实际值

  2. 拟合的趋势成分

  3. 拟合的季节性成分

(注:所有这些都是由 Prophet 提供的)

并将 damp_factor 改为‘auto’——方法将为你选择一个合适的参数。它通过评估拟合趋势的强度来实现这一点——这个趋势越强,我们就越愿意在长远中信任它。

注:这里的强度是通过 Wang, Smith, & Hyndman¹定义的趋势强度度量来确定的。

我们使用的这个‘经验法则’将根据算法的不同而有所不同。如果我们将逻辑应用于 ETS 方法,那么按照这些定义,趋势将会异常强烈,我们会过度信任它。但对于像 Prophet 这样的方法,或任何其他有某种确定性趋势的方法,它在实践中效果良好。

所以让我们看看这个方法如何处理我们简单的线。首先,我们需要一些拟合趋势和季节性的度量。

trend = y_train
seasonality = np.zeros(len(y_train))

那么我们将把这些值传递给 dampen 方法。

dampener = TrendDampen(damp_factor='auto',
                       damp_style='smooth')
dampened_trend = dampener.dampen(future_trend,
                                 y=y_train,
                                 trend_component=trend,
                                 seasonality_component=seasonality)
plt.plot(future_trend, color='black', alpha=.5,label='Actual Trend', linestyle='dotted')
plt.plot(dampened_trend, label='Damped Trend', alpha=.7)
plt.legend()
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由作者提供

正如你所见,‘auto’决定完全不抑制趋势!

这是一个很好的检查——虽然我相信你可以很容易地打破这个逻辑。

让我们看看这个序列的各种 damp_factors:

for damp_factor in [.1, .3, .5, .7, .9, 'auto']:
    dampener = TrendDampen(damp_factor=damp_factor,
                           damp_style='smooth')
    dampened_trend = dampener.dampen(future_trend,
                                     y=y_train,
                                     trend_component=trend,
                                     seasonality_component=seasonality)
    plt.plot(dampened_trend, label=damp_factor, alpha=.7)

plt.plot(future_trend, color='black', alpha=.5,label='Actual Trend', linestyle='dotted')
plt.legend()
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由作者提供

很酷!

再次,我们有一个很好的、易于理解的参数,并且有一些‘auto’逻辑看起来还不错。

现在进入主要内容:修复 Prophet 的趋势。

在这个示例中,我使用了 M4 的每周数据集中第 52 个时间序列,并将预测 100 期。这个预测期虽然较长,但能够很好地展示问题。此外,我将使用一个函数,该函数接收 Prophet 的输出并抑制趋势,以返回这些修正后的预测。这个函数看起来像这样:

def dampen_prophet(y, fit_df, forecast_df):
    """
    A function that takes in the forecasted dataframe output of Prophet and
    constrains the trend based on it's percieved strength'

    Parameters
    ----------
    y : pd.Series
        The single time series of actuals that are fitted with Prophet.
    fit_df : pd.DataFrame
        The Fitted DataFrame from Prophet.
    forecast_df : pd.DataFrame
        The future forecast dataframe from prophet which includes the predicted trend.

    Returns
    -------
    forecasts_damped : np.array
        The damped trend forecast.

    """
    predictions = forecast_df.tail(len(forecast_df) - len(fit_df))
    predicted_trend = predictions['trend'].values
    trend_component = fit_df['trend'].values
    if 'multiplicative_terms' in forecast_df.columns:
        seasonality_component = fit_df['trend'].values * \
                                fit_df['multiplicative_terms'].values
        dampener = TrendDampen(damp_factor='auto',
                                damp_style='smooth')
        dampened_trend = dampener.dampen(predicted_trend,
                                         y=y,
                                         trend_component=trend_component,
                                         seasonality_component=seasonality_component)
        forecasts_damped = predictions['additive_terms'].values + \
                           dampened_trend + \
                           (dampened_trend * \
                           predictions['multiplicative_terms'].values)
    else:
        seasonality_component = fit_df['additive_terms'].values
        dampener = TrendDampen(damp_factor='auto',
                                damp_style='smooth')
        dampened_trend = dampener.dampen(predicted_trend,
                                         y=y,
                                         trend_component=trend_component,
                                         seasonality_component=seasonality_component)
        forecasts_damped = predictions['additive_terms'].values + dampened_trend
    return forecasts_damped

但这个函数也在 TSUtilities 中,我们可以直接使用以下方法导入:

from TSUtilities.functions import dampen_prophet

现在让我们看看正常的 Prophet 输出和抑制后的输出有什么不同:

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由作者提供

我们可以看到,Prophet 在预测范围内的最近趋势变化点远低于任何先前的值,而抑制后的趋势则显得更加保守。

从趋势的角度来看,我们可以看到究竟发生了什么,以及方法改变了什么。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由作者提供

注意,我们的方法并没有改变拟合的值。我们仅仅是旨在收敛未来的趋势成分。

最后,我们看到两种方法的 SMAPE——使用衰减趋势的准确性有了相当显著的提高:

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由作者提供

在 M4 上的基准测试

到目前为止,我们只看了一个单一的例子,但现在让我们看看 M4 的两个数据集——每周每日数据集。提醒一下,这些数据集都是开源的,并且作为还可以的基准,适合在时间序列领域进行尝试。

让我们继续导入每周的数据集:

import matplotlib.pyplot as plt
import numpy as np
from tqdm import tqdm
import pandas as pd
from prophet import Prophet
import seaborn as sns
sns.set_style('darkgrid')

train_df = pd.read_csv(r'm4-weekly-train.csv')
test_df = pd.read_csv(r'm4-weekly-test.csv')
train_df.index = train_df['V1']
train_df = train_df.drop('V1', axis = 1)
test_df.index = test_df['V1']
test_df = test_df.drop('V1', axis = 1)

接下来,让我们定义用来评估预测的 SMAPE 函数:

def smape(A, F):
    return 100/len(A) * np.sum(2 * np.abs(F - A) / (np.abs(A) + np.abs(F)))

现在我们已经定义了数据和度量指标,让我们遍历数据集,使用 Prophet 生成预测——一个标准预测和一个带有衰减趋势的预测。

seasonality = 52
no_damp_smapes = []
damp_smapes = []
naive_smape = []
j = tqdm(range(len(train_df)))
for row in j:
    y = train_df.iloc[row, :].dropna()
    y = y.iloc[-(3*seasonality):]
    y_test = test_df.iloc[row, :].dropna()
    #create a random datetime index to pass to Prophet
    ds = pd.date_range(start='01-01-2000',
                       periods=len(y) + len(y_test),
                       freq='W')
    ts = y.to_frame()
    ts.columns = ['y']
    ts['ds'] = ds[:len(y)]
    j.set_description(f'{np.mean(no_damp_smapes)}, {np.mean(damp_smapes)}')
    prophet = Prophet()
    prophet.fit(ts)
    fitted = prophet.predict()

    # create a future data frame
    future = prophet.make_future_dataframe(freq='W',periods=len(y_test))
    forecast = prophet.predict(future)

    #get predictions and required data inputs for auto-damping
    predictions = forecast.tail(len(y_test))
    predicted_trend = predictions['trend'].values
    trend_component = fitted['trend'].values
    seasonality_component = fitted['additive_terms'].values
    forecasts_no_dampen = predictions['yhat'].values
    forecasts_damped = dampen_prophet(y=y.values,
                                      fit_df=fitted,
                                      forecast_df=forecast)

    #append smape for each method
    no_damp_smapes.append(smape(y_test.values, forecasts_no_dampen))
    damp_smapes.append(smape(y_test.values, forecasts_damped))
    naive_smape.append(smape(y_test.values, np.tile(y.iloc[-1], len(y_test))))
print(f'Standard Prophet {np.mean(no_damp_smapes)}')
print(f'Damped Prophet {np.mean(damp_smapes)}')
print(f'Naive {np.mean(naive_smape)}')

结果如何?

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由作者提供

衰减方法有效!但 SMAPE 仅减少了约 2.5%。虽然不足以大书特书,但还是一个不错的提升。值得记住的是,我们并没有真正失去 Prophet,只是改变了它的趋势。如果 Prophet 是你预测过程中的一个重要部分,那么这可以看作是锦上添花。

不过,对于每日数据集,结果确实有所改善。重新运行过程得到:

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由作者提供

,SMAPE 减少了 7.8%!

与其他在 M4 上进行基准测试的模型比较这些结果时,我们发现即使有这些改进,Prophet 仍然无论如何都产生较差的结果。我想一个理由是我们没有调整超参数。不管怎样——我不会使用 Prophet。

我不会在集成模型中使用 Prophet。

我不会使用 Prophet 进行多变量预测。

我绝对不会输入:import Prophet

但这一系列的重点是尝试增强Prophet。下次——如果有下次的话——我们会有另一个技巧来实现这个目标。

结论

我们展示了基于经验法则的趋势衰减过程在提高 Prophet 预测准确性方面的有效性。该过程非常灵活,可以与任何输出一起使用(尽管 dampen_prophet 函数本身不灵活)。我们在 M4 系列的两个数据集上观察到 SMAPE 减少了约 2.4%和约 7.8%。

如果你喜欢这篇文章,你可能会喜欢我其他的一些文章!

## LazyProphet: Time Series Forecasting with LightGBM

一切都与特征有关

[towardsdatascience.com ## Time Series Forecasting with ThymeBoost

一种梯度提升时间序列分解方法

towardsdatascience.com ## 梯度增强 ARIMA 用于时间序列预测

提升 PmdArima 的 Auto-Arima 性能

towardsdatascience.com

剧透警告

本系列的终点将完全摆脱 Prophet —— 思路类似于 Nixtla 使用的 StatForecast 包。与 Nixtla 不同的是,我们将保留相同的时间序列分解和易于理解的线性变点趋势。这将通过 ThymeBoost 实现。

当然,您也可以通过我的推荐链接注册 Medium:

[## 通过我的推荐链接加入 Medium - Tyler Blume

阅读 Tyler Blume 的每一个故事(以及 Medium 上成千上万其他作家的故事)。您的会员费将直接支持……

medium.com](https://medium.com/@tylerblume/membership?source=post_page-----b473afe2cc70--------------------------------)

参考文献

  1. Wang, X., Smith, K. A., & Hyndman, R. J. (2006). 基于特征的时间序列数据聚类。数据挖掘与知识发现, 13(3), 335–364.

Flapjax: 使用 Plotly 和 Flask 进行网络数据可视化

原文:towardsdatascience.com/flapjax-data-visualization-on-the-web-with-plotly-and-flask-465090fa3fba

使用 Plotly 和 Flask 构建一个数据可视化网页,并用一些 UI 组件使其互动

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传 Alan Jones

·发表于 Towards Data Science ·17 分钟阅读·2023 年 11 月 17 日

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由 Mae Mu 提供,Unsplash

构建数据可视化应用程序的最佳框架是什么?是 Streamlit 还是 Dash?或者你可以用 MercuryVoilá 将 Jupyter Notebook 转换为网络应用程序?

所有这些都是创建应用程序的好方法,而且相对容易入门。但通常,容易入门的东西会随着你变得更有冒险精神而变得稍微复杂一些。因此,我将试图说服你,回到基础,使用 Python 服务器代码和 HTML 页面作为用户界面,其实并没有看起来那么令人生畏。

我们可以使用相当多的模板代码和模板构建引人入胜的交互式应用程序,这意味着你可以仍然集中精力在你的 Python 代码上,对 HTML 和 Javascript 的接触是最小的。我称这种方法为Flapjax——稍后我会解释原因。

创建一个简单的 Python 网络应用程序的一种最简单方法是使用 Flask,这正是我们要做的,我们将创建一个看起来像下面图片中的应用程序。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

一个示例互动应用程序

Flask 框架

Flask 是一个用于开发 Web 应用的极简框架。在 Flask 应用中,网页通常是由模板和 Python 代码提供的数据构建的——这些数据可以是形成网页内容的文本或图形。结果会被发送给用户以在浏览器中显示。

下图展示了一个交互式应用的基本结构。当应用运行时,Python 部分在服务器上执行,并将数据传递给在浏览器中运行的 HTML。网页上的用户输入会传回 Python 代码,Python 代码可能会发送更多数据以更新 HTML 内容,例如用户选择的新图表。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

这是创建 Web 应用程序的最简单方法吗?

在我看来,将用户界面设计与程序逻辑分离确实使生活变得更轻松。但如果你习惯于在 Streamlit 或 Jupyter Notebooks 中构建应用,你可能会发现有一定的学习曲线。然而,一旦你采用了一个基本应用的模式,创建新的应用会容易得多。

因此,我们将使用 Flask 开发一个数据可视化应用,并且我们还将使用 Jinja 模板来定义我们的 HTML 页面——虽然实际出现在这些页面中的数据将由我们的 Python 代码定义。

要制作一个交互式用户界面,我们需要一些 UI 组件和一点 JavaScript,但我们会看到这基本上是可以在未来应用中重用的样板代码。

我们还将使用 Bootstrap 5 UI 组件,因为你为什么要满足于一个看起来像左侧的网页呢,而稍加努力就可以让它像右侧的那个网页一样?

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

本教程分为两个部分:首先,我们创建一个静态网页,了解 Flask 和 HTML 如何协同工作;接下来,我们将处理回调,以创建一个交互式页面。

本文的所有代码将存储在我的 GitHub 仓库中。我会在文章发布后的不久提供链接。

Bootstrap UI

使用 Bootstrap 创建吸引人的网页并不需要太多努力。包含 Bootstrap 5 文件并为 HTML 元素添加一些属性,可以轻松改善基本 HTML。

这不是一个 Bootstrap 教程,但让我快速向你展示上面两个网页中标题的基本 HTML 代码之间的区别。

基础 HTML

<header>
   <h1>Title</h1>
   <p>subtitle</p>
</header>

添加了 Bootstrap 属性

<header class="bg-primary text-white text-center py-2">
    <h1 class="display-4">Title</h1>
    <p class="lead">subtitle</p>
</header>

你可以看到,头部由两个元素组成,一个 <h1>,即顶级标题,以及一个段落 <p>。在 Bootstrap 版本中,这些元素有了额外的属性:头部本身具有 primary 背景色和白色文本,文本居中,上下边距设置为 2 像素;标题标签使用 display-4 字体,而段落的字体设置为 lead —— 这些字体在 Bootstrap 中定义,其中 display 字体大而粗体,而 lead 字体用于需要突出的普通文本。

这些特性是在 HTML class 属性中设置的。我们将在接下来的代码中看到更多这些特性,它们应该相当容易理解。我不会详细描述这些特性,但你可以在 他们的网站 上找到 Bootstrap 文档 —— 这将告诉你所需了解的一切。

一个 Flask 项目

Flask 框架使得创建基于 web 的应用程序变得简单。一个 Flask 应用通常由至少两个文件组成:一个 Python 应用程序和一个 HTML 模板。

Python 部分包含应用逻辑:例如,在数据可视化应用中,它可能将数据加载到 Pandas 数据框中,进行一些分析,并在 Plotly 中创建图表。HTML 模板定义了网页的布局,并由 Python 程序提供显示的数据。

/project_name
    |--- app.py
    |
    |--- /templates
            |
            |--- index.html

一个简单应用的目录结构,如我们将要创建的那样,应该类似于上面的示意图。主 Python 应用程序位于项目文件夹中,模板文件位于名为 templates 的子文件夹中。

当然,你还需要安装 Flask:

pip install flask

在这样做之前,你可能需要创建一个虚拟环境。

Flask 应用的 Python 部分定义了一个或多个应用将响应的路由。通常其中一个路由是‘/’,即项目的根目录。

下面是一个使用模板的最小 Python 应用。模板名为 index.html,必须位于 templates 文件夹中,并由 Flask 使用 render_template() 库函数呈现为网页。请注意,我们为网页标题创建了一个值并将其传递给函数。

from flask import Flask, render_template

app = Flask(__name__)

@app.route('/')
def index():
    title = "This is the title"
    return render_template('index.html', title=title)

下面是 index.html 模板,它期望 title 值被纳入其中。你可以看到,在 <h1> 标签中,标识符 title 被包围在双大括号中。

<!DOCTYPE html>
<html>
<body>
   <h1>{{title}}</h1>
</body>
</html>

Flask 使用 Jinja 模板引擎来用传递给 render_template() 的值替换 HTML 模板中的占位符。

你可以通过在终端中输入 flask run 来运行应用程序,你应该会得到类似于下面的响应。(这假设你将应用程序命名为 app.py —— 如果你用其他名字命名,需要输入 flask --app app_name run。将 app_name 更改为你的应用程序名称,不包含 .py 扩展名。)

flask run
 * Debug mode: off
WARNING: This is a development server. Do not use it in a production deployment. Use a production WSGI server instead.
 * Running on http://127.0.0.1:5000
Press CTRL+C to quit

将浏览器指向 http://127.0.0.1:5000localhost:5000,你将看到一个简单的网页,展示了在 Python 代码中定义的文本。

想要了解更多关于 Flask 的信息,你可以从查看他们的快速入门教程开始,但我希望我会在这里涵盖你开始所需的所有内容。

一个静态数据可视化应用程序

我们的第一个应用程序将基于迄今为止看到的内容,创建一个包含 Plotly 图表和一些辅助文本的网站。稍后,我们将添加一些交互功能。

让我们看看 Python 的部分。下面的代码列表中,请暂时关注以#### Simple template ####开头的部分。在这里,你可以看到我们定义了一个名为/simple的路由,这意味着我们通过将浏览器指向localhost:5000/simple来调用应用程序,而装饰器下方的simpleindex(),函数将会被执行。

在这个函数中,我们设置了一些文本和我们希望在网页上显示的图形。我们首先设置了一些变量,然后使用这些变量创建一个 HTML 模板将使用的参数字典。变量的名称清楚地表明了它们的使用方式。

get_graph()函数设置了图形参数。首先,它加载了 1881 年至 2022 年的全球温度异常数据,并追踪气候变化如何影响这段时间的温度(有关详细信息,请参见新数据表明 2023 年是有记录以来最热的夏天)。这些数据在下表中显示(在图表中会更清晰!)。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

全球温度异常

对于静态应用程序,我们将使用单列JJA,它指的是 6 月、7 月和 8 月的温度。交互式应用程序也会使用其他一些列,因此我们当前应用程序的period参数有一个默认值,之后可以被交互式应用程序更改。

数据用于创建一个 Plotly 条形图,生成的图形被转换为 JSON,网页将使用这些数据。因此,这些 JSON 数据被返回以设置字典参数中的graph条目。

回到之前的函数,我们现在需要使用我们设置的参数调用render_template。为了节省输入时间,我创建了一个名为template的辅助函数,它提取模板参数并将所有数据传递到网页上。

from flask import Flask, request, jsonify, render_template
import json
import pandas as pd
import plotly.express as px

app = Flask(__name__)

def get_graph(period = 'JJA'):
    df = pd.read_csv('GlobalTemps1880-2022.csv')
    fig = px.bar(df, x='Year', y = period, 
                 color=period, title = period, 
                 color_continuous_scale='reds', 
                 template='plotly_white', width=1000, height=500)

    graphJSON = fig.to_json()
    return json.dumps(graphJSON)

def template(params):
    return render_template(params['template'], params=params)

#### App ####

@app.route('/')
def index():
    return render_template('index.html')

#### Simple template ####

@app.route('/simple')
def simpleindex():
    header = "Global Temperature"
    subheader = "Global Temperature changes over the last few centuries"
    description = """The graph shows the increase in temperature year on year.
    The data spans the years 1881 to 2022 and includes temperature anomalies 
    for the months June through August.
    """
    params = {
        'template':'simpleindex.html',
        'title'   : header,
        'subtitle': subheader,
        'content' : description,
        'graph'   : get_graph()
    }
    return template(params)

#### Main ####

if __name__ == '__main__':
    app.run(debug=True)

现在来看 HTML 模板。下面的代码列表中,好消息是你可以忽略除<body>...</body>标签内的代码外的所有内容。其余的都是你需要的 Bootstrap 和 Plotly 网页的样板代码,可以剪切并粘贴到任何类似的网页中。最后的<script>标签也是样板代码,包含了 Bootstrap 的 Javascript 代码,可以安全忽略。

<!DOCTYPE html>
<html lang="en">
<head>
    <meta charset="UTF-8">
    <meta name="viewport" content="width=device-width, initial-scale=1.0">
    <link href="https://cdn.jsdelivr.net/npm/bootstrap@5.0.2/dist/css/bootstrap.min.css" rel="stylesheet" integrity="sha384-EVSTQN3/azprG1Anm3QDgpJLIm9Nao0Yz1ztcQTwFspd3yD65VohhpuuCOmLASjC" crossorigin="anonymous">
    <script src='https://cdn.plot.ly/plotly-latest.min.js'></script>
</head>
<body>

    <header class="bg-primary text-white text-center py-4">
        <h1 class="display-4">{{ params.title }}</h1>
        <p class="lead">{{ params.subtitle }}</p>
    </header>

    <div id = 'content' class="container mt-4">
        <div id='chart'></div>
        <div class="lead">{{params.content}}</div>
        </div>
    <script type='text/javascript'>
      var figure = JSON.parse({{params.graph | safe}})
      Plotly.newPlot('chart', figure, {});
    </script>

    <script src="https://cdn.jsdelivr.net/npm/bootstrap@5.0.2/dist/js/bootstrap.bundle.min.js" integrity="sha384-MrcW6ZMFYlzcLA8Nl+NtUVF0sA7MsXsP1UyJoMp4YLEuNSfAP+JcXn/tWtIaxVXM" crossorigin="anonymous"></script>

</body>
</html>

如果我们忽略使页面看起来漂亮的 Bootstrap 属性,我们得到的代码如下,这要简单得多,这就是我从这里开始引用的内容。

 <header>
        <h1>{{ params.title }}</h1>
        <p>{{ params.subtitle }}</p>
    </header>

    <div">
        <div id='chart'></div>
        <div>{{params.content}}</div>
    </div>
    <script type='text/javascript'>
      var figure = JSON.parse({{params.graph | safe}})
      Plotly.newPlot('chart', figure, {});
    </script>

我们之前见过使用 Jinja 参数,这次唯一的不同是我们将几个参数打包成一个名为 params 的字典。因此,我们通过在参数名称前加上字典的名称来引用它们。因此,<h1>{{params.title}}</h1> 只是将 title 参数放入一对标题标签中。其他三个文本参数的标签也是类似的。

为了绘制图表,我们需要一个可以放置图表的元素,并且该元素必须有一个 id 属性 (<div id='chart'></div>)。这个元素被放置在标题下方和描述上方。

下面的脚本元素再次是调用 Plotly Javascript 绘制图表的样板代码。唯一需要注意的是,当我们包含 graph 参数时,我们使用了 safe 关键字。这指示 Jinja 不要尝试解释 graph 中的任何特殊字符,而是将它们按字面意思对待。因此,代码如下:

var figure = JSON.parse({{params.graph | safe}})

现在,请记住模板必须在项目目录中的 templates 文件夹中,并且为了使此代码正常工作,数据文件必须位于项目目录本身(当然,你可以移动它,但在 Python 程序中打开时必须更改路径)。

所以运行应用程序并在浏览器中指向 localhost:5000/simple,你应该能看到一个如下图所示的网页。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

一个静态应用

这就是网页应用的静态版本。

一个互动数据可视化应用

但是为什么要限制在夏季?如果我们能够选择除 JJA 以外的其他时期不是很好吗?数据还包含全年的列,J-D,以及三个时间段:DJFMAMJJASON。(这些字母代表英文中的月份名称:十二月、一月、二月;三月、四月、五月;等等。)

为此,我们需要整合一个用户控件,用于选择适当的时间段。我选择了一个下拉菜单,用于显示各种时间段。它将与之前的网页非常相似(参见下图)。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

一个互动应用

代码最初也非常相似。主要区别出现在处理新图表的选择时。

当选择一个新值时,这将调用一个 Javascript 函数,该函数将值发送到服务器上的回调函数并等待响应。这个回调函数将返回一个新的图表,然后由调用的 Javascript 显示。

让我们先处理熟悉的内容。下面是实现新端点 /ddsimple 的函数。

@app.route('/ddsimple')
def ddsimpleindex():
    # The root endpoint builds the page
    header = "Global Temperature"
    subheader = "Global Temperature changes over the last few centuries"
    description = """The graph shows the increase in temperature year on year.
    The data spans the years 1881 to 2022 and includes temperature anomalies for periods of each year as indicated.
    """
    menu_label = "Select a period"

    params = {    
        'template': 'ddsimpleindex.html',
        'title': header,
        'subtitle': subheader,
        'content' : description,
        'menu_label': menu_label,
        'options' : [{'code':'J-D', 'desc':'Whole year'},
                     {'code':'DJF','desc':'Winter (North)'},
                     {'code':'MAM','desc':'Spring (North)'},
                     {'code':'JJA','desc':'Summer (North)'},
                     {'code':'SON','desc':'Autumn/Fall (North)'}],
        'graph'   : get_graph()
    }
    return template(params)

如你所见,它与/simple端点非常相似。区别(除了名称)在于额外的参数:菜单的标签和一个表示菜单项列表的字典,首先是一个对应于数据框列的值,其次是该值的文本描述,这些描述将显示在菜单中。

HTML 稍有不同,因为它展示了 Jinja 的更复杂用法。代码如下。

<form id="userForm" name="form1" onChange="getFormValues('form1')">
    <div class="mb-3">
        <label for="dropdown" class="form-label lead">{{params.menu_label}}</label>
        <select class="form-select" id="dropdown" name="dropdown">
            {% for opt in params.options %}
                <option value="{{opt.code}}">{{opt.desc}}</option>
            {% endfor %}
        </select>
    </div>
</form>
<div>
    <!-- Main Content Area -->
    <p class="lead">{{params.content}}</p>
</div>
<div id="graph"></div>

在这里,我们构建了一个包含下拉菜单的表单。该表单还包含一个菜单标签,这个标签包含在我们之前见过的双花括号中。

与之前的示例相比,主要区别在于菜单的构建方式。在<select>元素内,我们需要放置一系列代表菜单选项的<option>标签。<option>标签具有一个值和一个描述,这些值和描述是在params.options字典中定义的。我们通过执行类似{% for opt in param.options %}的 Jinja 循环来包含这些值和描述,该循环遍历字典,将每个元素放入本地变量opt中。然后,我们利用这些值和描述,通过使用opt.codeopt.desc将它们插入到<option>标签中。

你可以在Flask 教程:Pythonbasics.org 网站上的模板中找到 Jinja 模板的简单示例和解释。

表单标签中还有另一个对我们目的至关重要的部分。表单标签内有一个名为onChange的属性,它接受某种类型的动作值,在这种情况下,它是一个 Javascript 函数,每当表单中的值发生变化时——在本例中,当选择菜单中的选项时,就会调用这个函数。

这就是有趣的开始。

回调

为了用新图表更新网页,我们使用回调机制,下面的图表展示了浏览器与服务器之间的事务。请注意,响应请求新图表的方式是不是重新加载页面,而是更新页面——这种方式更快,并且避免了重新加载时出现的瞬时空白屏幕,从而提供了更好的用户体验。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

使用回调更新网页

回调是通过网页上表单的变化来调用的,正如我们之前提到的。这种机制是一个在表单的onChange属性中标识的 Javascript 函数。

我将详细解释 Javascript 函数的工作原理,但基本上,它从表单中获取值,并将这些值发送到 Flask 应用中的回调端点。

现在,如果想到编写 Javascript 让你感到不安,不用担心,你不需要真正了解这些内容,你可以直接复制它,它适用于你可能想在网页上包含的任何表单。因此,对冒险者而言,解释如下,对于其他人——直接跳到下一部分。

实际上有两个函数:从表单中获取值的函数如下所示。

 function getFormValues(f) {
            const form = document.forms.namedItem(f);
            const formData = new FormData(form);
            const value = Object.fromEntries(formData.entries());
            postJSON(value);
        }

所有代码都使用内置的 Javascript 函数:第一行从文档(即网页)中获取表单;第二行检索包含表单数据的数据结构;第三行提取该结构中的所有值。

最后,这些值被传递给另一个函数 postJson

 async function postJSON(data) {
            try {
                const response = await fetch("/callback", {
                    method: "POST",
                    headers: {
                        "Content-Type": "application/json",
                    },
                    body: JSON.stringify(data),
                });

                const result = await response.json();
                console.log("Success:");//, result);

                drawGraph(result);
            }
            catch (error) {
                console.error("Error:", error);
            }
        }

这个函数是一个异步函数,这意味着在调用后,程序执行会立即返回到调用代码,而异步函数会在一个单独的执行线程中继续运行。也就是说,它会与网页的执行并行地继续执行所需的操作。

postJSON 函数接收需要发送到 Python 回调代码的数据,并使用异步 fetch 函数发送。fetch 接收数据应该传递到的端点和数据本身作为参数——我们使用 HTTP POST 机制来发送数据。postJSON 等待 fetch 完成,即等待服务器返回一些数据。然后将这些数据传递给 drawGraph 函数,更新网页上的图表。

请注意,代码被包含在 try... catch... 块中。这与在 Python 程序中看到的基本相同:如果 try 块中的代码失败——没有响应或其他通信故障——那么该故障将被记录到控制台中。

Python 回调

为了使这一切工作,我们需要一个 Flask 代码中的回调函数,该函数将接收数据,对其进行处理(即创建新图表),并返回结果。如下所示:

@app.route('/callback', methods=['POST'])
def callback():
    # The callback updates the page
    if request.is_json:
        data = request.get_json() 
        return get_graph(period=data['dropdown'])
    else:
        return jsonify({"error": "Invalid JSON data"}), 400

首先需要定义回调的端点,你可以看到我们还指定了端点期望使用 POST 方法发送数据。

我们也期望数据为 JSON 格式,如果不是,则返回错误。

如果数据 JSON,我们从下拉菜单中提取值,并将其传递给get_graph函数,该函数绘制图表并以 Plotly 期望的 JSON 格式返回图表数据。这个图表数据由网页上的 Javascript 函数接收,页面得到更新。

这段代码将为你提供上面显示的交互式网页。

Flapjax —— 名字中有什么

诚然,这个名字有些刻意:它代表了Flask,Python,Javascript 和 ax,这些代表了使这种技术得以实现的异步通信。

我希望你能看到,通过使用预先编写的模板和一段模板代码,你可以创建有用的交互式网页,同时主要集中于 Python 代码的逻辑。

这个应用程序仅包含一个图表,当用户选择新选项时,该图表会更新,但可以通过这种方式从 HTML 表单中收集任意数量的值,这些值可以由 Flask 应用程序处理,然后用新信息更新网页(也许我们会在未来的文章中探讨这个问题)。

所有在这里展示的代码和数据可以从我的 GitHub 仓库 下载(查看 jinja-article 文件夹)。

更新:我在 GitHub 仓库的 reuse 文件夹中编写了一个新应用程序,该应用程序使用一组新的数据——HTML 和 Javascript 保持不变,只有 Python 代码和数据发生了变化: 重用 Flapjax 模板和代码

我希望你觉得这些内容有用。如果你想查看更多我的作品,请访问我的 网站,并通过订阅我的免费通讯来获取我发布的更新,点击这里:

[## 数据可视化、数据科学和 Python | Alan Jones | Substack

关于数据科学、数据可视化及主要使用 Python 进行的动手编码的教程和其他文章。点击…

technofile.substack.com

如果 Flask 不是你的首选,我根据我在 Medium 上的文章编写了一本电子书 从头开始学习 Streamlit

说明和参考文献

本文和应用程序中使用的数据来源于下面第 1 和第 2 条说明中的描述。

  1. GISTEMP 团队,2023: GISS 地表温度分析 (GISTEMP),第 4 版。NASA 戈达德空间研究所。数据集访问日期 2023–09–19,网址:data.giss.nasa.gov/gistemp/。请注意,NASA 数据集的使用没有特定的许可证。NASA 将这些数据集免费提供用于非商业目的,但应给予归属(如上所述)。

  2. Lenssen, N., G. Schmidt, J. Hansen, M. Menne, A. Persin, R. Ruedy, 和 D. Zyss, 2019: GISTEMP 不确定性模型的改进。J. Geophys. Res. Atmos., 124, no. 12, 6307–6326, doi:10.1029/2018JD029522。

除非另有说明,所有的图像、图表、截图和代码均由作者创建。

关注 TDS 列表,发现我们的最佳文章

原文:towardsdatascience.com/follow-tds-lists-to-discover-our-best-articles-760fced78c0f?source=collection_archive---------18-----------------------#2023-07-18

一种方便的方式来跟进我们作者的最强和最受欢迎的作品

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传 TDS Editors

·

关注 发表在 Towards Data Science ·2 分钟阅读·2023 年 7 月 18 日

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

照片由Alisa Anton拍摄,发布在Unsplash

我们每周发布数十篇精心制作的文章,主题涵盖最新的机器学习工具和 Python 库到前沿研究。现实中,阅读我们分享的每一篇文章非常困难,但这并不意味着你应该错过你最有可能喜欢的那些文章。

为了帮助你跟上数据科学、人工智能和机器学习的最新动态,我们创建了一些Medium 列表,在这些列表中汇集了我们推荐的阅读内容。我们将定期更新每个列表,因此你只需保存一次,就可以轻松访问你 Medium 图书馆的收藏列表标签

(不确定如何保存列表?只需点击以下任意列表上的查看列表按钮,然后点击列表标题下的书签符号。)

寻找 TDS 的最佳内容?我们创建了数据科学和机器学习必读列表,让你随时能接触到我们各个话题和格式的最强文章。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

TDS 编辑

数据科学和机器学习必读

查看列表275 篇故事外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

如果你喜欢较长、深思熟虑的文章,我们的数据科学和机器学习深度解析列表就是你最好的去处——从教程到解释,这些帖子因其深入和易懂而脱颖而出。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

TDS 编辑

数据科学和机器学习深度解析

查看列表248 篇故事外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

坚信群体智慧的人应该收藏TDS 热门列表,我们将在此突出展示流行和广泛传播的故事。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

TDS 编辑

TDS 热门

查看列表197 篇故事外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

我们希望您喜欢使用我们的新列表,并且这些列表能激发您去探索那些您之前没有关注过的作者的作品。如果您有其他您认为有用的 TDS 列表的想法,请留言告知我们。

祝您阅读愉快!

遵循此数据验证过程以提高数据科学准确性

原文:towardsdatascience.com/follow-this-data-validation-process-to-improve-your-data-science-accuracy-99422dfbee72

当训练数据和推断数据来自不同来源时

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传 Matt Przybyla

·发表于 Towards Data Science ·阅读时间 8 分钟·2023 年 9 月 1 日

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片由 NordWood Themes 提供,出处:Unsplash [1]。

目录

  1. 介绍

  2. 启用数据收集

  3. 设置基准

  4. 检测异常值

  5. 摘要

  6. 参考资料

介绍

本文旨在为数据科学家提供数据验证过程的指南,特别是那些刚开始或希望改进现有数据验证流程的人,提供了一些示例作为一般性概述。首先,我想在这里定义数据验证,因为它对其他类似职位的意义可能不同。为了本文的目的,我们将数据验证定义为确保用于模型的训练数据与推断数据一致的过程。对于一些公司和一些用例,如果数据来自同一来源,则无需担心此问题。因此,这个过程必须发生,并且只有在数据来自不同来源时才有用。数据不会来自相同来源的一些原因包括,如果你的训练数据是历史数据和自定义数据(例如:从现有数据中提取的特征),和/或你的推断数据来自实时表格,而训练数据是快照数据。总之,这种不匹配的原因有很多,制定一个规模化的过程来确保你在推断时提供给模型的数据是你——即训练好的模型数据所期望的,将极有益处。

启用数据收集

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

照片由 Dennis Kummer 提供,来自 Unsplash [2]。

有很多方法可以启用数据收集。但再次强调,我们首先要定义收集的数据,即推理数据。我们预计我们的训练数据(包括训练集和测试集)已经存储在某处,也许是在 S3、文件存储工具、数据库中的临时表,甚至是 CSV 文件等。

然而,推理数据在预测时使用更为常见,但不一定会被存储。因此,你需要启用其收集。类似于训练数据的示例,你可以将推理数据存储在数据库表、数据湖、CSV 文件、JSON 文件、S3、已存储的 Kafka 流等。选择哪种收集方式将取决于你公司提供的选择。

还有一些工具、库和 Python 函数可以启用推理数据存储。

下面是一些可以用来收集推理数据的工具:

  • Kafka 流入数据库:调用模型的服务端点可以将数据加载到 Kafka 流中,然后从那里,你可以将数据加载到数据库表中

  • Amazon Web Services 捕获内容 [3]:这个工具非常用户友好,并允许你在模型代码中添加一个参数。你将需要使用 Amazon SageMaker 来使用这个工具。

代码示例如下,你可以看到如何简单地将参数添加到 AWS SageMaker 模型代码中,如果这是你正在使用的工具:

# Configuration object passed in when deploying Models to SM endpoints
data_capture_config = DataCaptureConfig(
    enable_capture = enable_capture, 
    sampling_percentage = sampling_percentage, # Optional
    destination_s3_uri = s3_capture_upload_path, # Optional
    capture_options = ["REQUEST", "RESPONSE"],
)

# Example code from AWS [3]

对于具体工具的详细介绍,我将在下一篇文章中展开,并重点介绍代码。目前,在你的数据验证计划中,你需要确保有一种方法来捕获和存储推理数据,以及你的训练数据。

设定基准线

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

变量示例。截图由作者提供 [4]。

基准线有助于理解你的训练数据。理解你的数据应该如何在视觉上呈现并记住它可能会让人感到困惑,因此最好创建一个从数据中导出描述性统计的自动化过程。

基准数据示例

对于分类列字符串、对象、类别等),你可以为已知数据和预期缺失的数据创建基准线:

  • 唯一的 X 列值,例如:“八月、九月

对于数值列整数、浮点数):

  • 列 X 平均值

  • 列 Y 标准差

  • 列 X、Y 最小值、最大值等,例如:20,000.60

尽管这看起来相当简单,但当你的模型和数据变得更大时,尽早跟踪数据将是有益的,这样你可以理解何时以及为何你的推断数据中会出现异常值。

上述是创建基线统计的一般示例,下面将介绍一个非常易于使用的具体工具

Pandas Profiling:

这个 Python 库 [5] 为你执行了大量的探索性工作,这将作为你的数据基线信息。总体而言,这个库会查看以下类型的基线统计。

分位数

  • 最小值

  • 第 5 百分位数

  • Q1,中位数,Q3,第 95 百分位数

  • 最大值

  • 范围

  • 四分位距(IQR)

描述性

  • 标准差

  • 变异系数(CV)

  • 峰度

  • 平均值,中位数绝对偏差(MAD)

  • 偏度

  • 总和

  • 方差

  • 单调性

# install library 
#!pip install pandas_profiling
import pandas_profiling
import pandas as pd
import numpy as np

# create data 
df = pd.DataFrame(np.random.randint(0,200,size=(15, 6)), columns=list('ABCDEF'))

# run your report!
df.profile_report()

检测异常值

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

Will MyersUnsplash 拍摄的照片 [6]。

现在你已经收集了训练数据和推断数据,并且有了预期的基线值,你可以知道当推断中的值变为异常值时。

需要记住的是,训练数据本身被假设为 100%准确,因为它被用作真相源。

有了收集的推断数据,你可以使用大量工具、库、模块等,将推断数据与训练数据进行比较。

为了进一步说明这一点,我会给出两个例子,一个没有问题,另一个是异常值。

这里是一个分类数据和数值数据的非异常值示例:

  • 训练数据分类特征(months)→ “August

  • 推断数据分类特征(months)→ “September

  • 训练数据数值特征(temperature)→ 100

  • 推断数据数值特征(temperature)→ 110

这里是一个分类数据和数值数据的异常值示例:

  • 训练数据分类特征(months)→ “August

  • 推断数据分类特征(months)→ “AUgust

  • 训练数据数值特征(temperature)→ 100

  • 推断数据数值特征(temperature)→ 1000000000

再次看似显而易见,但想象一下你有 50+个特征和 100,000 行数据。几乎不可能检查所有数据是否符合预期。

此外,为了更突出这个例子,虽然“August”与“AUgust”看起来相似,但模型当然不会将它们视为相同的值。对于数值数据,我们也不希望我们的温度读数如此之高。

虽然这些值在推理时显然不正确,但模型仍然可以对它们进行预测(取决于你的模型),这使得知道何时发生这种情况变得更加困难。

作为额外内容,你需要设置一个过程以确保你的预测与你的期望一致。

举个具体的例子,你可以使用Amazon SageMaker Model Monitoring工具,或者如果你在其他平台,可以使用像TensorFlow Data Validation [7]这样的开源库。

 stats = tfdv.generate_statistics_from_tfrecord(data_location=path)

tfdv.visualize_statistics(stats)

other_stats = tfdv.generate_statistics_from_tfrecord(data_location=other_path)

anomalies = tfdv.validate_statistics(statistics=other_stats, schema=schema)

# You cannot run this code without importing the libraries, however, 
# it is the outline for using the anomlay detection or outlier detection. 

# Code from AWS [7]

以下是你为什么要使用这个特定库的主要用例[7],而且正如我之前所说,AWS 在这里很好地解释了原因:

  • 验证新数据以确保我们没有突然开始接收到不良特征

总结

了解何时发生离群值可能非常重要,因为你的模型可能会预测出由于数据不准确而可能造成伤害的结果。你不希望在涉及安全的欺诈检测中出现这种问题。

这篇文章作为一个快速大纲和提醒,以确保你的数据不仅对训练模型有用,也适用于预测时。

总结一下,我们在这篇文章中讨论了以下内容:

* Enabling Data Collection
* Setting a Baseline
* Detecting Outliers

我希望你觉得我的文章既有趣又有用。如果你的数据验证经验相同或不同,请随时在下方评论。为什么会这样?还有哪些其他话题你认为应该更多讨论?这些确实可以进一步澄清,但我希望我能对数据科学家数据验证的常见大纲提供一些见解。如果你希望看到关于这个过程中特定工具的文章,也可以随时评论。

我与这些公司没有任何关联。

请随时查看我的个人资料, Matt Przybyla以及其他文章,并通过以下链接订阅以接收我的博客的电子邮件通知,或者通过点击屏幕顶部的订阅图标**,如有任何问题或评论,请在 LinkedIn 上联系我。*

感谢阅读!

订阅链接: datascience2.medium.com/subscribe

推荐链接: datascience2.medium.com/membership

(如果你在 Medium 上注册会员,我会获得佣金)

参考文献

[1] 图片由 NordWood Themes 提供,在 Unsplash 上,(2018)

[2] Matt Przybyla,变量示例的截图,(2020)

[3] 2023, 亚马逊网络服务公司或其附属公司。版权所有,从实时端点捕获数据,(2023)

[4] 由 Stephen Dawson 拍摄于 Unsplash,(2018)

[5] pandas-profiling,GitHub 上的文档和所有贡献者,(2020)

[6] 由 Will Myers 拍摄于 Unsplash,(2018)

[7] TensorFlow,除非另有说明,本页面内容按 知识共享署名 4.0 国际许可证 许可,代码示例按 Apache 2.0 许可证 许可。有关详细信息,请参阅 Google 开发者网站政策。Java 是 Oracle 及/或其附属公司的注册商标,TensorFlow 数据验证,(2023)

预测多个视野:以天气数据为例

原文:towardsdatascience.com/forecast-multiple-horizons-an-example-with-weather-data-8d5fa4321e07

使用预测视野作为特征预测瑞士的降水量。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传 Davide Burba

·发表于 Towards Data Science ·阅读时长 8 分钟·2023 年 8 月 6 日

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

天气预报,作者:Giulia Roggia。已获许可使用。

  • 介绍

  • 示例:瑞士的降水量

介绍

传统方法

当我们想要预测时间序列的未来值时,我们通常对多个未来视野感兴趣,例如 1、2 或 3 个月后的情况。预测这些不同视野的传统方法是为每个目标视野训练一个单独的模型。

常见替代方案

一种常见的替代方法是训练一个短期视野的单一模型,然后通过递归应用(即将之前的预测作为输入来产生后续预测)将其扩展到多视野。然而,这种方法可能在生产系统中实施复杂,并且可能导致误差传播:在近视野上的误差可能对后续的视野产生不利影响。

另一种替代方法是用多变量模型同时预测所有视野。然而,支持多变量输出的模型种类有限,并且在数据处理和模型维护上需要额外的努力。

视野作为特征

一种更简单的方法是将为每个视野准备的数据拼接在一起,并添加一个新的“视野”特征。这种方法具有几个优点

  1. 它易于理解和实施,因为它只需要训练和维护一个模型。

  2. 它可能提高预测准确性,因为模型是在更大的数据集上训练的。它甚至可以作为一种“数据增强”技术使用:如果你只对几个视野感兴趣,你仍然可以在训练阶段添加额外的视野来改善模型估计。

  3. 该模型可以用于预测其未经过训练的范围,这在你有多个范围需要预测时可能很有帮助。

这种方法是 全球模型 的另一种表现形式,但在多个预测范围的背景下,而不是多个时间序列。因此,它有类似的 缺点

  1. 如果特定预测范围的性能开始下降,更新模型很难而且可能会影响其他范围的预测。

  2. 当出现新的预测范围时,可能需要完全重新训练(尽管这不是强制性的)。

  3. 你不能使用传统的预测模型(如 ARIMA、指数平滑等)。

示例:瑞士的降水量

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

卢加诺湖,图片由 Makalu 提供,来自 Pixabay

在本节中,我们展示了“范围作为特征”技术的具体 Python 实现,并将其与每个范围训练一个模型的传统方法进行比较。

我们的目标是预测卢加诺的降水量,这是一座位于瑞士的城市。此示例中使用的数据由 MeteoSwiss 提供,数据可以在 这里 获得(数据使用已获许可)。

数据加载

让我们首先导入必要的库,以便处理数据、可视化数据和训练 LightGBM 模型(如果你对模型的选择感兴趣,可以查看这篇文章)。

import pandas as pd
import plotly.graph_objects as go
from lightgbm import LGBMRegressor

让我们将数据加载到 Pandas 数据框中:

def load_climate_data(path):
    # Load data in a dataframe.
    with open(path,"r", encoding="ISO-8859-1") as f:
        data = f.readlines()
    columns = data[27].split()
    data = [v.split() for v in data[28:]]
    data = pd.DataFrame(data, columns = columns)

    # Fix time
    data["time"] = pd.to_datetime(data.Year + "-" + data.Month.astype(str))
    data = data.drop(columns = ["Year","Month"])
    data = data.set_index("time")

    # Fix types
    data = data.replace("NA",None)
    return data.astype(float)

DATA_PATH = "./data/climate-reports-tables-homogenized_LUG.txt"
data = load_climate_data(DATA_PATH)

快速数据探索

现在我们可以查看数据:

def show_data(data,title=""):
    trace = [go.Scatter(x=data.index,y=data[c],name=c) for c in data.columns]
    go.Figure(trace,layout=dict(title=title)).show()

# Let's visualize the data.
show_data(data,"Weather Data in Lugano")

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

卢加诺的温度和降水量。图片来源:作者。

我们可以看到,从 1864 年开始,我们有每月的温度和降水量数据(这几乎是 160 年!)。让我们利用 Plotly 的交互功能更仔细地查看。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

过去 10 年卢加诺的温度。图片来源:作者。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

过去 10 年卢加诺的降水量。图片来源:作者。

我们可以看到数据具有季节性,温度比降水量更规律,这符合预期。

由于降水量更难预测,它们也更有趣!让我们尝试预测未来 1、2 和 3 个月的值。

数据工程

让我们准备数据以训练预测模型。我们将使用滞后的降水量和温度值来预测未来的降水量:

  • 对于降水量,我们考虑最近 3 个滞后值,以及目标月份之前 1 年、2 年和 3 年的滞后值。

  • 对于温度,我们考虑最近 2 个滞后值,以及目标月份之前 1 年和 2 年的滞后值。

我们还添加了一个“month_horizon”特征,指示我们正在预测哪个月份(例如 1 月或 4 月),以帮助模型学习季节性。注意这与“作为特征的时间跨度”不同,后者取值为 1、2 或 3。

def build_features(data, horizon):
    """Build lagged features.

    We depend on horizon due to relative lags shift. 
    E.g, if the horizon is equal to 1, the target value 
    of 12 months before corresponds to a lag of 11.
    """
    # Here we hardcode values for simplicity, but everything could 
    # (and should) be parametrized.
    precipitation_lags = [0, 1, 2, 12 - horizon, 24 - horizon, 36 - horizon]
    temperature_lags = [0, 1, 12 - horizon, 24 - horizon]

    # Concatenate precipitation and temperature features.
    features = pd.concat(
        [
            build_lagged_features(data.Precipitation, lags=precipitation_lags),
            build_lagged_features(data.Temperature, lags=temperature_lags),
        ],
        axis=1,
    )

    # Add horizon_month as a feature.
    features["horizon_month"] = (features.index.month + horizon - 1) % 12 + 1

    # Trick to later allow concatenation of features for different 
    # target horizons.
    features = features.rename(
        columns={
            f"Precipitation_lag_{12-horizon}": "Precipitation_lag_12_before_target",
            f"Precipitation_lag_{24-horizon}": "Precipitation_lag_24_before_target",
            f"Precipitation_lag_{36-horizon}": "Precipitation_lag_36_before_target",
            f"Temperature_lag_{12-horizon}": "Temperature_lag_12_before_target",
            f"Temperature_lag_{24-horizon}": "Temperature_lag_24_before_target",
        }
    )

    return features

def build_lagged_features(series, lags):
    return pd.concat(
        [series.shift(lag).rename(f"{series.name}_lag_{lag}") for lag in lags],
        axis=1,
    )

我们现在可以为每个时间跨度构建目标和特征。

def build_target_features(data, horizon):
    targ = build_target(data.Precipitation, horizon)
    feat = build_features(data, horizon)

    # Drop missing values generated by lags/horizon.
    idx = ~(feat.isnull().any(axis=1) | targ.isnull())
    feat = feat.loc[idx]
    targ = targ.loc[idx]

    return targ, feat

def build_target(series, horizon):
    return series.shift(-horizon)

# Let's build the targets and features for each horizon.
HORIZONS = [1,2,3]
target_features = {h: build_target_features(data, h) for h in HORIZONS}

拆分训练集和测试集

我们将数据中最近 10 年的部分作为测试集:

def split_train_test(target_features, test_size):
    targ_feat_split = {}
    for horizon, (targ,feat) in target_features.items():
        targ_train = targ.iloc[:-test_size]
        feat_train = feat.iloc[:-test_size]
        targ_test = targ.iloc[-test_size:]
        feat_test = feat.iloc[-test_size:]

        targ_feat_split[horizon] = targ_train, feat_train, targ_test, feat_test

    return targ_feat_split

TEST_SIZE = 10 * 12
targ_feat_split = split_train_test(target_features, test_size=TEST_SIZE)

模型训练

由于我们只关心比较这两种不同的方法,因此两种情况下我们都保持默认超参数。

让我们开始传统方式训练模型,为每个时间跨度训练一个:

def train_models_by_horizon(targ_feat_split, model_params=None):
    if model_params is None:
        model_params = {}

    # Train one model for each horizon
    models_by_horizon = {}
    for horizon, (targ_train,feat_train,_,_) in targ_feat_split.items():
        model = LGBMRegressor(**model_params)
        model.fit(feat_train, targ_train)
        models_by_horizon[horizon] = model

    return models_by_horizon

models_by_horizon = train_models_by_horizon(targ_feat_split)

现在让我们运行另一种方法:使用时间跨度作为特征的单一模型:

def train_model_across_horizons(targ_feat_split, model_params=None):
    if model_params is None:
        model_params = {}

    # Concatenate data across horizons.
    targ_train_all = []
    feat_train_all = []
    for horizon, (targ_train,feat_train,_,_) in targ_feat_split.items():
        # Add horizon as a feature.
        feat_train = feat_train.copy()
        feat_train["target_horizon"] = horizon

        targ_train_all.append(targ_train)
        feat_train_all.append(feat_train)

    targ_train_all = pd.concat(targ_train_all)
    feat_train_all = pd.concat(feat_train_all)

    # Train a single model.
    model = LGBMRegressor(**model_params)
    model.fit(feat_train_all, targ_train_all)

    return model

model_shared = train_model_across_horizons(targ_feat_split)

在测试集上进行预测

让我们用这两种方法对测试集进行预测:

def predict_models_by_horizon(targ_feat_split, models_by_horizon):
    preds = {}
    for horizon, (_,_,_,feat_test) in targ_feat_split.items():
        preds[horizon] = models_by_horizon[horizon].predict(feat_test)
    return preds

def predict_model_across_horizons(targ_feat_split, model):
    preds = {}
    for horizon, (_,_,_,feat_test) in targ_feat_split.items():
        # Add horizon as a feature.
        feat_test = feat_test.copy()
        feat_test["target_horizon"] = horizon

        preds[horizon] = model.predict(feat_test)
    return preds

preds_by_horizon = predict_models_by_horizon(targ_feat_split, models_by_horizon)
preds_model_shared = predict_model_across_horizons(targ_feat_split, model_shared)

错误分析

现在让我们评估预测性能。首先,我们将输出合并成便于比较的格式。

# Let's combine the output in a convenient format.
output = {}
for horizon in HORIZONS:
    df = targ_feat_split[horizon][2].rename("target").to_frame()
    df["pred_model_by_horizon"] = preds_by_horizon[horizon]
    df["pred_model_shared"] = preds_model_shared[horizon]
    output[horizon] = df

下面您可以看到时间跨度为 1 个月的输出数据框的开头部分。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

时间跨度为 1 个月的输出数据框。图片由作者提供。

我们现在计算并打印整体和每个时间跨度的平均绝对误差(MAE)。

def print_stats(output):
    output_all = pd.concat(output.values())
    mae_by_horizon = (output_all.target - output_all.pred_model_by_horizon).abs().mean()
    mae_shared = (output_all.target - output_all.pred_model_shared).abs().mean()
    print("                 BY HORIZON     SHARED")
    print(f"MAE overall    :    {mae_by_horizon:.1f}         {mae_shared:.1f}\n")

    for h,df in output.items():   
        mae_by_horizon = (df.target - df.pred_model_by_horizon).abs().mean()
        mae_shared = (df.target - df.pred_model_shared).abs().mean()
        print(f"MAE - horizon {h}:    {mae_by_horizon:.1f}         {mae_shared:.1f}")

# Let's show some statistics.
print_stats(output)

这将产生:

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

我们看到跨时间跨度共享的模型总是导致较低的误差。这在某种程度上是预期的,因为它是在更大的数据集上训练的。

让我们看看一些预测结果:

# Let's have a look at the predictions.
for horizon, df in output.items():
    show_data(df,f"Predictions at Horizon {horizon}")

这里为了简洁,只显示时间跨度为 3 个月的预测:

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

时间跨度为 3 个月的预测。图片由作者提供。

我们看到模型仍然无法捕捉到极端事件(如 2014 年 8 月)。这并不令人惊讶,因为:

  • 我们对系统的实际状态了解非常有限。

  • 极端事件很难通过机器学习进行预测,因为定义上在训练集中观察到的实例有限。

下一步

我们可以尝试几个方法来提高预测性能。以下是一些非详尽的列表:

  • 通过使用 L2 损失来惩罚大误差。

  • 目标转换,例如取平方根。

  • 包括月内数据。包括额外的数据来源。

  • 使用在更多目标上训练的全局模型。

  • 更频繁地进行再训练,例如每年。

  • 超参数调整。

要评估这些变化是否提高了预测准确性,您应该依赖于一个可靠的回测策略

本示例中使用的完整代码可以在 这里 找到

喜欢这篇文章吗? 查看我的其他文章 并关注我以获取更多内容! 点击这里 阅读无限制的文章,并在不增加您额外成本的情况下支持我 ❤️

像大师一样预测多个时间序列

原文:towardsdatascience.com/forecast-multiple-time-series-like-a-master-1579a2b6f18d

从局部到全局算法

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传 Bartosz Szabłowski

·发表于Towards Data Science ·阅读时间 30 分钟·2023 年 4 月 26 日

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

图片来源:Jesús RochaUnsplash

我在商业中处理多个时间序列的预测(准确来说,是需求预测)。在我之前的文章中Sell Out Sell In Forecasting,我介绍了我在雀巢实施的需求预测方法。在这篇文章中,我想向你介绍目前用于预测多个时间序列的通用(这并不意味着理想)算法——例如最先进的时间序列算法。对于零售商或制造商来说,预测需求对业务至关重要。它允许他们制定更准确的生产计划并优化库存。不幸的是,许多公司(并非雀巢 😃 )并未意识到这个问题,他们仍然使用简单统计的电子表格。如果他们能改变这种情况,他们可以显著降低成本。毕竟,仓储和过时产品——这也是额外的成本。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

如何预测多个时间序列,作者提供的图片

很难找到一个数据科学领域的人不熟悉Scikit-learn。对于数据框架,你可以使用Scikit-learn来完成机器学习中涉及的大部分元素——从预处理到超参数选择、评估和模型预测。我们可以将线性回归、决策树或支持向量机(SVM)分配给变量model,并每次使用相同的方法,如fitpredict。我们有很大的灵活性,但也有一种简单的方式来实现解决方案。

对于时间序列来说,情况是不同的。如果你在实验并想比较不同的算法,算法本身不仅仅是一个问题。

如果你开始处理时间序列,你需要对它们进行处理,例如重新采样或填补缺失值——Pandas对此非常有用。

如果你想进行分解、可视化 ACF/PACF,或检查平稳性测试,那么Statsmodels库将非常有用。

对于可视化,你可能会使用Matplotlib,即使不是这个库,也有许多建立在它之上的库。

当你想使用不同的算法时,乐趣才开始。当你想使用ARIMA时,你可能会使用pmdarimaProphet是另一个库。典型的Machine Learning 算法可以在之前提到的Scikit-learn中找到,但你也可能想使用像LightGBMCatBoost这样的提升模型。对于Deep Neural Networks 和最新论文中的架构,PyTorch Forecasting值得使用。

WOW🤯 你可能需要的库非常多。如果你想能够使用上述提到的库,这将是一个大量的工作,因为大多数库使用不同的 API、数据类型,并且对于每个库中的模型,你都必须准备自己的回测和超参数选择函数。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

Library Darts,由作者提供的图片,灵感来自于图书馆文档

这里我们得到帮助的是Darts,它试图成为时间序列的Scikit-learn,其目的是简化时间序列的工作。它的功能通常基于其他库,例如,它使用Statsmodels进行分解。如果某些功能在其他库中没有实现,Darts也能很好地与其他库协作,你可以将其与MatplotlibSeaborn互相配合使用进行比较。

在必要的地方,它们有自己的实现,但它们不想重新发明轮子,而是使用其他流行时间序列库中已经存在的东西。

这一概念在时间序列领域并不新鲜,还有其他很好的库,例如sktimeGluonTSnixtla,但在我看来,Darts的入门门槛最低,功能也更为完善。这不是对这个库的广告,归根结底,你的预测应该为你工作的企业带来价值。你也可以从头开始用代码编写这些模型。我将使用 Darts 进行以下示例,但你也可以在上述提到的库中找到这些模型(全部或部分)。如果我们想训练多个本地模型,我认为 Darts 库在优化计算方面还有改进的空间——可以尝试nixtla库,它提供与 Spark、Dask 和 Ray 的兼容性。

从我的角度来看,Darts 已经是一个成熟的库,并且仍在不断开发中,只需查看变更日志即可。现在你可以按照标准方式进行安装:

pip install darts

一旦我们在环境中安装了库,就可以导入并在实践中使用它。

import darts

单变量与多变量时间序列

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

单变量与多变量时间序列,作者提供的图像

如上图所示,基于Walmart 数据集,你可以看到单变量多变量时间序列。现在,许多问题涉及同时处理多个点。这些数据可以来自各种过程,可以是这个例子和我日常工作的需求预测,也可以是能源消耗预测、公司股市收盘价、从车站租用的自行车数量等等许多其他问题。

除了时间序列本身,我们还可能有其他变量,其中一些可能已知未来值,而其他仅在过去可用——稍后会详细讲解。

在这篇文章中,我想向你展示预测多个时间序列的不同方法,但我希望它是实用的——以便你不仅仅停留在理论层面。所以让我们导入所有后续使用的库——包括 Darts 和其他数据科学家熟知的库。

# multiprocessing
from joblib import Parallel, delayed

# data manipulation
import numpy as np
import pandas as pd
from darts import TimeSeries
from darts.utils.timeseries_generation import datetime_attribute_timeseries

# data visualization
import matplotlib.pyplot as plt
import seaborn as sns
from tqdm import tqdm

# transformers and preprocessing
from darts.dataprocessing.transformers import Scaler

# models
from darts.models import NaiveSeasonal, StatsForecastAutoARIMA, ExponentialSmoothing, Prophet #local
from darts.models import LightGBMModel, RNNModel, NBEATSModel, TFTModel #global

# likelihood
from darts.utils.likelihood_models import GaussianLikelihood

# evaluation
from darts.metrics import mape

# settings
import warnings
warnings.filterwarnings("ignore")
import logging
logging.disable(logging.CRITICAL)

现在让我们加载数据集,这个数据集涉及需求预测,并来自Kaggle。如果你同意 Kaggle 上的比赛条款,那么你也可以下载这个数据集。

dataset = pd.read_csv('store_item_demand.csv')

# set the column type for column with date
dataset['date'] = pd.to_datetime(dataset['date'], format='%Y-%m-%d')
# sort values and reset index
dataset.sort_values(by=['date', 'store', 'item'], inplace=True)
dataset.reset_index(drop=True, inplace=True)
# creation of an auxiliary table with hierarchy and aggregated sales totals
hierarchy_df = dataset.groupby(['store', 'item'])[['sales']].sum()
hierarchy_df = hierarchy_df.reset_index(drop=False).sort_values(by=['sales'],
  ascending=False).reset_index(drop=True)

共有 10 家商店,每家商店有 50 种商品,总计 500 个时间序列。

让我们来看一下在所有商店商品组合中销售最少、中等和最多的 10 种商品。要找出时间序列中的关系,通常只需观察它,因为这已经能告诉我们很多信息,比如趋势或季节性,但往往还有更多。

fig, ax = plt.subplots(figsize=(30, 10))

for single_ts in list(np.arange(0, 10)) + list(np.arange(245, 255)) + list(np.arange(490, 500)):
    single_ts_df = pd.merge(dataset, hierarchy_df.loc[[single_ts], ['store', 'item']], how='inner', on=['store', 'item'])
    ax.plot(single_ts_df['date'], single_ts_df['sales'], color='black', alpha=0.25)
ax.set_xlim([dataset['date'].min(), dataset['date'].max()])
ax.set_ylim([0, dataset['sales'].max()])
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

这些时间序列之间有很多相似性。我们将稍后检查是否存在季节性(周几或一年中的周/月)或趋势,但你可以直观地想象,从所有时间序列中学习的模型将比仅从其历史中学习的模型更好地预测单个时间序列。

我为EDA(探索性数据分析)创建了数据集的副本。然后,我使用MinMAXScaler对时间序列进行缩放,使所有时间序列可以互相比较。最后,我创建箱型图,以检查是否存在趋势和季节性。

# make copy of df
dataset_scaled_EDA = dataset.copy()

# min max value calculation
dataset_scaled_EDA['min_sales'] = dataset_scaled_EDA.groupby(['store', 'item'])['sales'].transform(lambda x: x.min())
dataset_scaled_EDA['max_sales'] = dataset_scaled_EDA.groupby(['store', 'item'])['sales'].transform(lambda x: x.max())
# scale
dataset_scaled_EDA['sales_scaled'] = (dataset_scaled_EDA['sales'] - dataset_scaled_EDA['min_sales'])/(dataset_scaled_EDA['max_sales'] - dataset_scaled_EDA['min_sales'])
# add info about year, week of year and day of week
dataset_scaled_EDA['year'] = dataset_scaled_EDA['date'].dt.year
dataset_scaled_EDA['month'] = dataset_scaled_EDA['date'].dt.month
dataset_scaled_EDA['day_of_week'] = [d.strftime('%A') for d in dataset_scaled_EDA['date']]
dataset_scaled_EDA['day_of_week'] = pd.Categorical(dataset_scaled_EDA['day_of_week'], 
  categories=['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday', 'Sunday'], 
  ordered=True)

# visualize
fig, ax = plt.subplots(1, 3, figsize=(30, 10))
sns.boxplot(x='year', y='sales_scaled', data=dataset_scaled_EDA, ax=ax[0]).set(
    xlabel='Year', 
    ylabel='Scaled Sales'
)
ax[0].set_title('Box plot for years (trend)')
sns.boxplot(x='month', y='sales_scaled', data=dataset_scaled_EDA, ax=ax[1]).set(
    xlabel='Month', 
    ylabel='Scaled Sales'
)
ax[1].set_title('Box plot for months (seasonality)')
sns.boxplot(x='day_of_week', y='sales_scaled', data=dataset_scaled_EDA, ax=ax[2]).set(
    xlabel='Day of week', 
    ylabel='Scaled Sales'
)
ax[2].set_title('Box plot for day of week (seasonality)')
ax[2].set_xticklabels(ax[2].get_xticklabels(), rotation=30)
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

是的,有趋势和两种季节性。如果你认为这些时间序列很容易预测——你是对的。本文的目的是向你展示预测多个时间序列的最流行的方法。数据并不总是像股票市场指数那样简单,但那是另一个话题。

我认为这种探索足以理解模型应该学习哪些关系。

在这里,我们有多个时间序列(多个商店中的多个商品销售)。对于每一个时间序列,我们从 Pandas DataFrame 中创建一个TimeSeries对象。这种类型是 Darts 库中的模型所需的。然后将这些时间序列保存到列表中。

dataset_ts = dataset.copy()
dataset_ts = TimeSeries.from_group_dataframe(df=dataset_ts, 
                                             group_cols=['store', 'item'],
                                             time_col='date', 
                                             value_cols='sales')
dataset_ts

处理多个时间序列可能很有帮助,但通常也会带来问题。当我们只有一个时间序列时,我们有很多时间来处理它。查看它,验证趋势和季节性,并处理异常。我们可以优化我们的预测。对于多个序列,这种方法变得不可行。我们希望方法尽可能自动化,但这样可能会错过细节,比如异常,或者我们不应该以相同的方式处理每个时间序列。可能还有更多典型的问题:缺失数据、数据漂移和稀有事件(黑天鹅事件)。更多的序列可能对我们有帮助,因为我们的模型可以使用更多的数据,因此会有更多的代表性观察数据来识别特定模式。以我的工作为例——为了预测产品 X 在下周因促销造成的需求,我们的模型还可以利用其他促销的历史效果。

所以问题来了——如何预测多个时间序列?

你可能经常会问自己我是否有一个 多变量 多元 时间序列? 这些问题是合理的,但答案并不总是明确的。当你的时间序列来自单一过程相互关联相关,并且相互作用时,答案将是我的时间序列是 多元

当你在多个商店预测产品 X 的销售时,你有一个多时间序列,但当你有一个额外的产品 Y 时,对于单个商店来说,你有一个多变量时间序列,因为一个产品在商店中的销售可能影响另一个产品的销售,或者这只是一个假设。

如何评估预测结果?

在我们深入讨论不同的方法和模型之前,让我们先讨论如何衡量预测的质量。这是一个回归问题,因此我们仍然会将预测结果与真实值进行比较,这一点毫无疑问。你可以使用回归问题中熟悉的指标,如RMSE(均方根误差), MSE(均方误差), MAE(平均绝对误差),或者更典型的时间序列指标,如MAPE(平均绝对百分比误差), MARRE(平均绝对范围相对误差), 或MASE(平均绝对缩放误差)。进一步的讨论将使用MAPE作为评估指标。本文不是关于需求预测的,但由于数据和我的经验,有很多相关参考。因此,你可以考虑为该问题选择什么指标。始终选择能够反映业务目标的指标。在这篇文章中,Nicolas Vandeput 描述了需求预测中使用的 KPI 指标。

我们可以将这种方法扩展到多个序列,然后一次性计算所有序列的指标,或者分别对每个序列计算指标后再进行汇总。所以让我们继续探讨如何评估单个序列,然后再将其扩展到多个序列。

是的,这是一个回归问题,你可能会想为什么我解释这一点。在时间序列中,时间扮演了关键角色。数据是相对于时间排序的,观测值是相互关联的。因此,分割训练/测试集并使用交叉验证时无法使用随机化,因为这样会导致数据泄漏。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

预测模型的评估,图像由作者提供

首先,将数据分为训练集和测试集,很简单,对吧?模型在训练集上进行拟合,在测试集上进行测试。我们可以按比例划分,例如,测试集包括最后 20% 的数据,或者指定测试集的起始日期——可能由于业务考虑,测试集为过去一年数据会更重要。

在我们的例子中,测试集将是最后一年。

first_test_date = pd.Timestamp('2017-01-01')
train_dataset_ts, test_dataset_ts = [], []

for single_ts in tqdm(dataset_ts):
    # split into train and test tests
    single_train_ts, single_test_ts = single_ts.split_before(first_test_date)
    train_dataset_ts.append(single_train_ts)
    test_dataset_ts.append(single_test_ts)

训练集分离验证子集,用于选择超参数。对于在 Darts 中实现的模型,我们可以使用gridsearch方法,但对于基于神经网络的模型,推荐使用Optuna。gridsearch 方法实际上具备了我们所需的一切,无论选择哪个模型都能正常工作。

重要的是:

  • parameters ➜ 一个包含待检查超参数的字典

  • series ➜ TimeSeries 对象或包含 TimeSeries 对象的列表(如果算法是全局的)

  • start ➜ Pandas Timestamp 定义第一次预测发生的时间,或者 float 作为第一次预测之前观察值的比例。

  • forecast_horizon ➜ 模型预测的时间范围数量(int

  • stride ➜ 下一次预测之间的偏移量。为了使一切符合数据科学的艺术,并最能反映算法的实际操作,stride 应该等于 1。但请记住,在每一步之后你的算法都会重新训练,而且你的超参数网格可能有很多组合,这通常需要很长时间。因此,出于常识考虑,stride 可以大于 1,尤其是因为此目的在于选择最佳超参数,结果可能(但不一定)在 stride 等于 1 时与 stride 等于 5 时相同。

  • metric ➜ 函数,用于比较真实值和预测值。然后根据这个指标选择最佳超参数。

model.gridsearch(
  parameters={}, 
  series= ,
  start= , 
  forecast_horizon= ,
  stride= ,
  metric= 
)

那么如何在测试集上评估模型呢?我们有两种方法。

第一种方法是我们在训练集上拟合模型,并做一个覆盖测试集的预测。在以下示例中,我们将使用这个选项——计算速度更快,因为我们为每个时间序列生成一个预测。

第二种方法是我们测试不同的预测时间范围。我们拟合模型,做一个预测,重新训练,再做一个预测,以此类推,但我们必须在训练集结束之前开始训练算法,这样我们才能在相同的数据上比较时间范围——我们在比较相同的苹果。如果你从整个训练集开始训练,那么较远的时间范围有较少的观察值。在这种方法中,我们应该重新训练(特别是在局部模型中),因此这种方法中我们肯定会有更多的计算。

这也很简单,因为 darts 中的模型具有根据上述可视化返回历史预测的能力。你已经学习了一些在 gridsearch 方法中使用的变量,这些变量在这里也适用。然而,这里会有新的变量:

  • retrain ➜ 如果等于 True,则在每一步之后模型会重新训练,这最能反映现实情况。

  • overlap_end ➜ 如果等于 True,则预测可能超出测试集中的日期。如果我们对多个时间范围进行预测,而较远的时间范围超出测试集,而较近的时间范围不超出,这样的设置很有用。

  • last_points_only ➜ 如果等于 False,则返回所有时间范围的预测。

backtests_results = model.historical_forecasts(
  series= , 
  start= ,
  forecast_horizon= ,
  retrain=True,
  overlap_end=True,
  last_points_only=False,
)

然而,我们希望从感兴趣的时间范围中获取预测,并且在变量 backtests_results 中,有来自不同时间点的预测。要从特定的时间范围中获取预测,你可以使用我的函数 take_backtest_horizon

backtests_5W ➜ 有针对测试集每个点的预测,这些预测是在 5 周前做出的。

def take_backtest_horizon(backtests, horizon, last_horizon, observations):
    backtests_horizon_dates = [i.time_index[-1+horizon] for i in backtests[last_horizon-horizon:][:observations]]
    backtests_horizon_values = [i.values()[-1+horizon] for i in backtests[last_horizon-horizon:][:observations]]

    backtests_horizon = TimeSeries.from_times_and_values(times=pd.DatetimeIndex(backtests_horizon_dates), values=np.array(backtests_horizon_values))
    return backtests_horizon

backtests_5W = take_backtest_horizon(
  backtests=backtests_results, 
  horizon=5, # for which horizon we want a forecast, the number must be less than forecast_horizon
  last_horizon=forecast_horizon, # forecast_horizon from historical_forecasts method
  observations=len(test_series))

使用的变量

实际上,时间序列本身并不能完全自我解释。时间序列通常依赖于其他变量。这里我们没有这些变量,但了解如何在项目中遇到这些变量时进行区分是很有帮助的。如果我们没有通知模型即将到来的促销/降价,那么它将无法预测销售额的增加。如果你想预测股票公司的价格变化,添加来自技术分析或基本面分析的变量可能会有所帮助。这些变量依赖于价格,即你只能知道这些指标的过去值,毕竟,我们无法知道未来的公司财务报告或价格——我们希望进行预测 😃

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

时间序列算法使用的变量,图源作者

我们可以有也是时间序列的变量,即在不同时间点有不同的值,同时我们也可以有静态协变量(随时间保持不变),通常是分类变量。在我们的例子中,这将是商店 ID 和商品 ID。它们对全局模型非常重要。因为在每 100 个时间序列中,可能会有不同的关系,借助这些变量,你的模型可以区分不同的时间序列。

至于时间序列变量,我们可以区分协变量,其中有些已知未来(例如,我们可以知道未来的促销机制,也知道过去的促销机制)和仅已知过去的协变量(我们可以知道竞争产品的价格,但不知道它们未来的价格)。

使用的模型

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

局部算法与全局算法对比,图源作者

我们可以将机器学习分为监督学习、无监督学习和强化学习。进入详细内容后,我们可以将监督学习分为回归和分类。我们可以对时间序列预测进行类似的划分,即使用局部全局算法来预测这些时间序列。

局部算法是针对单一时间序列进行拟合的,模型仅能预测这一时间序列。更多的时间序列意味着更多的模型。在这里我们看到优缺点,模型简单,但对于许多时间序列来说,这种方法变得难以维护。

全局算法则是一个模型可以拟合多个时间序列。所以如果我们有多个时间序列,我们可以有一个模型来预测所有这些序列。这种方法显然更灵活,例如,你可以使用迁移学习。对于时间序列,这意味着你在一个不同的时间序列上拟合模型,而不是进行预测。这里有一个使用示例。另一个与全局模型相关的重要点,因为我可能会忘记,并且这非常重要——那就是时间序列缩放。最常见的方法是MinMaxScaler,但你可以使用更适合你数据的东西。然而,我不会在这里详细说明如何缩放时间序列,这肯定是另一个文章的话题。我们来考虑一下为什么我们应该缩放时间序列。答案可能很简单,许多全局算法是神经网络,这就是我们缩放数据的原因,就像我们对卷积神经网络的像素进行处理一样。然而,这还不是全部。然而,我们可以使用像随机森林这样的模型(非参数模型),我们仍然应该缩放它们。但停下来,为什么?毕竟,对于这些类型的模型,你不需要缩放变量。我们应该缩放时间序列的原因是为了让模型学习关系,而不是尺度。例如,对于季节性关系,可能是在夏季月份,值比冬季月份平均高 150%。另一个例子是,3 次显著增加后跟随一次下降。如果我们不缩放时间序列,模型很难学习这些关系。这是一种与表格数据中变量缩放略有不同的方法,因为在这里我们单独缩放每个时间序列。如果我们使用前面提到的MinMaxScaler,那么对于训练集中的每个时间序列,最大值是 1。所以让我们缩放我们的数据,这将被全局模型使用。

scaler = Scaler() # MinMaxScaler
train_dataset_ts_prepared = scaler.fit_transform(train_dataset_ts)
test_dataset_ts_prepared = scaler.transform(test_dataset_ts)
dataset_ts_prepared = scaler.transform(dataset_ts)

一会儿你将阅读到最受欢迎的局部全局算法。虽然不可能描述所有可能的算法,但有一些算法常被专家使用,并且通常能满足预期。

没有免费的午餐定理

没有一个答案——这个模型是最好的,不要使用其他模型。然而,如果你正在创建一个 MVP——最好从简单的东西开始。

在以下示例中,我没有在验证集上选择最佳超参数,而是使用了默认模型。所以如果你告诉我模型可能有更好的结果——我已经同意你的观点

局部模型

在我们深入具体的局部模型之前,我为你准备了函数,以便使用多处理来加快所有核心的计算速度。

forecast_horizons = len(test_dataset_ts[0])

def _backtests_local_estimator(_estimator, _ts_set, _split_date, _horizons, _single_forecast):
    model = _estimator
    if _single_forecast:
        model.fit(_ts_set.split_before(_split_date)[0])
        backtests_single_ts = model.predict(_horizons)

    else:
        backtests_single_ts = model.historical_forecasts(series=_ts_set, 
                                                         start=_split_date - np.timedelta64(_horizons-1, 'D'), 
                                                         verbose=False, 
                                                         overlap_end=False,
                                                         last_points_only=True, 
                                                         forecast_horizon=_horizons,
                                                         retrain=True)

    return backtests_single_ts

def backtests_multiple_local_estimators(estimator, multiple_ts_sets=dataset_ts, split_date=first_test_date, horizons=forecast_horizons, single_forecast=True):
    backtests_multiple_ts = Parallel(n_jobs=-1,
                                     verbose=5, 
                                     backend = 'multiprocessing',
                                     pre_dispatch='1.5*n_jobs')(
            delayed(_backtests_local_estimator)(
                _estimator=estimator,
                _ts_set=single_ts_set,
                _split_date=split_date,
                _horizons=horizons,
                _single_forecast=single_forecast
            )
        for single_ts_set in multiple_ts_sets
    )

    return backtests_multiple_ts

我将使用这个函数为本地模型生成测试集的单个预测。此外,你可以用它生成多个历史预测(第二种方法,变量single_forecast应设置为False)。不过,我在这里不这样做,因为这会花费很多时间。

如果你使用 Cluster 和 Spark,那么你可以使用 Spark UDF 来显著加快计算速度。

我知道,你可能想跳到模型部分。最后但同样重要的是——一个评估我们预测结果的函数。我将使用MAPE作为评估指标,然而,如果你在做需求预测项目,WMAPEMAE无疑更接近业务预期。

def get_overall_MAPE(prediction_series, test_series=test_dataset_ts):
    return np.round(np.mean(mape(actual_series=test_series, 
                                 pred_series=prediction_series, n_jobs=-1)),
                    2)

基准线

好吧,但为什么要做神经网络,如果一个更好的主意是从一年前预测值呢?这正是我们首先创建这样一个模型的原因。当你处理实际数据时,你最好也从这样的方式开始(它也可以是训练集中的最后一个值,NaiveDrift 如果有趋势,或几种简单方法的组合)。然后如果你转向更高级的方法,你可以评估它比简单方法好多少,因为例如,如果你从神经网络开始,其 MAPE 为 10%,那么我的问题(可能也是利益相关者的问题)是——这好吗?

我们的模型(一个时间序列=一个模型)将重复一年前的值。

backtests_baseline_model = backtests_multiple_local_estimators(estimator=NaiveSeasonal(K=365))
print(f'overall MAPE: {get_overall_MAPE(backtests_baseline_model)}%')

整体 MAPE: 22.42%

现在让我们在测试集中可视化时间序列的预测值和实际值,这些时间序列具有最高的总销售额。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

这看起来还不错。

ARIMA

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

视觉化 ARIMA 算法的工作原理,图像来源于作者

ARIMA 是一个统计模型,以其简单性而广受欢迎且强大。当你听到ARIMA时,它可能指的是这个模型,但也可能是ARIMA的扩展集合。这个集合包括ARIMAX(考虑额外变量)、SARIMA(考虑季节性)或VARIMA(用于多变量时间序列)。但让我们回到ARIMAAutoRegressive Integrated Moving Average),这就是一切的起点。如果你理解得很好,那么使用之前提到的模型应该不会有问题。

许多文章已经对此算法进行了阐述。我想给你这个模型背后的直观理解。我希望最终你能轻松地在代码中实现它,并理解它是如何工作的。我打算从最后开始。我自己记得在招聘过程中有过一个关于这个问题的疑问,当时我还没有理解它。ARMA 模型只能与平稳时间序列一起使用,因此我们有一个组件——积分I),它通常(但不总是)将非平稳时间序列转换为平稳时间序列ARMA是模型,而I部分负责为建模准备数据(当然,如果需要的话)。你应该知道几个问题的答案,那就是时间序列平稳或非平稳的含义,以及积分I)组件进行的转换类型。所以让我们从平稳性开始。

如果值的分布(均值和方差)在时间上是不变的,那么时间序列就是平稳的。

因此,如果存在趋势和/或季节性,那么时间序列就是非平稳的。要检查时间序列是否平稳,最简单的方法是将其可视化,并在图上添加移动平均和移动标准差。如果它们随时间保持不变(或接近不变),那么你可以得出你的时间序列是平稳的结论。这种方法可能显得天真,并且并不总是有效,因为如果对滚动统计数据使用过大的窗口,你可能会认为时间序列是平稳的,而实际上并非如此。另一种方法是将时间序列拆分成随机分区,对每个分区计算上述统计数据。最后一种方法是计算扩展的迪基-福勒ADF)检验。如果我们的时间序列仍然不是平稳的,我们需要使用积分I)组件怎么办?它通过差分来使时间序列平稳,即计算观察值之间的差异。如果我们的时间序列仍然不是平稳的呢?我们可以选择d的阶数,这表示我们对时间序列进行差分的次数。

这个长段落是关于积分I),它准备数据以供自回归AR)和移动平均MA)组件使用。AR是对最后p个值的线性回归,这些值被称为滞后。当前值与最后的值相关联,并且依赖于这些值。MA是补充性的,考虑了预测中q个最后的误差(假设为白噪声),以更好地预测当前的时间点。

选择AR的阶数p时,我们使用PACFPartial AutoCorrelation Function),而选择MA的阶数q时,我们使用ACFAutoCorrelation Function)。在大学课程之外,我们在实际操作中不太可能这样做,因为我们有AutoARIMA可以为我们选择pdq

让我们回到实践中,按照与基线模型相似的方式来实现。正如你可能已经了解到的,由于 Darts 库,这个过程非常简单。

backtests_arima = backtests_multiple_local_estimators(estimator=StatsForecastAutoARIMA())
print(f'overall MAPE: {get_overall_MAPE(backtests_arima)}%')

总体 MAPE: 28.18%

fig, ax = plt.subplots(figsize=(30, 10))
test_dataset_ts[0].plot(label='True value', color='black')
backtests_arima[0].plot(label='Forecast', color='green')
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

如果我选择参数m(季节性差分的周期),结果可能会更好。你可以尝试一下,并告诉我结果的变化情况。

指数平滑

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

可视化指数平滑(ETS)算法的工作原理,图片由作者提供

指数平滑是另一种用于单变量时间序列的模型家族。你可以在术语ETSE-误差,T-趋势,S-季节性)下找到“这个家族”。在这种方法中,观察值被赋予权重,较旧的观察值权重较低,因为它们按指数衰减。我们可以区分三种类型:一种简单的类型假设未来将类似于近期值,一种扩展类型处理趋势,最后一种还处理季节性。我将会在稍后描述这三种类型,不过现在先插一个小插曲。在M-4 比赛Makridakis 比赛,这是最著名的时间序列预测比赛)中,Slawek Smyl 获胜,他提出了ES-RNN,这是指数平滑递归神经网络的混合体。

现在我们回到话题的第一个类型,即简单指数平滑。作为基线模型,我们可以选择一个总是预测训练集中的最后一个值的模型,这是一种稍显天真的方法,但可以给我们带来不错的结果。另一种方法是计算整个训练集的平均值,但这样的话,最近的观察值和最古老的观察值会被赋予同等的重要性。指数平滑结合了这两种方法,赋予最近的观察值更大的权重,权重会随着观察值的古老程度呈指数下降,这意味着最古老的观察值将拥有最小的权重。它使用α参数,其范围在 0 到 1 之间。值越高,最新值对预测的影响越大。请查看上面的图形,那里也有公式。这些公式非常容易理解,通常这些模型能给出不错的结果。在我们进入更高级的模型之前,我们应该在这里稍作停留,因为如果一个简单的模型给出的结果与非常先进的模型(例如深度神经网络)相同,那么我们应该保持使用更简单的模型,因为它们的操作对我们来说更可预测,并且更多的人能够理解它们的运作。

SES 本质上无法处理数据中的趋势。如果存在上升趋势,则预测会低估,因为它没有包含这种增加。因此,我们有另一种模型,即双重指数平滑。它有一个额外的因素,用于考虑趋势的影响。我们使用β参数,它控制趋势变化的影响。因此我们有两个公式,一个用于水平(水平方程),另一个用于趋势(趋势方程)。

三重指数平滑也考虑了季节性。你可以将它称为 Holt-Winters 的季节性方法。这里另一个参数γ进入了公式。这种方法允许水平趋势季节性的模式随时间变化。像趋势一样,季节性可以是加法的或乘法的,但在这里我不会描述详细信息,假设你知道区别或可以轻松找到它们。我只是不想把这篇文章写成一本书,现在我希望你能顺利阅读整个文章。😉

backtests_exponential_smoothing  = backtests_multiple_local_estimators(estimator=ExponentialSmoothing())
print(f'overall MAPE: {get_overall_MAPE(backtests_exponential_smoothing)}%')

总体 MAPE: 31.88%

fig, ax = plt.subplots(figsize=(30, 10))
test_dataset_ts[0].plot(label='True value', color='black')
backtests_exponential_smoothing[0].plot(label='Forecast', color='green')
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

至于 ARIMA,我没有添加季节性信息。现在轮到你了,使用seasonal_periods参数为指数平滑添加这些信息。

Prophet

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

可视化 Prophet 算法的工作原理,作者提供的图片

Prophet 是 Facebook 在 2017 年的 大规模预测 论文中提出的。它既是一个模型,又是一个同名库。与之前的模型一样,你可以在 Darts 中找到它。该算法是Generalized Additive Model,因此预测是各个组件的总和。这些组件是g(t) — 趋势,s(t) — 季节性(每年、每周和每日),以及h(t) — 假期效应。

y(t) = g(t) + s(t) + h(t) + error(t)

第一个是趋势,它可以随着时间的变化而变化,并且不必始终保持不变。当学生开始学习时间序列分析课程时,他们通常处理的是简单的时间序列。在时间序列中,他们可以看到一个连续增长的趋势。然而,在实际数据中,趋势可能会发生多次变化。Prophet 实现了拐点(可以将其视为超参数,例如它们的数量、范围和先验尺度)。这些点是趋势的变化,例如,增加趋势 -> 拐点 -> 减少趋势 -> 拐点 -> 更强的减少趋势,等等。这种方法更接近我们通常在数据中看到的情况。这些拐点的位置是 Prophet 在你之前设置的。拐点之间的趋势函数可以是简单的回归。

接下来,我们有季节性函数,它是傅里叶级数

另一个功能是假期效应,它会对我们的预测值进行加减调整。你可以使用 Prophet 库提供的日期,或者定义你自己的事件。你可以想象,黑色星期五效应会显著影响销售额。此外,你可以考虑假期影响预测的日期范围,例如,圣诞节不会影响假期当天的销售,但会影响之前的几天(很多天)。

backtests_prophet = backtests_multiple_local_estimators(estimator=Prophet())
print(f'overall MAPE: {get_overall_MAPE(backtests_prophet)}%')

总体 MAPE:14.38%

fig, ax = plt.subplots(figsize=(30, 10))
test_dataset_ts[0].plot(label='True value', color='black')
backtests_prophet[0].plot(label='Forecast', color='green')
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

全球模型

现在我们转向一个模型用于所有时间序列的方法。这也被称为跨学习,因为该模型为了对时间序列 A 做出良好的预测,从时间序列 A 以及 B、C、D 等中学习了关系。

监督模型 ~ 时间序列作为回归问题

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

可视化如何为监督学习模型创建特征,图像来源于作者

现在让我们尝试将监督学习模型应用于时间序列预测。这并不新鲜,但它们往往能给出很好的结果,且比神经网络效果更佳(参见 M-5 竞赛中的最佳解决方案)。然而,回归模型并非专门针对时间序列,因此如果我们想使用它们,就需要将时间序列问题转换为机器学习问题。我在之前的文章中已经详细介绍了这一点,卖出与买入预测中有更多内容,但我也将描述如何使用众所周知的回归算法来解决这个问题。

在我开始之前提到的转换之前,我们首先需要缩放数据。之前我描述了这种转换对于全局模型的必要性。在这种情况下,我使用了MinMaxScaler

下一步是特征工程,这是一个重复了几次的转换过程。基于时间序列的历史,我们创建有助于模型更好地预测未来的特征。这些变量可以指代所选时间序列的近期历史,例如滞后值(对于每周数据,t-1W,t-2W,t-3W,等等)。另一个例子是滚动统计的计算,如中位数(可以是最近 4 周的中位数)、均值、最小值、最大值、标准差以及你将在窗口上计算的其他内容。如果数据中存在季节性,那么最好给模型一个关于时间点 t 的提示。我经常使用周期变量的正弦和余弦(在上面的可视化中是年份中的一天和一周中的一天)。

最后一步是选择一个模型,你有很多选择,包括线性回归、线性混合效应模型、随机森林、LightGBM 等。模型的选择取决于时间序列的性质和问题的复杂性。另一个问题可能是你想要一个模型还是和预测期数一样多的模型。当你选择一个模型时,需要考虑其弱点。例如,当你选择随机森林时,记住叶子节点计算的是均值,因此它不能超出训练范围。LightGBM 没有这个问题,因为它不计算朴素均值,而是在后台进行回归计算。

现在是时候回到实际操作部分并在代码中实现模型了。我选择了LightGBM作为模型。使用它在 Darts 中比我希望在没有这个库的情况下使用要简单得多。正如你在代码中看到的,我们使用了过去 14 天的滞后数据。我还添加了编码器,添加了模型使用的协变量,这些变量会自动添加,并且所有计算都在内部完成。

  • 周期性 — 添加 2 列,基于周期变量(如月份)的正弦余弦编码

  • 日期时间属性 — 基于日期时间变量添加标量

  • position — 基于时间序列索引,将相对索引位置添加为整数值,其中 0 设置在预测点。

model_lightGBM = LightGBMModel(lags=14, 
                               output_chunk_length=365, 
                               random_state=0,
                               multi_models=False, 
                               add_encoders={"cyclic": {"future": ["month"]}, 
                                             'datetime_attribute': {'future': ['dayofweek']},
                                             'position': {'past': ['relative'], 'future': ['relative']}, 
                                             'transformer': Scaler()
                 })

model_lightGBM.fit(series=train_dataset_ts_prepared)
backtests_lightGBM = model_lightGBM.predict(n=forecast_horizons, series=train_dataset_ts_prepared)
# The model predicts scaled values, but then we have to reverse the transformation
backtests_lightGBM = scaler.inverse_transform(backtests_lightGBM) # 
print(f'overall MAPE: {get_overall_MAPE(backtests_lightGBM)}%')

总体 MAPE: 15.01%

fig, ax = plt.subplots(figsize=(30, 10))
test_dataset_ts[0].plot(label='True value', color='black')
backtests_lightGBM[0].plot(label='Forecast', color='green')
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

结果非常有前景,尤其是因为这是一个适用于所有时间序列的模型。然而,基于我的经验,我想提醒你。这些类型的模型在特征工程上表现良好,这既是优势也是大缺陷。假设你使用了滞后和移动平均。现在你将预测一个时间序列的值,但它中有异常值——在预测点之前有几个大值。你的模型肯定会高估。当你创建变量时,尽量想象它们对模型的影响。

DeepAR

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

基于论文的 arxiv 和模型架构的截图

DeepAR 是由 Amazon 团队开发的一种深度学习算法。它旨在使用递归神经网络(RNNs)对时间序列数据中的复杂依赖关系进行建模。

正如我们在摘要中阅读到的(这与我非常接近):

例如,在零售业务中,预测需求对于在正确的时间和地点提供合适的库存至关重要。

该模型是自回归的,并使用蒙特卡洛样本生成概率预测。神经网络架构基于 LSTM 层。通过概率方法,我们不关心单一的良好预测,而是整个预测分布,来确定真实值可能出现的位置。DeepAR 不是直接使用 LSTMs 进行预测,而是利用 LSTMs 来参数化高斯似然函数。即,估计高斯函数的均值和标准差(θ = (μ, σ) 参数)。

DeepAR 支持未来已知协变量,我们没有这样的数据,但可以创建它们。作为这些特征,我创建了带有星期几和月份的独热编码(OHE)。可能更好的方法是使用正弦和余弦函数,我鼓励你进行实验,并将反馈意见告诉我。

day_series = datetime_attribute_timeseries(
    dataset_ts[0], attribute="weekday", one_hot=True, dtype=np.float32
)
month_series = datetime_attribute_timeseries(
    dataset_ts[0], attribute="month", one_hot=True, dtype=np.float32
)

day_month_series = day_series.concatenate(month_series, axis=1, ignore_static_covariates=True)

model_deepar = RNNModel(
    model="LSTM",
    hidden_dim=10,
    n_rnn_layers=2,
    dropout=0.2,
    batch_size=32,
    n_epochs=5,
    optimizer_kwargs={"lr": 1e-3},
    random_state=0,
    training_length=21,
    input_chunk_length=14,
    likelihood=GaussianLikelihood(),
    pl_trainer_kwargs={
        "accelerator": "gpu",
        "devices": [0]
    }
)

model_deepar.fit(series=train_dataset_ts_prepared, future_covariates=[train_day_month]*len(train_dataset_ts_prepared), verbose=True)

backtests_deepar = model_deepar.predict(n=forecast_horizons, series=train_dataset_ts_prepared, 
                                        future_covariates=[day_month_series]*len(train_dataset_ts_prepared), 
                                        num_samples=1000, verbose=True)
backtests_deepar = scaler.inverse_transform(backtests_deepar)
print(f'overall MAPE: {get_overall_MAPE(backtests_deepar)}%')

总体 MAPE: 19.35%

fig, ax = plt.subplots(figsize=(30, 10))
test_dataset_ts[0].plot(label='True value', color='black')
backtests_deepar[0].plot(label='Forecast', color='green')
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

N-BEATS

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

基于论文的 arxiv 和模型架构的截图

N-BEATS(用于可解释时间序列预测的神经基础扩展分析)是一种深度学习算法,但它不包含递归层,如 LSTM 或 GRU。

架构可能看起来复杂,但一旦你深入了解细节,它实际上非常简单,是由块组合而成,所有层都是前馈的。

我们从最小的元素——块开始,每个块有一个输入并生成两个输出。输入是回顾期。输出是预测和回溯。我想你对预测的概念很容易理解。回溯是预测,但对于回顾期,它是拟合值,并展示了块在回顾期窗口上的关系好坏。

让我们转到堆栈,或多个块的组合。如你所读,每个块有两个输出和一个输入。接下来的块负责预测残差——类似于提升森林模型中发生的情况,如 AdaBoost。在每一步中,块生成的回溯从前一个块的输入中减去。最后,所有块的预测结果都会聚合。此外,它是一个可解释的模型,你可以分解并查看趋势和季节性的影响。

现在让我们转到组合堆栈。这部分增加了模型的深度,并提供了更多了解复杂性的机会。

model_nbeats = NBEATSModel(
    input_chunk_length=178,
    output_chunk_length=356,
    generic_architecture=False,
    num_blocks=3,
    num_layers=4,
    layer_widths=512,
    n_epochs=10,
    nr_epochs_val_period=1,
    batch_size=800,
    model_name="nbeats_interpretable_run",
    random_state=0,
    pl_trainer_kwargs={
      "accelerator": "gpu",
      "devices": [0]
    }
)

model_nbeats.fit(series=train_dataset_ts_prepared, verbose=True)

backtests_nbeats = model_nbeats.predict(n=forecast_horizons, series=train_dataset_ts_prepared, verbose=True)
backtests_nbeats = scaler.inverse_transform(backtests_nbeats)
print(f'overall MAPE: {get_overall_MAPE(backtests_nbeats)}%')

整体 MAPE:13.18%

fig, ax = plt.subplots(figsize=(30, 10))
test_dataset_ts[0].plot(label='True value', color='black')
backtests_nbeats[0].plot(label='Forecast', color='green')
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

TFT

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

arxiv 和模型架构的截图,基于论文

时间融合变换器(TFT)是谷歌开发的用于时间序列预测的深度学习算法。它旨在通过结合变换器网络和自回归建模来建模时间序列数据中的复杂依赖关系和关系。

TFT 是最复杂的架构,使用了各种底层技术。它像洋葱一样,由多层组成。此外,根据我的经验,它相比上述模型学习时间最长。TFT 使用多头注意力块来寻找长期模式,但使用 LSTM 序列到序列编码器/解码器来寻找这些较短的模式。

model_tft = TFTModel(
    input_chunk_length=28,
    output_chunk_length=356,
    hidden_size=16,
    lstm_layers=1,
    num_attention_heads=3,
    dropout=0.1,
    batch_size=32,
    n_epochs=5,
    add_encoders={"cyclic": {"future": ["month"]}, 
                  'datetime_attribute': {'future': ['dayofweek']},
                  'position': {'past': ['relative'], 'future': ['relative']}, 
                  'transformer': Scaler()
                 },
    add_relative_index=False,
    optimizer_kwargs={"lr": 1e-3},
    random_state=0,
    pl_trainer_kwargs={
      "accelerator": "gpu",
      "devices": [0]
    }
)

model_tft.fit(series=train_dataset_ts_prepared, verbose=True)

backtests_tft = model_tft.predict(n=forecast_horizons, series=train_dataset_ts_prepared, 
                                  num_samples=1000, verbose=True)
backtests_tft = scaler.inverse_transform(backtests_tft)
print(f'overall MAPE: {get_overall_MAPE(backtests_tft)}%')

整体 MAPE:13.37%

fig, ax = plt.subplots(figsize=(30, 10))
test_dataset_ts[0].plot(label='True value', color='black')
backtests_tft[0].plot(label='Forecast', color='green')
plt.show()

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

总结

在这篇文章中,我想向你展示你可以选择哪些方法来预测多个时间序列。我提供了完全实用的代码,随意使用它们,不要犹豫与我联系。

这只是对主题的介绍。我认为数据科学家在供应链公司工作的相关主题还包括以下内容:

  • 层级预测及随后合并来自不同层级的预测,即层级调整。我们可以在商店级别进行预测,也可以在国家级别进行预测,但当我们将商店级别的预测聚合时,作为总和,我们希望得到与国家级别预测显示的相同结果。这就是层级调整重要的原因。

  • 另一个话题是 库存优化,即我们应该在库存中拥有多少产品,以避免没有产品库存的情况,但另一方面,我们也不希望某一产品库存数月。

有问题吗?

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

问题,图片由作者提供

我意识到在这篇文章中我触及了许多主题。我想给你一个可以采取的方向的指示。也许其中一些应该在这里更详细地描述,而其他的则在新文章中详细描述。请不要犹豫,你可以在 Linkedin 找到我。在未来的文章中,我希望它们能够涵盖详细的主题,并展示如何使用 PyTorch 库从头开始实现时间序列模型。

感谢你的时间!

数据集来源:

杂项{仅限需求预测内核,

作者 = {inversion},

标题 = {商店商品需求预测挑战},

出版商 = {Kaggle},

年份 = {2018},

网址 = {https://kaggle.com/competitions/demand-forecasting-kernels-only},

许可证 = CC

}

预测 API:一个使用 Django 和 Google Trends 的示例

原文:towardsdatascience.com/forecasting-api-an-example-with-django-and-google-trends-9b55046bd578

构建一个 web 应用程序以预测 Google Trends 的发展趋势。

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传 Davide Burba

·发布于 Towards Data Science ·阅读时间 14 分钟·2023 年 8 月 1 日

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

“Google Trends”,由 Giulia Roggia。经许可使用。

  • 简介

  • Django 模型

  • 服务:数据源,预处理,机器学习,任务

  • 交互层:序列化器,视图,端点

  • 结论

简介

什么是 Django?

Django 是一个高级 Python 网络框架。它设计得快速、安全、可扩展,因此是开发预期会增长复杂性的强大 web 应用程序的热门选择。有关 Django 的介绍,请参阅 这个教程

在这个示例中,我们将使用 Django Rest Framework(DRF),这是 Django 的一个扩展,旨在简化 REST API 的开发。有关 DRF 的介绍,请参阅 这个教程

需求

我们将通过列出一些假设的需求来开始设计我们的应用:

  • 总体目标:实现一个系统来预测未来时间序列的值。

  • 数据Google Trends 的周频数据,包含特征和目标,未来可能会扩展。数据应根据需求下载。

  • 预处理:仅使用滞后值。

  • 机器学习模型:一个全球 LightGBM 模型(如果你想了解更多关于全球与本地模型的区别,可以查看 这篇文章)。

  • 推断:生成在线预测(与批量预测相对),但不需要提供输入特征。

本教程中使用的完整代码可在 这里获取

设置环境

让我们从列出所需的依赖项开始。

python = "³.8"
Django = "⁴.2.1"
lightgbm = "³.3.5"
pandas = "².0.1"
djangorestframework = "³.14.0"
pytrends = "⁴.9.2"
drf-extensions = "⁰.7.1"

我们将使用poetry来管理依赖项,使用Docker来容器化项目。你可以在这里查看本项目中使用的 poetry 和 docker 文件。

快速开始

如果你想跳过前面的内容,直接开始使用应用程序,你可以运行以下命令:

# Clone the project.
git clone git@github.com:davide-burba/code-collection.git

# Move to the right folder.
cd code-collection/examples/api-example-django

# Launch the app.
docker compose up -d

# Apply the migrations.
docker compose exec django ./manage.py migrate

# Interactively create a (super)user for your app.
docker compose exec django ./manage.py createsuperuser 

并连接到localhost:8000/gtrends

Django 模型

在本节中,我们列出了示例中使用的 Django 模型。

时间序列

对于一个预测系统,我们需要处理时间序列数据。通常我们只需要两个模型:一个用于识别每个时间序列,另一个用于存储其值。但由于Google Trends 的历史数据可能会因为归一化而每天发生变化,我们还需要对数据进行版本控制,这会导致额外的一个模型。

我们还创建了一个模型,用于列出不同的数据源(目前只有 Google Trends)。

class TimeSeries(models.Model):
    name = models.CharField(unique=True, max_length=64)
    source = models.CharField(max_length=32, choices=DataSource.choices)

class TSVersion(models.Model):
    timeseries = models.ForeignKey(TimeSeries, on_delete=CASCADE)
    created_at = models.DateTimeField(auto_now_add=True)
    expired = models.BooleanField(default=False)

class TSValue(models.Model):
    version = models.ForeignKey(TSVersion, on_delete=CASCADE)
    time = models.DateTimeField()
    value = models.FloatField()

class DataSource(models.TextChoices):
    GOOGLE_TRENDS = "GOOGLE_TRENDS"

配置

要训练一个监督模型,我们需要一组特征和目标。我们可以将这些信息存储在一个“数据配置”中。这个过程通过以下模型完成。

class DataConfig(models.Model):
    name = models.CharField(unique=True, max_length=64)

class DataFeatures(models.Model):
    config = models.ForeignKey(
        DataConfig, on_delete=CASCADE, related_name="features"
    )
    timeseries = models.ForeignKey(TimeSeries, on_delete=PROTECT)

class DataTargets(models.Model):
    config = models.ForeignKey(
        DataConfig, on_delete=CASCADE, related_name="targets"
    )
    timeseries = models.ForeignKey(TimeSeries, on_delete=PROTECT)

同样,我们需要存储预处理和机器学习(ML)模型的配置。为了简化,我们将这些存储在一个JSONField中:

class PreprocessingConfig(models.Model):
    name = models.CharField(unique=True, max_length=64)
    params = models.JSONField()

class MLConfig(models.Model):
    params = models.JSONField()

ML 模型

在这里,我们区分了 ML 配置和 ML 模型。ML 配置包含关于 LightGBM 参数的数据,而 ML 模型包含所有关于数据、预处理和 ML 配置的信息。

一个 ML 模型可以在不同的数据集上进行估计(由于新数据的到来),因此我们还需要对 ML 模型进行版本控制。这会导致以下两个模型:

class MLModel(models.Model):
    name = models.CharField(unique=True, max_length=64)
    ml_config = models.ForeignKey(MLConfig, on_delete=PROTECT)
    data_config = models.ForeignKey(DataConfig, on_delete=PROTECT)
    preprocess_config = models.ForeignKey(
        PreprocessingConfig, on_delete=PROTECT
    )

class MLModelVersion(models.Model):
    ml_model = models.ForeignKey(MLModel, on_delete=CASCADE)
    ml_file = models.FileField(upload_to="ml_models")
    created_at = models.DateTimeField(auto_now_add=True)
    metadata = models.JSONField()

MLModelVersion保存了指向 LightGBM 工件的链接在ml_file中,以及关于它所训练数据的信息在metadata中。工件存储在 django 设置模块中指定的位置:例如,可能是文件系统中的一个文件夹或云中的一个 S3 桶。

服务

在本节中,我们描述了包含应用逻辑的服务。根据Django-StyleGuide,将其与视图分开是最佳实践。

数据源

由于数据源预计将来会增长,我们将使用灵活的设计:

  • 一个抽象的DataSource类,定义了一个接口。

  • 一个GTrendSource类,继承自DataSource并实现下载 Google Trends 数据的细节。

  • 一个download_data工厂用来构建DataSource子类。

这会产生以下模块:

import datetime as dt
from abc import ABC, abstractmethod

import pandas as pd
from gtrends.models import TimeSeries
from pytrends.request import TrendReq

def download_data(timeseries: TimeSeries) -> pd.DataFrame:
    return DATASOURCE_MAPtimeseries.source.download()

class DataSource(ABC):
    def __init__(self, timeseries: TimeSeries):
        self.timeseries = timeseries

    @abstractmethod
    def download(self) -> pd.DataFrame:
        """Returns a dataframe with time as index and value as column."""

class GTrendSource(DataSource):
    START_DATE = "2022-01-01"

    def download(self) -> pd.DataFrame:
        # Download data.
        name = self.timeseries.name
        data = self.download_interest_over_time(name)
        # Format new data.
        data = (
            data[~data.isPartial]
            .reset_index()
            .rename(columns={name: "value", "date": "time"})
            .set_index("time")
        )[["value"]]
        data["value"] = data["value"].astype(float)
        return data

    @classmethod
    def download_interest_over_time(cls, search_term: str) -> pd.DataFrame:
        """Download Google Trends data."""
        pytrends = TrendReq()
        timeframe = (
            cls.START_DATE + " " + dt.datetime.now().strftime("%Y-%m-%d")
        )
        pytrends.build_payload([search_term], timeframe=timeframe)
        return pytrends.interest_over_time()

DATASOURCE_MAP = {
    "GOOGLE_TRENDS": GTrendSource,
}

预处理

一旦数据加载完成,我们希望预处理数据以便后续训练模型或推断未来值。让我们从定义一个接口类开始。

class BasePreprocessor(ABC):
    @abstractmethod
    def build_x_y(self, data: Dict) -> Tuple[pd.DataFrame, pd.Series]:
        """Return features and target ready for training."""

    @abstractmethod
    def build_x_latest(self, data: Dict) -> pd.DataFrame:
        """Return only latest values of features, useful for inference"""

目前我们只考虑用于特征工程的滞后值。让我们创建一个实现滞后逻辑的辅助函数。

def _build_lags(
    df: pd.DataFrame, column: str, lags: List[int], prefix: str
) -> pd.DataFrame:
    return pd.concat(
        [
            df[[column]]
            .shift(lag)
            .rename(columns={column: f"{prefix}_lag_{lag}"})
            for lag in lags
        ],
        axis=1,
    )

现在我们需要创建一个实现基类抽象方法的类。让我们首先定义它的属性。

@dataclass
class Preprocessor(BasePreprocessor):
    horizon: int
    target_lags: List[int]
    feature_lags: List[int]

请注意,目标本身可以作为特征使用。对于作为特征使用的目标,我们分配一个通用前缀“target”。当有多个目标时,这可能是有益的,因为它们将堆叠在相同的列中,从而减少特征的数量。

这里是构建滞后目标和特征的实现:

def _build_x_lags_targets(
        self, target_data: Dict
    ) -> Optional[pd.DataFrame]:
        if not self.target_lags:
            return None

        x = []
        for df in target_data.values():
            x.append(
                _build_lags(
                    df=df,
                    column="value",
                    lags=self.target_lags,
                    prefix="target",
                )
            )
        return pd.concat(x, axis=1)

def _build_x_lags_features(
        self, feature_data: Dict, target_data: Dict
    ) -> Optional[pd.DataFrame]:
        if not self.feature_lags:
            return None

        x = []
        for name, df in feature_data.items():
            x.append(
                _build_lags(
                    df=df,
                    column="value",
                    lags=self.feature_lags,
                    prefix=name,
                )
            )

        # Concat features on axis 1.
        x = pd.concat(
            [df.reset_index().drop(columns=["ts_name"]) for df in x], axis=1
        )

        # Use target to "reindex" on axis 0.
        for_reindex = pd.concat(target_data.values(), axis=1).reset_index()
        x = pd.merge(for_reindex, x, how="left", on="time")

        return x.drop(columns=["value"]).set_index(["time", "ts_name"])

现在我们可以将它们封装在一个build_x方法中:

def build_x(self, data: Dict) -> pd.DataFrame:
    target_data = data["targets"]
    feature_data = data["features"]

    # Build x_target and x_features.
    x_targ = self._build_x_lags_targets(target_data)
    x_feat = self._build_x_lags_features(feature_data, target_data)
    # Combine x_target and x_features.
    if x_feat is None and x_targ is None:
        raise ValueError("Cannot have no target lags and no feature lags.")
    elif x_feat is None:
        return x_targ
    elif x_targ is None:
        return x_feat
    return pd.merge(
        x_targ, x_feat, left_index=True, right_index=True, how="left"
    )

为了构建目标,我们只需将其移动horizon时间步:

def build_y(self, target_data: Dict) -> pd.DataFrame:
    y = {}
    for name, df in target_data.items():
        y[name] = (
            df["value"]
            .shift(-self.horizon)
            .rename(f"horizon_{self.horizon}")
        )
    return pd.concat(y.values())

现在我们可以为所需的抽象方法提供实现:

def build_x_y(self, data: Dict) -> Tuple[pd.DataFrame, pd.Series]:
    # Build x and y.
    x = self.build_x(data)
    y = self.build_y(data["targets"])
    # Align x indexes with y indexes.
    x = pd.merge(y, x, left_index=True, right_index=True, how="left")
    x = x.iloc[:, 1:]
    # Drop missing values generated by lags/horizon.
    idx = ~(x.isnull().any(axis=1) | y.isnull())
    x = x.loc[idx]
    y = y.loc[idx]
    return x, y

def build_x_latest(self, data: Dict) -> pd.DataFrame:
    x = self.build_x(data)
    return x[x.index == x.index.max()]

ML

让我们定义一个非常简单的 ML 模块,提供两个函数用于转储和加载 LightGBM 模型(如果你想知道我们为何选择使用 LightGBM,你可以查看这篇文章)。

import lightgbm

def save_engine(model, path):
    model.booster_.save_model(path)

def load_engine(path):
    return lightgbm.Booster(model_file=path)

任务

在定义 API 接口之前,让我们定义一些将在端点中使用的tasks

更新时间序列

我们首先需要存储一些时间序列数据。我们不希望每次下载一个新的数据点时都增加一个新版本到数据库中,因为这样会快速增长并产生大量重复数据。相反,我们希望将新下载的数据与最新版本进行比较,只有当数据历史不匹配时才创建新版本。这是在以下代码片段中完成的:

from typing import Tuple
import pandas as pd
from gtrends.services.data_sources import download_data
from gtrends.models import TimeSeries, TSValue, TSVersion

def update_timeseries(timeseries: TimeSeries) -> Tuple[bool, int]:
    """Update timeseries values.

    Either add the new values (if past values are the same) to the latest
    version, or create a new version.

    Args:
        timeseries: A timeseries object.

    Returns:
        A pair with:
            - bool: True if it created a new version, False otherwise
            - int: The number of new values added.
    """
    new_data = download_data(timeseries)

    # Assign version.
    new_version = True
    versions = timeseries.tsversion_set.order_by("created_at")
    if versions:
        version = versions.last()
        old_data = _build_old_data(version)

        if _is_old_data_in_new_data(old_data, new_data):
            # If old values match, just keep the new values.
            new_version = False
            new_data = new_data.loc[~new_data.index.isin(old_data.index)]
        else:
            # Else, set the old version to expired.
            version.expired = True
            version.save()

    if new_version:
        version = TSVersion(timeseries=timeseries)
        version.save()

    # Store new data.
    objs = [
        TSValue(version=version, time=d[0], value=d[1].value)
        for d in new_data.iterrows()
    ]
    TSValue.objects.bulk_create(objs)
    return new_version, len(objs)

def _build_old_data(version: TSVersion) -> pd.DataFrame:
    old_data = pd.DataFrame(
        version.tsvalue_set.values("time", "value")
    ).set_index("time")
    old_data.index = old_data.index.tz_localize(None)
    return old_data

def _is_old_data_in_new_data(
    old_data: pd.DataFrame, new_data: pd.DataFrame
) -> bool:
    if old_data.index.isin(new_data.index).all() and new_data.loc[
        old_data.index
    ].equals(old_data):
        return True
    return False

我们还可以添加一个包装器,以方便地更新数据库中的所有时间序列:

def update_all_timeseries():
    """Update all time-series values."""
    for ts in TimeSeries.object.all():
        update_timeseries(ts)

加载数据

现在我们已经存储了一些时间序列数据,我们希望能够加载它以便进行预处理。

from typing import Dict, List
import pandas as pd
from gtrends.models import TimeSeries

def load_data(
    target_ts: List[TimeSeries], feature_ts: List[TimeSeries]
) -> Dict[str, Dict[str, pd.DataFrame]]:
    targ_feat = {
        "targets": target_ts,
        "features": feature_ts,
    }
    data = {"targets": {}, "features": {}}
    metadata = {"targets": {}, "features": {}}
    for key, items in targ_feat.items():
        for item in items:
            ts = item.timeseries
            version = ts.tsversion_set.last()
            values = version.tsvalue_set.all().values("time", "value")

            df = pd.DataFrame(values)
            df["ts_name"] = ts.name
            df = df.set_index(["time", "ts_name"])

            data[key][ts.name] = df
            metadata[key][ts.name] = version.id

    return data, metadata

预处理

让我们将预处理封装成几个任务。为了简化,我们直接使用Preprocessor类;如果以后想添加另一个预处理类,我们可以使用与数据源相同的工厂模式。

def preprocess(data: Dict, prep_params: Dict) -> Tuple[pd.DataFrame, pd.Series]:
    return Preprocessor(**prep_params).build_x_y(data)

def build_x_latest(data: Dict, prep_params: Dict) -> pd.DataFrame:
    return Preprocessor(**prep_params).build_x_latest(data)

训练

现在我们已经有了预处理的数据,我们终于可以训练我们的 LightGBM 模型了。请注意,这是更简单的步骤!

from lightgbm import LGBMRegressor

def train(x, y, model_params):
    model = LGBMRegressor(**model_params)
    model.fit(x, y)
    return model

存储 ML 模型

现在我们可以创建一个任务,将 LightGBM 引擎保存以创建一个新的MLModelVersion。注意,我们首先将引擎转储到一个临时文件中。这只是一个技巧,以避免在任务中硬编码存储类型,而是通过 Django 设置动态处理。有关此主题的更多信息,请查看FileField 文档

import os
from typing import Dict
from uuid import uuid4

from django.core.files import File
from gtrends import models
from gtrends.ml import save_engine
from lightgbm import LGBMRegressor

def save_mlmodelversion(
    engine: LGBMRegressor, ml_model: models.MLModel, metadata: Dict
) -> models.MLModelVersion:
    filename = str(uuid4()).replace("-", "")[:10] + ".txt"
    tmp_path = "/tmp/" + filename
    save_engine(engine, tmp_path)

    with open(tmp_path, "r") as f:
        ml_model_file = File(f, filename)
        ml_model_version = models.MLModelVersion(
            ml_model=ml_model,
            ml_file=ml_model_file,
            metadata=metadata,
        )
        ml_model_version.save()

    os.remove(tmp_path)
    return ml_model_version

管道

我们现在可以完成训练和推断管道。由于这些管道速度较快,可以将它们简单地封装在一个函数中。然而,对于长时间运行的管道或高流量的情况,建议使用任务管理器,如 AirFlowCelery

这是训练管道:

from gtrends.models import MLModel, MLModelVersion
from gtrends.services.tasks import (
    load_data,
    preprocess,
    save_mlmodelversion,
    train,
)

def train_pipeline(ml_model: MLModel) -> MLModelVersion:
    target_ts = ml_model.data_config.targets.all()
    feature_ts = ml_model.data_config.features.all()
    prep_params = ml_model.preprocess_config.params
    model_params = ml_model.ml_config.params

    data, metadata = load_data(target_ts, feature_ts)
    x, y = preprocess(data, prep_params)
    engine = train(x, y, model_params)
    ml_model_version = save_mlmodelversion(engine, ml_model, metadata)

    return ml_model_version

这是推断管道:

from typing import Dict
from gtrends.models import MLModel
from gtrends.services.ml import load_engine
from gtrends.services.tasks import build_x_latest, load_data

def inference_pipeline(ml_model: MLModel) -> Dict:
    if ml_model.mlmodelversion_set.count() == 0:
        raise ValueError("No model has been trained yet!")

    target_ts = ml_model.data_config.targets.all()
    feature_ts = ml_model.data_config.features.all()
    prep_params = ml_model.preprocess_config.params

    data, _ = load_data(target_ts, feature_ts)
    x = build_x_latest(data, prep_params)
    engine = load_engine(ml_model.mlmodelversion_set.last().ml_file.path)
    y_pred = engine.predict(x)

    # Format predictions.
    horizon = ml_model.preprocess_config.params["horizon"]
    preds = {
        name: {"last_date": time, "prediction": pred, "horizon": horizon}
        for time, name, pred in zip(
            x.index.get_level_values("time"),
            x.index.get_level_values("ts_name"),
            y_pred,
        )
    }
    for k, v in preds.items():
        last_value = data["targets"][k].loc[v["last_date"]]["value"].item()
        preds[k]["last_value"] = last_value
        preds[k]["predicted_delta"] = v["prediction"] - last_value

    return preds

交互层

序列化器、视图和端点代表应用程序的“交互”层,这是最外部的层,并定义了用户如何与应用程序交互。

序列化器

序列化器负责序列化和反序列化数据。以下是一个简单的示例:

class TSVersionSerializer(serializers.ModelSerializer):
    class Meta:
        model = models.TSVersion
        fields = "__all__"

为了简洁起见,我们只讨论那些非平凡的视图。

时间序列

对于时间序列,我们重写默认的创建方法,以便在创建 TimeSeries 对象时下载第一个版本。

class TimeSeriesSerializer(serializers.ModelSerializer):
    class Meta:
        model = models.TimeSeries
        fields = "__all__"

    def create(self, validated_data):
        ts = super().create(validated_data)
        # Download data for the newly created timeseries.
        update_timeseries(ts)
        return ts

数据配置

对于数据配置,我们添加了特征和目标的链接,以便以后从单个端点创建数据配置。

class DataConfigSerializer(serializers.ModelSerializer):
    features = DataFeaturesSerializer(many=True)
    targets = DataTargetsSerializer(many=True, allow_null=False)

    class Meta:
        model = models.DataConfig
        fields = ("id", "name", "features", "targets")

    def validate_targets(self, value):
        if len(value) == 0:
            raise serializers.ValidationError("Must have at least one target.")
        return value

    def create(self, validated_data):
        config = models.DataConfig.objects.create(**validated_data)

        features_data = validated_data.pop("features")
        targets_data = validated_data.pop("targets")
        for feature_data in features_data:
            models.DataFeatures.objects.create(config=config, **feature_data)
        for target_data in targets_data:
            models.DataTargets.objects.create(config=config, **target_data)

        return config

预处理

让我们在预处理配置上添加一个验证步骤,以检查是否具有预期的参数:

class PreprocessingConfigSerializer(serializers.ModelSerializer):
    class Meta:
        model = models.PreprocessingConfig
        fields = "__all__"

    def validate_params(self, value):
        expected_keys = {"horizon", "target_lags", "feature_lags"}
        if set(value.keys()) != expected_keys:
            raise serializers.ValidationError(f"Expected keys: {expected_keys}")

        if not isinstance(value["horizon"], int):
            raise serializers.ValidationError("Horizon must be an int.")

        for key in ["target_lags", "feature_lags"]:
            if not isinstance(value[key], list) or not all(
                isinstance(lag, int) for lag in value[key]
            ):
                raise serializers.ValidationError(f"{key} must be a list[int].")

        return value

ML 配置

同样,让我们检查为 ML 配置提供的参数是否有效:

class MLConfigSerializer(serializers.ModelSerializer):
    class Meta:
        model = models.MLConfig
        fields = "__all__"

    def validate_params(self, value):
        # At the moment we only support LightGBM model.
        valid_params = LGBMRegressor().get_params()
        invalid_params = [k for k in value if k not in valid_params]
        if invalid_params:
            raise serializers.ValidationError(f"Invalid: {invalid_params}")
        return value

视图

现在让我们定义用户如何与 API 交互。我们将定义一些动作来触发数据检索或 ML 模型训练。同样,这种做法是可以接受的,因为这些动作执行速度很快,并且我们希望在项目开始时保持简单(如果你对这种哲学感兴趣,可以查看 YAGNI)。然而,随着项目的增长,可能需要将这些动作异步化并使用任务管理器。

请注意,为了保持版本控制,我们不允许使用 putpatch 方法。

时间序列

让我们为时间序列建立一个视图。我们希望能够:

  • 创建时间序列并列出它们

  • 获取最新版本的值

  • 更新值

class TimeSeriesViewSet(viewsets.ModelViewSet):
    queryset = models.TimeSeries.objects.prefetch_related("tsversion_set").all()
    serializer_class = serializers.TimeSeriesSerializer
    http_method_names = ["get", "post", "head", "delete"]

    @action(detail=True, url_path="latest-values")
    def latest_values(self, request, pk):
        """Get timeseries values for the last version."""
        versions = self.queryset.get(pk=pk).tsversion_set.order_by("created_at")
        if not versions:
            return Response([])
        return Response(versions.first().tsvalue_set.values("time", "value"))

    @action(detail=True, url_path="update-values")
    @transaction.atomic
    def update_values(self, request, pk):
        """Update values for one timeseries."""
        timeseries = self.queryset.get(pk=pk)
        try:
            new_version, how_many = update_timeseries(timeseries)
        except ResponseError:
            return Response(
                "Data download failed!", status=status.HTTP_400_BAD_REQUEST
            )
        msg = f"Data updated, {how_many} new values. "
        if new_version:
            msg += " New version created."
        else:
            msg += " Updated latest version."
        return Response(msg)

    @action(detail=False, url_path="update-all-values")
    @transaction.atomic
    def update_all_values(self, request):
        """Update values for all timeseries."""
        try:
            for timeseries in self.queryset:
                update_timeseries(timeseries)
        except ResponseError:
            return Response(
                "Data download failed!", status=status.HTTP_400_BAD_REQUEST
            )
        msg = f"Data updated, {len(self.queryset)} timeseries updated."
        return Response(msg)

我们还希望列出一个时间序列的不同版本,并检查其值。

class TSVersionViewSet(viewsets.ModelViewSet):
    queryset = models.TSVersion.objects.prefetch_related("tsvalue_set").all()
    serializer_class = serializers.TSVersionSerializer
    http_method_names = ["get", "head"]

    @action(detail=True)
    def values(self, request, pk, **kwargs):
        """Get timeseries values."""
        values = self.queryset.get(pk=pk).tsvalue_set.values("time", "value")
        return Response(values)

配置

我们希望列出、添加和删除配置。它们的视图是直接的。

class MLConfigViewSet(viewsets.ModelViewSet):
    queryset = models.MLConfig.objects.all()
    serializer_class = serializers.MLConfigSerializer
    http_method_names = ["get", "post", "head", "delete"]

class DataConfigViewSet(viewsets.ModelViewSet):
    queryset = models.DataConfig.objects.all()
    serializer_class = serializers.DataConfigSerializer
    http_method_names = ["get", "post", "head", "delete"]

class PreprocessingConfigViewSet(viewsets.ModelViewSet):
    queryset = models.PreprocessingConfig.objects.all()
    serializer_class = serializers.PreprocessingConfigSerializer
    http_method_names = ["get", "post", "head", "delete"]

ML 模型

对于 ML 模型,我们希望能够创建和列出它们,并执行训练和推断。

class MLModelViewSet(viewsets.ModelViewSet):
    queryset = (
        models.MLModel.objects.select_related("preprocess_config", "ml_config")
        .prefetch_related(
            "data_config__targets",
            "data_config__features",
            "mlmodelversion_set",
        )
        .all()
    )
    serializer_class = serializers.MLModelSerializer
    http_method_names = ["get", "post", "head", "delete"]

    @action(detail=True)
    def train(self, request, pk, **kwargs):
        ml_model = self.queryset.get(pk=pk)
        ml_model_version = train_pipeline(ml_model)
        return Response(
            serializers.MLModelVersionSerializer(ml_model_version).data
        )

    @action(detail=True)
    def predict(self, request, pk, **kwargs):
        ml_model = self.queryset.get(pk=pk)
        if ml_model.mlmodelversion_set.count() == 0:
            return Response(
                "No model has been trained yet!",
                status=status.HTTP_404_NOT_FOUND,
            )

        predictions = inference_pipeline(ml_model)
        return Response(predictions)

我们还希望列出特定 MLModel 的所有版本。

class MLModelVersionViewSet(viewsets.ModelViewSet):
    queryset = models.MLModelVersion.objects.all()
    serializer_class = serializers.MLModelVersionSerializer
    http_method_names = ["get", "head"]

端点

最后,我们可以将刚创建的视图链接到一组端点:

  • timeseries 用于时间序列

  • timeseries/<id>/versions 用于时间序列版本

  • model-configdata-configpreprocessing-config 用于配置对象

  • model 用于 ML 模型

  • model/<id>/versions 用于 ML 模型版本

这在以下模块中完成:

from django.urls import include, path
from gtrends import views
from rest_framework_extensions.routers import ExtendedDefaultRouter

router = ExtendedDefaultRouter()

router.register("model-config", views.MLConfigViewSet)
router.register("data-config", views.DataConfigViewSet)
router.register("preprocessing-config", views.PreprocessingConfigViewSet)

models = router.register("model", views.MLModelViewSet)
models.register(
    "versions",
    views.MLModelVersionViewSet,
    basename="version",
    parents_query_lookups="model_id",
)

timeseries = router.register("timeseries", views.TimeSeriesViewSet)
timeseries.register(
    "versions",
    views.TSVersionViewSet,
    basename="version",
    parents_query_lookups="timeseries_id",
)

urlpatterns = [
    path("", include(router.urls)),
]

请注意,我们不需要为数据检索或训练定义端点,因为这些会自动添加为“额外操作”。例如,训练模型的端点是 model/<id>/train

结论

我们在 Django 中看到了一个预测系统的示例实现。你现在可以开始玩玩这个系统了!要开始,请按照文章开头的说明进行操作,以便你可以连接到 localhost:8000/gtrends 并玩转这个应用。

完成设置并训练好模型后,你应该会看到类似于这样的结果(这里我们预测的是术语 “forecasting”):

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

当你预测未来值时,你应该会看到类似这样的结果:

外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传

就这些,你现在可以根据需要实现不同的变体!

这里是可能改进的非详尽列表:

  • 实现任务调度

  • 添加更多数据源

  • 添加一个 回测 端点

  • 添加更多预处理选项

  • 集成 MLFlow 用于跟踪和版本控制

  • 允许用一个模型训练多个预测视角

  • 设置适当的网络服务器(生产环境所需)

喜欢这篇文章吗? 查看我的其他文章 并关注我以获取更多内容! 点击这里 以无限阅读文章,并且对你没有额外费用的情况下支持我❤️

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值