回到未来：使用马尔可夫转移矩阵分析时间序列数据
概念概述及实际应用
·
关注 发表在 Towards Data Science ·11 分钟阅读·2023 年 2 月 1 日
–
图片由 Oto Godfrey 和 Justin Morton 提供，来自 Wikimedia Commons：根据 CC-BY-SA-4.0 许可自由使用
在本文中，我们将探讨如何通过使用马尔科夫转移矩阵重新构造时间序列数据，来获得有趣的描述性洞察以及优雅的预测、回溯和收敛分析方法。就像科幻经典《回到未来》中 Doc 改装的 DeLorean 时间机器一样，前进和后退时间。
注意： 以下章节中的所有方程式和图示均由本文作者创建。
基本构建块
设 E 定义组成时间序列数据的 k 个独特事件。例如，一个时间序列可能由以下三个基本且独特的事件组成，这些事件代表在离散时间步上绘制数据时观察到的路径轨迹类型：下跌、平稳 和 上涨。设 S 定义一个长度为 n（表示离散时间步）的序列，包含在 E 中定义的事件，代表部分或全部数据。例如，序列 [上涨、下跌、上涨、平稳、上涨] 代表五个时间步的数据。
现在可以定义一个维度为k²的马尔科夫转移矩阵 M，其中每个元素 M(i, j) 描述了从时间步 t 的事件 E(i) 转移到时间步 t+1 的事件 E(j) 的概率。换句话说，M(i, j) 表示在连续时间步之间转移的条件概率。从图论的角度看，事件 E(i) 和 E(j) 可以被视作通过有向边 E(i) → E(j) 连接的节点，如果 E(i) 在时间序列数据中被 E(j) 跟随；那么马尔科夫转移矩阵 M 本质上表示了图中节点所描绘事件的邻接矩阵（或共现矩阵）的归一化版本。
接下来，让我们看看这些基本构建块可以做些什么。
转移矩阵的应用：一个简单的例子
假设我们有以下涵盖 11 个连续时间步的原始时间序列数据：[1, 2, -2, -1, 0, 0, 2, 2, 1, 2, 3]。使用上述描述的路径轨迹简化视图，我们可以将数据转换为描述相邻时间步之间转移的以下 10 个事件序列：[上涨、下跌、上涨、上涨、平稳、上涨、平稳、下跌、上涨、上涨]。
我们现在可以构造以下邻接矩阵，以捕捉事件序列中共现的模式：
元素 A(i, j) 表示在我们的事件序列中，在某个时间步 t 上发生事件 i 后，接着发生事件 j 在时间步 t+1 的次数；i 和 j 分别是行和列索引。请注意，行表示事件的顺序为 up、flat 和 down，从上到下，而列表示从左到右的顺序相同。例如，A 的左上角元素表示在给定的事件序列中，一个 up 事件后面跟随另一个 up 事件的次数为两次。A 的中心右侧元素表示在事件序列中，一个 flat 事件后面跟随一个 down 事件的次数为一次。依此类推。
我们可以通过行或列标准化矩阵 A 来得到转移矩阵。如果我们使用基于行的标准化，则元素 M(i, j) 将描述在时间步 t 给定事件 E(i) 的情况下，在时间步 t+1 看到事件 E(j) 的概率。因此，每一行中的概率应加起来等于 1。在我们的示例中，行标准化矩阵如下所示：
类似地，如果我们使用基于列的标准化，则元素 M(i, j) 将描述在时间步 t 给定事件 E(j) 的情况下，在时间步 t-1 发生事件 E(i) 的概率。现在每一列中的概率应加起来等于 1。在我们的示例中，列标准化矩阵如下所示：
请注意，行标准化的条件概率（名义上是向前看）可能与列标准化的条件概率（向后看）不同。
Python 代码示例
要尝试这些概念，这里有一些基本的 Python 代码，它捕捉了上述示例中的内容。
首先确保你已经安装了 Pandas 包：
pip install pandas==0.25.2
然后运行以下代码：
import pandas as pd
# Define helper functions
def get_transition_tuples(ls):
''' Converts a time series into a list of transition tuples
'''
return [(ls[i-1], ls[i]) for i in range(1, len(ls))]
def get_transition_event(tup):
''' Converts a tuple into a discrete transition event
'''
transition_event = 'flat'
if tup[0] < tup[1]:
transition_event = 'up'
if tup[0] > tup[1]:
transition_event = 'down'
return transition_event
# Generate raw time series data
ls_raw_time_series = [1, 2, -2, -1, 0, 0, 2, 2, 1, 2, 3]
# Derive single-step state transition tuples
ls_transitions = get_transition_tuples(ls_raw_time_series)
# Convert raw time series data into discrete events
ls_events = [get_transition_event(tup) for tup in ls_transitions]
ls_event_transitions = get_transition_tuples(ls_events)
# Create an index (list) of unique event types
ls_index = ['up', 'flat', 'down']
# Initialize Markov transition matrix with zeros
df = pd.DataFrame(0, index=ls_index, columns=ls_index)
# Derive transition matrix (or co-occurrence matrix)
for i, j in ls_event_transitions:
df[j][i] += 1 # Update j-th column and i-th row
''' Derive row-normalized transition matrix:
- Elements are normalized by row sum (fill NAs/NaNs with 0s)
- df.sum(axis=1) sums up each row, df.div(..., axis=0) then divides each column element
'''
df_rnorm = df.div(df.sum(axis=1), axis=0).fillna(0.00)
''' Derive column-normalized transition matrix:
- Elements are normalized by column sum (fill NAs/NaNs with 0s)
- df.sum(axis=0) sums up each col, df.div(..., axis=1) then divides each row element
'''
df_cnorm = df.div(df.sum(axis=0), axis=1).fillna(0.00)
这应该生成以下转移矩阵：
>>> df # Transition matrix with raw event co-occurrences up flat down
up 2 2 1
flat 1 0 1
down 2 0 0 >>> df_rnorm # Row-normalized transition matrix up flat down
up 0.4 0.4 0.2
flat 0.5 0.0 0.5
down 1.0 0.0 0.0 >>> df_cnorm # Column-normalized transition matrix up flat down
up 0.4 1.0 0.5
flat 0.2 0.0 0.5
down 0.4 0.0 0.0
一种很好的可视化转移矩阵的方法是将它们描绘为有向加权图，使用像 Graphviz 或 NetworkX 这样的绘图包。
我们将在这里使用 Graphviz，因此你需要安装该包以便跟随操作：
pip install graphviz==0.13.2
值得浏览简短而精炼的 官方安装指南，以确保你正确安装了该包，特别是对于可能需要执行一些额外安装步骤的 Windows 用户。
一旦 Graphviz 设置好了，创建一些用于绘图的辅助函数：
from graphviz import Digraph
# Define functions to visualize transition matrices as graphs
def get_df_edgelist(df, ls_index):
''' Derive an edge list with weight values
'''
edgelist = []
for i in ls_index:
for j in ls_index:
edgelist.append([i, j, df[j][i]])
return pd.DataFrame(edgelist, columns=['src', 'dst', 'weight'])
def edgelist_to_digraph(df_edgelist):
''' Convert an edge list into a weighted directed graph
'''
g = Digraph(format='jpeg')
g.attr(rankdir='LR', size='30')
g.attr('node', shape='circle')
nodelist = []
for _, row in df_edgelist.iterrows():
node1, node2, weight = [str(item) for item in row]
if node1 not in nodelist:
g.node(node1, **{'width': '1', 'height': '1'})
nodelist.append(node1)
if node2 not in nodelist:
g.node(node2, **{'width': '1', 'height': '1'})
nodelist.append(node2)
g.edge(node1, node2, label=weight)
return g
def render_graph(fname, df, ls_index):
''' Render a visual graph and saves it to disk
'''
df_edgelist = get_df_edgelist(df, ls_index)
g = edgelist_to_digraph(df_edgelist)
g.render(fname, view=True)
现在你可以生成每个转移矩阵。默认情况下，输出的图形将存储在你的工作目录中。
# Generate graph of transition matrix (raw co-occurrences)
render_graph('adjmat', df, ls_index)
# Generate graph of row-normalized transition matrix
render_graph('transmat_rnorm', df_rnorm, ls_index)
# Generate graph of column-normalized transition matrix
render_graph('transmat_cnorm', df_cnorm, ls_index)
原始共现：
行标准化的转移概率：
列标准化的转移概率：
实际应用
描述性见解
我们可以用转移矩阵做的第一件最明显的事情是通过检查矩阵及其可视化图形表示来获得描述性见解。例如，利用我们前一节的示例输出，我们可以得到如下高层次的见解：
-
在 9 种可能的事件转移中，有 3 种在我们的样本中从未发生过（平稳 → 平稳、下降 → 下降 和 下降 → 平稳）。连续的平稳事件的低概率可能表明时间序列数据所跟踪的系统中的波动性。
-
上升 事件是唯一一个有非零概率（0.4）连续发生的事件类型。事实上，这个转移概率在我们的数据中是最高的之一，可能指向数据背后的系统中的强化效应。
-
在我们的案例中，基于行和列的归一化产生了不同的矩阵，尽管有些重叠。这告诉我们我们的时间序列在时间上本质上是不对称的，即我们观察到的模式在回顾或前瞻时会有所不同。
预测与回溯
通过将转移矩阵的副本连接起来，我们可以生成未来和过去的事件发生概率；这可以分别称为 预测 和 回溯。这里的一个核心假设是“历史无关紧要”；无论我们选择哪个时间步 t 作为参考点，我们都假设转移矩阵为 t+1（如果是行归一化）和 t-1（如果是列归一化）提供相关的概率。结果是，我们可以使用转移矩阵从任何任意时间步进行预测和回溯。特别地，我们可以使用行归一化的转移矩阵进行预测，使用列归一化的转移矩阵进行回溯。
以上述示例中计算的矩阵为例，假设我们在时间步 t = 25 观察到一个 上升 事件，我们希望预测在时间步 t = 27 时最可能发生的事件是什么。通过检查行归一化转移矩阵的顶行，我们可以直接看到，在紧接着的时间步 t = 26，观察到 上升、平稳 和 下降 事件的概率分别为 0.4、0.4 和 0.2。为了推导时间步 t = 27（即距离我们参考点两个时间步）的类似事件概率，我们需要将转移矩阵自身相乘，如下所示：
注意观察事件概率相对于我们的参考时间步长的变化。例如，假设在t = 25 时出现一个上升事件，那么在t = 26 时观察到另一个上升事件的概率为 0.4（向未来进一步一步），并且在t = 27 时增加到 0.56（向未来两步）。与此同时，在t = 26 时观察到平稳事件的概率也为 0.4，但在t = 27 时减少到 0.16。关键是，矩阵乘法支持预测和回溯。通常来说，向前或向后推测事件发生的概率n步，我们可以分别计算该转移矩阵的 n 次幂的行规范化或列规范化。
转移矩阵也可以用于预测原始的基础时间序列数据。假设上升或下降事件对时间序列数据造成了一单位的变化。现在假设时间序列从 1 到 2（一个上升事件）在t = 25 时发生，并且我们希望预测t = 26 和t = 27 时时间序列的进展。在上升事件之后，上升和平稳事件在t = 26 时都有最高的发生概率（0.4）。因此，我们可以预测，在t = 26 时，时间序列很可能是[1, 2, 3]或[1, 2, 2]，两者都有可能在t = 27 时产生两种更可能的情况：[1, 2, 3]可能变成[1, 2, 3, 4]或[1, 2, 3, 3]（概率均为 0.4），而[1, 2, 2]可能变成[1, 2, 2, 3]或[1, 2, 2, 1]（概率均为 0.5）。一般来说，我们期望所用于生成转移矩阵的数据集越大越丰富，能够捕捉更多事件链，因此每步的预测精度越高。
转移矩阵的乘法链导致了越来越复杂但完全可分解的原始事件转移概率的组合。这种可分解性可以帮助我们更深入地了解构成时间序列数据（或随机过程）的事件相互之间的关系。
收敛分析
连接转移矩阵的概念自然地引出了一个有趣的问题：转移矩阵 M 的概率会收敛吗？具体地说，是否存在一个稳定的转移矩阵M使得MM=M*？如果是，那么lim(n → ∞)Mⁿ = M，也就是说我们期望由矩阵乘法链Mⁿ表示的马尔可夫过程在某个时间点收敛到一个稳定状态M；在这种情况下，该过程是收敛的并因此稳定。假设我们的转移矩阵是行规范化的，元素M(i, j)给出了事件i之后紧随事件j的稳定长期概率。然而，如果找不到稳定矩阵M，那么该过程就不是收敛的也不是稳定的。
利用之前部分的运行示例，我们可以简要描述一下马尔可夫过程的收敛是如何通过分析得出的。
首先，让我们假设存在一个稳定的转移矩阵M，使得MM=M，并且M是行归一化的。由于我们知道M的样子，我们可以把矩阵乘法写成如下形式：
然后我们有下面的线性方程组：
如果这个方程组存在解（我们可以通过高斯消元等方法来检查），那么我们也可以得到一个收敛和稳定的转移矩阵。
包裹
一旦你掌握了这个技巧，使用马尔可夫转移矩阵来重构时间序列数据可以成为你数据科学工具箱中的一个有用部分。就像你通常用折线图来可视化时间序列数据以更好地把握整体趋势一样，转移矩阵提供了数据的一个互补表述，它高度压缩但在使用场景上非常多样化。当将其可视化为有向图时，转移矩阵已经可以用于获取高层次的描述性见解。当嵌入到更大的工作流程中时，转移矩阵可以成为预测和逆向预测更复杂方法的基础。此外，虽然我们在上述部分运行的简单例子将转移矩阵视为静态实体，但我们可以针对不同的时间间隔推导出不同的矩阵；这在分析显示数据中体现出明显的趋势逆转，反映在数据中具有显著 U 形或弯头形式模式的时间序列数据中尤为有用。显然，上面讨论的观点有几种可能的延伸，所以继续尝试吧——它们在你下一个数据科学项目中可能派上用场。
气候变化的时间序列：预测能源需求
原文：
towardsdatascience.com/time-series-for-climate-change-forecasting-energy-demand-79f39c24c85e
如何利用时间序列分析和预测来应对气候变化
·发布于 Towards Data Science ·阅读时间 7 分钟·2023 年 5 月 2 日
–
图片由 Matthew Henry 提供，来源于 Unsplash
这是系列文章 气候变化的时间序列 的第四部分。文章列表：
-
第一部分: 风能预测
-
第二部分: 太阳辐射预测
-
第三部分: 预测大型海浪
到目前为止，我们探讨了预测在将清洁能源源整合到电网中的重要性。
预测在能源系统的需求侧也发挥着关键作用。
平衡能源需求和供应
电力系统需要确保能源供应和需求之间的平衡。这一平衡对电网的可靠性至关重要。如果需求大于供应，会导致停电。当供应超过需求时，会出现多余的能源，往往被浪费掉。
电力系统使用预测模型来帮助预测能源需求。准确的需求预测有助于更高效地生产和使用能源。这直接影响到气候，因为它减少了浪费。
分析能源消耗在家庭中也很有价值。例如，个人可以检查哪些电器消耗更多的能源，并利用这些信息在高峰时段避免更高的费用。顺便提一句：据估计，约 8%的住宅电力需求来自待机功耗 [4]。
预测能源需求
预测能源需求是一个困难的问题。
能源消耗依赖于多个因素，其中一些可能不容易获得用于建模。例如，包括影响加热或制冷设备使用的天气和经济状况。天气的特征是高度变化的模式。这些模式使得预测天气对能源需求的影响程度变得困难。
能源需求数据在不同时间尺度上展示季节性模式，包括每日、每月或每年。例如，在冬季，能源需求因加热需求而增加。
按月分布的能源需求。图片由作者提供。
捕捉所有季节性效应在构建准确的能源需求预测模型中非常重要。
你怎么做呢？
实操：预测能源需求
在本文的其余部分，我们将开发一个能源需求预测模型。你将学到如何：
-
使用自相关和数据可视化分析多个季节性效应；
-
从日期和时间信息中提取特征，以处理多个季节性效应。
本教程中使用的完整代码可在 Github 上找到：
数据集
我们将使用一个代表美国肯塔基州每小时电力消耗（以兆瓦为单位）的数据集[1]。数据从 2013 年到 2018 年收集，总共有 45,344 个观测值。
这就是数据的样子：
美国肯塔基州的每小时能源需求时间序列。黄色线是每日平均需求。数据来源见参考文献 [1]。图片由作者提供。
可视化季节性
上述图表表明系列中存在年度规律模式。
另一种可视化季节性的方法是使用名为季节性图的图形：
每年的平均月能源需求的季节性图。图片由作者提供。
季节性图使得查看年度月度模式变得容易。例如，能源消耗在冬季和夏季增加，在春季和秋季减少。这可能与加热（冬季）或制冷（夏季）有关。
下图中的季节性子序列图对于分析数据在月份内部和跨月份的动态也很有帮助：
每年的平均月能源需求的季节性子序列图。图片由作者提供。
使用自相关函数（ACF）分析季节性
你还可以使用自相关来分析季节性。季节性时间序列在每个季节性滞后中会显示出更高的自相关性。
自相关图。波动表明每日季节性。图片由作者提供。
上图展示了高达 48 个滞后的自相关图。自相关显示了一个由每日季节性引起的振荡模式。
所以，除了年度规律性变化外，还有明显的每日季节性。特定小时观察到的值与前一天同一小时捕捉到的值相关。
自回归
作为起点，我们将开发一个用于预测电力消耗的自回归模型。你可以查看之前的文章了解这种建模的详细信息。
from sklearn.model_selection import train_test_split
from src.tde import time_delay_embedding
# Train / test split
train, test = train_test_split(series, test_size=0.2, shuffle=False)
# using past 12 observations as explanatory variables
N_LAGS = 12
# using the next 12 hours as the forecasting horizon
HORIZON = 12
# transforming time series into a tabular format for supervised learning
X_train, Y_train = time_delay_embedding(train, n_lags=N_LAGS, horizon=HORIZON, return_Xy=True)
X_test, Y_test = time_delay_embedding(test, n_lags=N_LAGS, horizon=HORIZON, return_Xy=True)
处理多重季节效应
大多数时间序列分解或预测方法设计为处理单一季节周期。
随着数据采集成本的降低，它使得以高采样频率（如每日或每小时）收集时间序列成为可能。高频时间序列提供了更多的数据，这是训练机器学习模型的重要因素。然而，这些数据也包含复杂的季节性模式，这可能很难建模。
那么，你如何处理多重季节效应？
我们将应用两个特征提取过程来完成这项工作：
-
基于日期和时间数据提取特征；
-
使用基于傅里叶项的三角函数表示。
基于日期和时间提取特征
建模多重季节模式的一个简单方法是总结每个时间步的日期和时间信息。你可以提取相关的周期，比如小时或一年中的天数，并将其作为解释变量。
库 sktime 提供了一个方便的类来实现这一点：
from sktime.transformations.series.date import DateTimeFeatures
hourly_feats = DateTimeFeatures(ts_freq='H',
keep_original_columns=False,
feature_scope='efficient')
dtime_train = hourly_feats.fit_transform(X_train)
这是输出的一个示例：
这些特征对于许多具有复杂季节性的时间序列可能是有效的方法。但是，它们可能会忽略时间的连续性。例如，假设你将月份信息编码为从 1（1 月）到 12（12 月）的整数：
将月份信息编码为整数。图像由作者提供。
在这种情况下，模型无法理解 12 月在 1 月之前。使用独热编码转换月份信息也会导致这个问题。
傅里叶项
傅里叶项是基于正弦和余弦波的周期性和确定性系列。这些项的平滑性使得时间连续性的建模成为可能。
使用傅里叶项编码月份信息。图像由作者提供。
这是一个用于提取傅里叶级数的 sklearn 兼容类：
from datetime import datetime
import numpy as np
import pandas as pd
class FourierTerms:
def __init__(self, period: float, n_terms: int, prefix=''):
self.period = period
self.n_terms = n_terms
self.prefix = prefix
def transform(self, index: pd.DatetimeIndex, use_as_index: bool = True):
t = np.array(
(index - datetime(1970, 1, 1)).total_seconds().astype(float)
) / (3600 * 24.)
fourier_x = np.column_stack([
fun((2.0 * (i + 1) * np.pi * t / self.period))
for i in range(self.n_terms)
for fun in (np.sin, np.cos)
])
col_names = [
f'{self.prefix}{fun.__name__[0].upper()}{i}'
for i in range(self.n_terms)
for fun in (np.sin, np.cos)
]
fourier_df = pd.DataFrame(fourier_x, columns=col_names)
if use_as_index:
fourier_df.index = index
return fourier_df
主要输入是季节周期（例如，12 用于每月时间序列）、日期时间信息和项数。项数影响表示的平滑性。最佳项数取决于输入数据。
几种流行的方法使用傅里叶项来建模复杂的季节性。这些方法包括Prophet、TBATS、greykite或时间变化回归。
傅里叶项作为解释性特征包括在内，如下所示：
fourier_daily = FourierTerms(n_terms=2, period=24, prefix='D_')
fourier_monthly = FourierTerms(n_terms=2, period=24 * 30.5, prefix='M_')
fourier_yearly = FourierTerms(n_terms=2, period=24 * 365, prefix='Y_')
dfourier_train = fourier_daily.transform(X_train.index)
mfourier_train = fourier_monthly.transform(X_train.index)
yfourier_train = fourier_yearly.transform(X_train.index)
feats_train = pd.concat([X_train, dtime_train, dfourier_train,
mfourier_train, yfourier_train],
axis=1)
model = RandomForestRegressor()
model.fit(feats_train, Y_train)
因此，自动回归与总结日期和时间数据的特征相结合。下面的特征重要性评分表明这些特征为模型提供了相关信息：
随机森林重要性评分。图片由作者提供。
就像预测能源生产的情况一样，能源需求模型在较长时间范围内的准确性会降低：
MASE 误差按预测时间范围。图片由作者提供。
主要结论
-
能源需求预测与气候变化相关。它们使电力系统能够做出明智的决策，并将清洁能源来源整合到电网中；
-
需求时间序列受到多个因素的影响，并具有复杂的季节性；
-
你可以通过总结日期和时间信息来处理多个季节性模式。傅里叶项是实现这一目标的常用方法；
-
能源需求在长期内难以预测，即超过几小时。提高长期预测的准确性对电力系统的效率至关重要。
感谢阅读，下次故事再见！
参考文献
[1] PJM 每小时能源消耗（许可证：CC0：公共领域）
[2] Rolnick, David, 等。“利用机器学习应对气候变化。” ACM 计算调查（CSUR）55.2（2022）：1–96。
[3] MacKay, David. 可持续能源-没有热空气。UIT cambridge，2008 年。
气候变化时间序列：大型海洋波浪预测
原文：
towardsdatascience.com/time-series-for-climate-change-forecasting-large-ocean-waves-78484536be36
如何使用时间序列分析和预测来应对气候变化
·发表在 Towards Data Science ·7 分钟阅读·2023 年 4 月 25 日
–
照片由 Silas Baisch 提供，来源于 Unsplash
这是系列气候变化时间序列的第三部分。文章列表：
-
第一部分: 风能预测
-
第二部分: 太阳辐射预测
海洋波浪能
波浪能转换器是从海洋波浪中获取能量的浮标。照片（浮标的照片，但不是 WEC 的照片）由 Emmy C 提供，来源于 Unsplash
海洋波浪是一个有前景的可再生能源来源。
为什么选择海洋波浪？
当谈到可再生能源时，人们通常会想到太阳能或风能。这些是最受欢迎的可再生能源来源。然而，由于其一致性，海洋波浪具有巨大潜力。
我们可以在约 90%的时间内从海洋波浪中获取能量。这个数字对于太阳能或风能来说大约是 20%-30%。详细信息请参见参考文献[1]。例如，太阳能技术仅在白天有效。
除了其一致性，波浪能也比上述两种替代能源更具预测性。海洋波浪能的主要限制是生产成本。目前，这些成本相对于太阳能或风能较高。
从波浪到电力
波浪能通过称为波浪能转换器（WECs）的设备转换为电力。这些设备是放置在海面上的浮标。
海洋波浪输出的能量取决于波浪的高度。这个量随时间变化。因此，预测波浪的高度是从海洋波浪中高效生产能量的关键任务。
波浪的高度是根据所谓的显著波高来量化的。该量值定义为最高三分之一波浪的平均波高，从波谷到波峰：
统计波浪分布。图片来源here（公有领域许可证）。
当波浪过大时——第二个动机
预测波浪高度对于除能源生产之外的其他因素也很重要。
预测即将到来的大浪有助于管理海洋操作的安全。准确的预测可以防止沿海灾害，并保护波能转换器。这些设备可能需要关闭以防止大浪造成损坏。
安全问题也与船只通行有关。一艘船只在移动时需要最低水深。大浪会减少水深，这可能无法达到最低要求。因此，预测可以提高船只移动的效率，从而降低成本，提高港口的可靠性。
总结来说，预测海洋波浪的高度很重要，原因有几个：
-
估计能源生产以管理电网；
-
管理海洋操作，包括船只通行。
教程：预测大浪
在本文的其余部分，我们将开发一个模型来预测即将到来的大浪。你将学习如何构建一个概率预测模型，该模型估计事件发生的可能性。
完整的代码可在 Github 上获得：
数据集
我们将使用一个来自爱尔兰海岸放置的智能浮标收集的实际数据集。除了其他信息，收集的数据包括显著波高，这是我们想要预测的变量。此数据在参考文献[2]的来源中提供。
下载数据后，你可以使用以下代码读取：
import pandas as pd
file = 'path_to_data/IrishSmartBuoy.csv'
# reading the data set
# skipping the second with skiprows
# parsing time column to datetime and setting it as index
data = pd.read_csv(file, skiprows=[1], parse_dates=['time'], index_col='time')
# defining the series and converting cm to meters
series = data['SignificantWaveHeight'] / 100
# resampling to hourly and taking the mean
series = series.resample('H').mean()
以下是数据的样子：
每小时显著波高时间序列。图片来源于作者
有几个时期的数据缺失。也许浮标正在维护。此外，还可以看到明显的季节性因素。波浪的高度在冬季通常较高。
时间序列的分布。图片来源于作者
数据分布显示出右偏态。重尾表示大浪，这对预测很重要。
关于如何定义大浪没有明确的共识。我们将其定义为至少 5 米高的波浪。这个阈值在上面的图中用黄色垂直线表示。
因此，目标是预测未来的波浪是否会超过 5 米。这个问题可以被框架化为超越概率预测任务。
超越概率预测概述
超越概率预测是预测时间序列在预定义的未来时间段内超越预定义阈值的概率的问题。
这个任务在极端值高度相关的领域是重要的。例如，在环境科学中，模型用于预测自然灾害如飓风或洪水的可能性。另一个例子是工程领域，专业人员使用模型来预测设备故障的可能性。
从机器学习的角度来看，主要挑战是缺乏关于超越事件的信息。根据定义，超越事件是稀有的。模型需要应对不平衡的类别分布，大多数观察值是非事件情况。
为什么用概率？
大浪的发生很难预测。因此，将预测结果以概率的形式呈现是有用的，以传达这些预测背后的不确定性和局限性。
一般来说，概率预测通过改善对每个可能行动相关风险的评估来帮助决策。
建立模型
我们可以使用自回归来处理这个任务。
目标是使用过去近期的海洋波高作为解释变量。目标是一个二元变量，表示是否会很快出现大浪。
from sklearn.model_selection import train_test_split
# https://github.com/vcerqueira/tsa4climate/tree/main/src
from src.tde import time_delay_embedding
# using past 24 observations as explanatory variables
N_LAGS = 24
# using the next 12 hours as the forecasting horizon
HORIZON = 12
# forecasting the probability of waves above 5 meters
THRESHOLD = 5
# leaving last 20% of observations for testing
train, test = train_test_split(series, test_size=0.2, shuffle=False)
# transforming time series into a tabular format for supervised learning
X_train, Y_train = time_delay_embedding(train, n_lags=N_LAGS, horizon=HORIZON, return_Xy=True)
X_test, Y_test = time_delay_embedding(test, n_lags=N_LAGS, horizon=HORIZON, return_Xy=True)
y_train = Y_train.apply(lambda x: (x > THRESHOLD).any(), axis=1).astype(int)
y_test = Y_test.apply(lambda x: (x > THRESHOLD).any(), axis=1).astype(int)
在这种情况下，我们将预测范围设置为 12 小时。因此，在每个时刻，我们的目标是预测在接下来的 12 小时内出现大浪的可能性。
目标变量的分布不平衡，大多数观察值指的是正常的波浪：
每个数据分区中按波浪类型（大/正常）的分布。图片作者提供。
由于目标变量是二元的，我们希望构建一个二元概率分类模型。对于这个案例研究，我们选择了随机森林。
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import roc_auc_score, roc_curve
model = RandomForestClassifier(max_depth=5)
model.fit(X_train, y_train)
probs = model.predict_proba(X_test)[:, 1]
roc_auc_score(y_test, probs)
fpr, tpr, thresholds = roc_curve(y_test, probs)
在训练模型后，我们可以使用*.predict_proba*方法获得概率预测。
二元概率预测可以使用 ROC（接收器操作特征）曲线来评估。其思想是绘制假阳性率（x 轴）与真正率（y 轴）在不同决策阈值下的关系。
这会生成如下的曲线：
预测模型的 ROC 曲线。图片作者提供。
曲线越接近左上角越好。对角线虚线是随机模型应该达到的水平。
ROC 曲线通常由其下方的面积（AUC）来总结。AUC 是用于评估二元概率分类器的指标，它量化了模型区分两类的能力。我们模型的 AUC 为 0.94，这是一个不错的分数。
关键要点
-
海洋波浪是一个有前景的可再生能源来源；
-
虽然这一来源比太阳能或风能更为稳定，但生产成本却阻碍了它的推广；
-
大型海洋波浪对海事操作的安全构成了担忧，包括沿海灾害或船只的通行。此外，大波浪还可能损坏波浪能转换器；
-
预测即将到来的大波浪是一项有用的任务。概率预测更为理想，因为这些提供了对决策至关重要的更多信息；
-
预测模型的进步可以帮助加速海洋波浪能的采纳。
感谢阅读，下次故事见！
参考文献
[1] Drew, Benjamin, Andrew R. Plummer, and M. Necip Sahinkaya. “波浪能转换器技术综述。” (2009): 887–902。
[2] 爱尔兰波浪浮标（许可证：创意共享署名 4.0）
气候变化中的时间序列：风力发电预测
原文：
towardsdatascience.com/time-series-for-climate-change-forecasting-wind-power-8ed6d653a255
如何利用时间序列分析和预测应对气候变化
·发布于 Towards Data Science ·阅读时间 6 分钟·2023 年 3 月 29 日
–
由 American Public Power Association 提供的照片，来源于 Unsplash
朝着清洁能源生产迈进
非可再生能源对我们的星球造成了沉重的生态足迹。这一问题促使了清洁能源的科学和技术进步，例如太阳能、风能和海洋波浪能。这些能源对环境友好，不像煤炭或石油。
延迟清洁能源广泛应用的原因之一是它们的不规则性。它们是高度可变的资源，这使得其行为难以预测。
因此，预测这些资源的条件是一个关键挑战。准确的预测对于高效生产清洁能源至关重要。
在本文中，我们将开发一个模型来预测风力发电。
风力发电
风力发电是日益普及的可再生能源之一。截至 2020 年，风力发电约占丹麦电力生产的 47%。其他国家也增加了电力网中的风力发电份额。
风力发电也有一些缺点。例如，风力涡轮机的视觉影响和噪音。此外，风力发电基础设施需要相当大的初期投资。
风力发电的电网整合也很困难。风力发电只能在风吹起时生成。这使得它成为一种间歇性和不可预测的能源。因此，它需要与其他替代能源配合使用。
风力突增
风力波动也是电力系统运营商面临的主要问题。这些是在短时间内（分钟到小时）风力发电的大幅变化。如果未及时检测，风力波动可能会影响电网的可靠性。
风力波动可以是向上或向下的变化。当发生突然的功率下降时，必须提高其他来源的能量来补偿损失。突然的上升变化可能促使操作员减少其他来源的输出，或者选择销售过剩的能源。
预测的作用
电力系统运营商依赖预测模型来预测风力条件。这些模型使运营商能够高效地平衡和整合多个能源来源。准确的预测对于电网的效率以及降低成本非常重要。
实践操作
在本文的其余部分，我们将构建一个预测风力发电的模型。目标是展示这个问题的挑战性以及未来的发展如何带来价值。
你可以在 Github 上找到这个项目的完整代码：
数据集
在这个教程中，我们将使用一个关于比利时风电场的公开数据集。
时间序列以 15 分钟为间隔从 2014 年到 2018 年收集。除了风力发电，我们还获取了关于已安装容量的信息（最大可能生成的电力）：
风力发电时间序列（兆瓦），包括已安装的容量。数据来源见参考文献[1]。图片由作者提供。
随着新风机的加入，已安装容量随时间增加。因此，我们将风力发电量标准化为风电场的容量。这就形成了一种风力发电量占总容量百分比的度量。
风力发电占已安装容量的百分比。前四个月的样本。图片由作者提供。
构建预测模型
我们将使用机器学习算法构建一个预测模型。这个想法是应用一种叫做自回归的建模技术。自回归涉及使用最近的观察数据（滞后）来预测未来的观察数据。你可以在以前的文章中阅读更多关于自回归的内容。
首先，我们需要将时间序列转换为表格格式。这可以通过一种叫做时间延迟嵌入的滑动窗口方法完成：
from sklearn.model_selection import train_test_split
# src module here: https://github.com/vcerqueira/tsa4climate/tree/main/src
from src.tde import time_delay_embedding
# number of lags and forecasting horizon
N_LAGS, HORIZON = 24, 24
# leaving last 20% of observations for testing
train, test = train_test_split(series, test_size=0.2, shuffle=False)
# transforming time series into a tabular format for supervised learning
X_train, Y_train = time_delay_embedding(train, n_lags=N_LAGS, horizon=HORIZON, return_Xy=True)
X_test, Y_test = time_delay_embedding(test, n_lags=N_LAGS, horizon=HORIZON, return_Xy=True)
我们将滞后数和预测视野都设置为 24。在每个时间步骤中，我们希望模型使用过去 24 小时的数据来预测接下来 24 小时的波动高度。
这是训练解释变量和目标变量的样本：
第一和最后的解释变量与目标变量的样本
接下来，我们使用训练集选择一个模型。在本教程中，我们进行随机搜索以选择和优化回归算法。此外，我们还测试是否应包含特征提取。特定的特征提取过程基于滞后汇总统计。
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import FunctionTransformer
from sklearn.model_selection import RandomizedSearchCV, TimeSeriesSplit
from src.model_selection import (MetaEstimator,
search_space_with_feature_ext)
# Create a pipeline for hyperparameter optimization
# 'feature' contains different possibilities for feature extraction
# 'model' contains different regression algorithms and respective hyperparameters
pipeline = Pipeline([('feature', FunctionTransformer()),
('model', MetaEstimator())])
# do random search optimization for model selection
search_mod = RandomizedSearchCV(estimator=pipeline,
param_distributions=search_space_with_feature_ext,
scoring='r2',
n_iter=30,
n_jobs=1,
refit=True,
verbose=2,
cv=TimeSeriesSplit(n_splits=3),
random_state=123)
search_mod.fit(X_train, Y_train)
print(search_mod.best_estimator_)
# Pipeline(steps=[('feature', FunctionTransformer()),
# ('model', RidgeCV(alphas=0.25))])
选择的模型是没有特征提取的岭回归。我们使用完整的训练集重新训练该模型，并在测试集上进行评估。
from sklearn.metrics import r2_score
# forecasting testing observations using the selected model
Y_hat_test = search_mod.predict(X_test)
Y_hat_test = pd.DataFrame(Y_hat_test, columns=Y_train.columns)
# evaluating the selected model over the forecasting horizon
r2_scores = {col: r2_score(y_true=Y_test[col], y_pred=Y_hat_test[col])
for col in Y_hat_test}
结果揭示了两个主要挑战：
-
在大范围内进行预测；
-
预测极端观察值。
预测范围内的不确定性增加
众所周知，长期预测（即超过几个时间步）比短期预测更困难。风力预测也遵循这一趋势：
在预测范围内的预测性能（R²得分）。图片由作者提供。
预测下一小时的波浪功率是容易的。R2 得分约为 0.98——几乎完美。然而，当我们进行较长时间范围的预测时，性能会大幅下降。
长期预测（在这种情况下，超过几个小时）对有效平衡能源需求和供应非常重要。
预测极端值
显示实际值和预测值关系的散点图。图片由作者提供。
我们还需要模型在预测极端值时的准确性。在这种情况下，这些极端值代表高或低风力。这样的值很重要，因为它们可能对电网操作产生影响。
上图显示了一个预测值和实际值的散点图。红色虚线是理想情况，其中预测值与观察值相符。
大部分数据点围绕红线分布。但是，对于极端观察值，数据点偏离了红线。因此，模型有时未能预测极端值。
克服这两个局限性对于将风能整合到电力网中至关重要。
主要收获
-
风力发电是越来越受欢迎的可再生能源；
-
预测风力条件是估算来自此来源的能量的关键任务；
-
这个预测任务可以通过基于自回归的模型来解决。
-
预测风力具有两个挑战：低长期预测性能和极端观察值的低性能。
感谢阅读。这篇文章是关于气候变化时间序列系列帖子的第一篇。敬请关注更多内容！
参考文献
[1] 比利时风力发电风电场（许可证 CC BY 4.0）
[2] Rolnick, David, 等. “利用机器学习应对气候变化。” ACM 计算机调查 (CSUR) 55.2 (2022): 1–96.
气候变化时间序列：起源-目的地需求预测
挖掘浮动车数据以应对气候变化
·发表于 Towards Data Science ·阅读时间 6 分钟·2023 年 6 月 15 日
–
图片来源: Denys Nevozhai 于 Unsplash
这是《气候变化的时间序列》系列的第八部分。文章列表：
-
第一部分: 预测风力
-
第二部分: 太阳辐射预测
-
第三部分: 预测大洋波浪
-
第四部分: 预测能源需求
-
第五部分: 预测极端天气事件
-
第六部分: 使用深度学习进行精准农业
-
第七部分: 通过聚类减少食品浪费
浮动车数据用于流动建模
挖掘浮动车数据是智能交通系统中的关键任务。浮动车数据指的是由配备 GPS 设备的车辆收集的数据，这些数据提供了关于车辆位置和速度的信息。
了解城市内部的流动模式是交通运输中的一项重要任务。例如，这有助于减少拥堵和整体交通活动。减少在交通中的时间意味着排放的温室气体减少。因此，准确的模型对气候变化有积极影响。
GPS 设备的普及产生了许多与移动性相关的数据集。但是，从 GPS 数据中学习是一个具有挑战性的问题。空间依赖性很复杂但非常重要。而且，还有时间依赖性，例如高峰时段。移动模式在工作日和非工作日也有所不同。
起始-目的地流量计数估计
旧金山的一次出租车行程。图片来源：作者
浮动车数据为移动建模提供了许多可能性。其中之一是起始-目的地（OD）流量计数问题。
OD 流量计数指的是在给定时间段内估计有多少辆车辆从一个子区域穿越到另一个子区域。这项任务有多个相关的原因。出租车公司可以根据特定区域的预期需求动态分配车队。
实操：预测旧金山的 OD 需求
在本文的其余部分，我们将预测旧金山的出租车乘客需求。我们将把这个问题作为 OD 流量计数任务来处理。
本教程中使用的完整代码可以在 Github 上找到：
数据集
我们将使用由美国加利福尼亚州旧金山的一家出租车车队收集的数据集。该数据集包含来自 536 辆出租车的 GPS 数据，时间跨度为 21 天。总共有 1.21 亿条 GPS 路径，分布在 464045 次行程中。有关更多细节，请参见参考文献 [1]。
数据集的样本。
在每个时间步长和每辆出租车上，我们都有其坐标信息和是否有乘客乘坐的信息。
问题定义
旧金山的一些出租车行程的结束位置。图片来源：作者。
我们的目标是根据起始点建模人们的去向。OD 流量计数估计可以分为四个子任务：
-
空间网格分解
-
起始-目的地对的选择
-
时间离散化
-
建模和预测
让我们逐个深入探讨每个问题。
空间网格分解
空间分解是 OD 流量计数估计的一个常见预处理步骤。其思想是将地图分割成网格单元，这些网格单元代表城市的一小部分。然后，我们可以统计每对网格单元之间有多少人穿越。
旧金山的两个示例网格单元。图片来源：作者。
在这个案例研究中，我们将城市地图分成 10000 个网格单元，如下所示：
import pandas as pd
from src.spatial import SpatialGridDecomposition, prune_coordinates
# reading the data set
trips_df = pd.read_csv('trips.csv', parse_dates=['time'])
# removing outliers from coordinates
trips_df = prune_coordinates(trips_df=trips_df, lhs_thr=0.01, rhs_thr=0.99)
# grid decomposition with 10000 cells
grid = SpatialGridDecomposition(n_cells=10000)
# setting bounding box
grid.set_bounding_box(lat=trips_df.latitude, lon=trips_df.longitude)
# grid decomposition
grid.grid_decomposition()
在上面的代码中，我们去除了异常地点。这些异常可能是由于 GPS 故障引起的。
获取最受欢迎的行程
在空间分解过程之后，我们获取了每次出租车行程的起点和终点，当它们被乘客占用时。
from src.spatial import ODFlowCounts
# getting origin and destination coordinates for each trip
df_group = trips_df.groupby(['cab', 'cab_trip_id'])
trip_points = df_group.apply(lambda x: ODFlowCounts.get_od_coordinates(x))
trip_points.reset_index(drop=True, inplace=True)
该想法是重建数据集，以包含以下信息：每次乘客旅行的起点、终点和起始时间戳。这些数据构成了我们的出发地-目的地 (OD) 流量计数模型的基础。
这些数据使我们能够计算从 A 单元格到 B 单元格的旅行次数：
# getting the origin and destination cell centroid
od_pairs = trip_points.apply(lambda x: ODFlowCounts.get_od_centroids(x, grid.centroid_df), axis=1)
为了简便起见，我们选择了旅行次数最多的前 50 个 OD 网格单元对。选择这个子集是可选的。然而，只有少量旅行的 OD 对在时间上会显示稀疏需求，这很难建模。此外，从车队管理的角度来看，需求低的旅行可能不太有用。
flow_count = od_pairs.value_counts().reset_index()
flow_count = flow_count.rename({0: 'count'}, axis=1)
top_od_pairs = flow_count.head(50)
时间离散化
在找到需求最高的 OD 对后，我们将其在时间上离散化。通过计算每小时每个给定顶级对的旅行次数来完成此操作。可以按如下方式进行：
# preparing data
trip_points = pd.concat([trip_points, od_pairs], axis=1)
trip_points = trip_points.sort_values('time_start')
trip_points.reset_index(drop=True, inplace=True)
# getting origin-destination cells for each trip, and origin start time
trip_starts = []
for i, pair in top_od_pairs.iterrows():
origin_match = trip_points['origin'] == pair['origin']
dest_match = trip_points['destination'] == pair['destination']
od_trip_df = trip_points.loc[origin_match & dest_match, :]
od_trip_df.loc[:, 'pair'] = i
trip_starts.append(od_trip_df[['time_start', 'time_end', 'pair']])
trip_starts_df = pd.concat(trip_starts, axis=0).reset_index(drop=True)
# more data processing
od_count_series = {}
for pair, data in trip_starts_df.groupby('pair'):
new_index = pd.date_range(
start=data.time_start.values[0],
end=data.time_end.values[-1],
freq='H',
tz='UTC'
)
od_trip_counts = pd.Series(0, index=new_index)
for _, r in data.iterrows():
dt = r['time_start'] - new_index
dt_secs = dt.total_seconds()
valid_idx = np.where(dt_secs >= 0)[0]
idx = valid_idx[dt_secs[valid_idx].argmin()]
od_trip_counts[new_index[idx]] += 1
od_count_series[pair] = od_trip_counts.resample('H').mean()
od_df = pd.DataFrame(od_count_series)
这会产生一组时间序列，每个顶级 OD 对一个时间序列。以下是四个示例对的时间序列图：
四个示例出发地-目的地对的流量计数时间序列。图片由作者提供。
时间序列显示出日常季节性，这主要由高峰时段驱动。
预测
由时间离散化产生的时间序列集可用于预测。我们可以构建一个模型来预测相对于给定 OD 对的乘客旅行需求量。
下面是如何针对一个示例 OD 对进行操作的步骤：
from pmdarima.arima import auto_arima
# getting the first OD pair as example
series = od_df[0].dropna()
# fitting an ARIMA model
model = auto_arima(y=series, m=24)
上述内容，我们基于 ARIMA 构建了一个预测模型。该模型根据最近的需求预测下一小时的乘客需求。我们使用 ARIMA 方法是为了简便，但也可以使用其他方法，如深度学习。
深入了解图神经网络
上述方法是一种简单而有效的解决 OD 流量计数问题的方法。但它将每个 OD 对视为一个独立的时间序列。
实际上，每对之间与相邻的 OD 对或周围道路相关联。因此，图神经网络在预测交通状况方面得到了越来越多的应用。道路网络被建模为图形，神经网络可以捕捉其中的复杂交互。你可以查看这个 Keras 示例来了解如何实现这种方法。
关键要点
-
移动建模是智能交通系统中的一项重要任务；
-
OD 流量计数模型可以帮助减少城市交通，从而减少温室气体的排放；
-
你可以通过基于空间分解和时间离散化的方法来解决 OD 流量计数问题。这将生成每个 OD 对的一组时间序列，可用于预测。
感谢阅读，下次故事见！
参考文献
[1] 旧金山，USA 的出租车移动轨迹数据集。（许可证 CC BY 4.0）
[2] Moreira-Matias, Luís 等人。“使用高速 GPS 数据流进行时变 OD 矩阵估计。” 应用专家系统 44 (2016): 275–288。
气候变化的时间序列: 通过聚类减少食物浪费
使用时间序列聚类进行更好的需求预测
·发表于 Towards Data Science ·6 分钟阅读·2023 年 6 月 7 日
–
这是系列气候变化的时间序列的第七部分。文章列表：
减少食物浪费
改善供应链是减少我们生态足迹的另一个关键步骤。在发达国家，通常会有大量的消费品，如食品。这些过剩的物品需要大量的能源和资源，通常会被浪费。
减少过度生产是减少温室气体排放的重要里程碑。我们可以通过更好地了解需求来解决这个问题。
以食品为例。每年我们损失约 13 亿公吨食品[1]。当然，这并非所有都是剩余食品或与供应链有关的。一部分食品在生产或运输过程中丧失，例如由于冷藏条件差。然而，更好的需求预测模型可以对减少过度生产产生显著影响。
食品需求时间序列的聚类
图片由Diego Marín提供，来源于Unsplash
我们可以使用聚类分析来改进需求预测。
聚类涉及根据相似性对观察值进行分组。在这种情况下，每个观察值是代表某个产品销售的时间序列。一般而言，您可以使用时间序列聚类来：
-
识别具有相似模式的时间序列，例如趋势或季节性；
-
将时间序列划分为不同的组。这在时间序列数量较大时特别有用。
对于需求时间序列，聚类可以用来识别销售模式相似的产品。然后可以根据每个簇的特征量身定制预测模型。最终，这将带来更好的预测。
聚类需求时间序列对企业也很有价值。识别相似的产品有助于制定更好的营销或促销策略。
实践
在本文的其余部分，我们将对食品需求时间序列进行聚类分析。您将学习如何：
-
使用特征提取总结一组时间序列；
-
使用 K-Means 和层次方法进行时间序列聚类。
完整代码可在 Github 上获取：
数据集
我们将使用美国农业部收集的每周食品销售时间序列数据集。该数据集包含按产品类别和子类别划分的食品销售信息。时间序列按州划分，但我们将使用每个时期的全国总销售额。
下面是数据集的一个示例：
美国不同产品子类别的销售额（以百万美元计）
整体数据如下所示：
不同食品子类别的销售金额（百万美元）。图像由作者提供。
基于特征的时间序列聚类
我们将使用基于特征的时间序列聚类方法。这个过程包括两个主要步骤：
-
将每个时间序列总结为一组特征，例如平均值；
-
对特征集应用传统的聚类算法，例如 K-means。
我们逐步进行每个步骤。
使用tsfel进行特征提取
我们从提取一组统计数据开始，以总结每个时间序列。目标是将每个系列转换为一小组特征。
有几种时间序列特征提取工具。我们将使用tsfel，它相对于其他方法提供了竞争力的性能[3]。
这是如何使用tsfel的方法：
import pandas as pd
import tsfel
# get configuration
cfg = tsfel.get_features_by_domain()
# extract features for each food subcategory
features = {col: tsfel.time_series_features_extractor(cfg, data[col])
for col in data}
features_df = pd.concat(features, axis=0)
这个过程会生成大量特征。其中一些可能是冗余的，因此我们进行特征选择过程。
下面，我们对特征集进行三项操作：
-
归一化：将变量转换为 0–1 的值范围；
-
方差选择：去除方差为 0 的变量；
-
相关性选择：去除与其他现有变量有高相关性的变量。
from sklearn.preprocessing import MinMaxScaler
from sklearn.feature_selection import VarianceThreshold
from src.correlation_filter import correlation_filter
# normalizing the features
features_norm_df = pd.DataFrame(MinMaxScaler().fit_transform(features_df),
columns=features_df.columns)
# removing features with 0 variance
min_var = VarianceThreshold(threshold=0)
min_var.fit(features_norm_df)
features_norm_df = pd.DataFrame(min_var.transform(features_norm_df),
columns=min_var.get_feature_names_out())
# removing correlated features
features_norm_df = correlation_filter(features_norm_df, 0.9)
features_norm_df.index = data.columns
K 均值聚类
在预处理数据集后，我们准备对时间序列进行聚类。我们将每个序列总结为一小组无序特征。因此，我们可以使用任何传统的聚类算法。一个流行的选择是 K 均值。
使用 K 均值时，我们需要选择所需的聚类数量。除非我们有一些领域知识，否则这个参数没有明显的先验值。但是，我们可以采用数据驱动的方法来选择聚类数量。我们测试不同的值，并选择最佳值。
下面，我们测试最多 24 个聚类的 K 均值算法。然后，我们选择最大化轮廓系数的聚类数量。这个指标量化了获得的聚类的凝聚度。
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
kmeans_parameters = {
'init': 'k-means++',
'n_init': 100,
'max_iter': 50,
}
n_clusters = range(2, 25)
silhouette_coef = []
for k in n_clusters:
kmeans = KMeans(n_clusters=k, **kmeans_parameters)
kmeans.fit(features_norm_df)
score = silhouette_score(features_norm_df, kmeans.labels_)
silhouette_coef.append(score)
如下图所示，轮廓系数在 5 个聚类时最大。
最多 24 个聚类的轮廓系数。图片由作者提供。
我们可以绘制平行坐标图来理解每个聚类的特征。以下是一个包含三个特征的样本示例：
带有特征样本的平行坐标图。图片由作者提供。
我们还可以利用聚类的信息来改进需求预测模型。例如，通过为每个聚类构建一个模型。参考文献[5]中的论文是这种方法的一个很好的例子。
层次聚类
层次聚类是 K 均值的替代方法。它通过迭代地合并聚类对，形成类似树的结构。scipy库提供了这种方法的实现。
import scipy.cluster.hierarchy as shc
# hierarchical clustering using the ward method
clustering = shc.linkage(features_norm_df, method='ward')
# plotting the dendrogram
dend = shc.dendrogram(clustering,
labels=categories.values,
orientation='right',
leaf_font_size=7)
层次聚类模型的结果最好通过树状图进行可视化：
使用树状图可视化层次聚类的结果。图片由作者提供。
我们可以使用树状图来理解聚类的特征。例如，我们可以看到大多数罐装食品被分组（橙色）。橙子也与煎饼/蛋糕混合物聚类。这两者经常在早餐中一起出现。
关键要点
-
发达国家有大量的食品资源盈余。这些盈余通常会被浪费。
-
减少食物浪费可以通过减少温室气体排放对气候变化产生强烈影响；
-
我们可以通过改进食品需求预测模型来实现这一点；
-
对需求时间序列进行聚类可以改善涉及多个时间序列的预测模型。一种方法是为每个聚类训练一个模型。
-
聚类也可以帮助理解数据集中不同群体的特征。
感谢阅读，下次故事见！
参考文献
[1] Jenny Gustavsson, Christel Cederberg, Ulf Sonesson, Robert Van Otterdijk, 和 Alexandre Meybeck. 2011. 全球食品损失和食品浪费. 联合国粮食及农业组织，罗马。
[2] 每周零售食品销售 由经济研究服务提供（许可证：公有领域）
[3] Henderson, Trent, 和 Ben D. Fulcher. “使用 theft 包的基于特征的时间序列分析。” arXiv 预印本 arXiv:2208.06146 (2022).
[4] Rolnick, David, 等. “用机器学习应对气候变化。” ACM Computing Surveys (CSUR) 55.2 (2022): 1–96.
[5] Kasun Bandara, Christoph Bergmeir, 和 Slawek Smyl. 利用递归神经网络对相似系列的时间序列数据库进行预测：一种聚类方法. 专家系统应用, 140:112896, 2020.
气候变化的时间序列：太阳辐射预测
原文：
towardsdatascience.com/time-series-for-climate-change-solar-irradiance-forecasting-a972dac7418f
如何利用时间序列分析和预测应对气候变化
·发布于 Towards Data Science ·8 分钟阅读·2023 年 4 月 5 日
–
图片由 Andrey Grinkevich 提供，来源于 Unsplash
这是系列文章 气候变化的时间序列 的第二部分。文章列表：
- 第一部分: 预测风能
太阳能系统
太阳能是一种越来越普及的清洁能源来源。
太阳光通过光伏设备转化为电能。由于这些设备不产生污染物，因此被认为是清洁能源的来源。除了环境效益外，太阳能因其低成本而具有吸引力。初期投资较大，但长期低成本是值得的。
生产的能源量取决于太阳辐射水平。然而，太阳条件可能会迅速变化。例如，云层可能会突然遮住太阳，降低光伏设备的效率。
因此，太阳能系统依赖于预测模型来预测太阳条件。像 风能预测的情况 一样，准确的预测直接影响这些系统的有效性。
超越能源生产
预测太阳辐射除了用于能源外，还有其他应用，例如：
-
农业：农民可以利用预测来优化作物生产。例如，估算何时种植或收获作物，或优化灌溉系统；
-
土木工程：预测太阳辐射对设计和建造建筑物也很有价值。预测可以用于最大化太阳辐射，从而减少取暖/制冷成本。预测还可以用于配置空调系统，这有助于建筑内能源的高效利用。
挑战及下一步
尽管其重要性，太阳条件高度可变且难以预测。这些条件依赖于几个气象因素，而这些信息有时无法获得。
在本文的其余部分，我们将开发一个太阳辐射预测模型。除此之外，你将学会如何：
-
可视化多变量时间序列；
-
转换多变量时间序列以进行监督学习；
-
根据相关性和重要性评分进行特征选择。
教程：预测太阳辐射
本教程基于由美国农业部收集的数据集。你可以在参考文献[1]中查看更多详细信息。本教程的完整代码可在 Github 上获取：
数据是一个多变量时间序列：在每个时刻，观测值由多个变量组成。这些变量包括以下天气和水文变量：
-
太阳辐射（每平方米瓦特）；
-
风向；
-
积雪深度；
-
风速；
-
露点温度；
-
降水量；
-
蒸汽压力；
-
相对湿度；
-
气温。
该系列从 2007 年 10 月 1 日到 2013 年 10 月 1 日，按小时收集，总计 52,608 次观测。
下载数据后，我们可以使用 pandas 读取它：
import re
import pandas as pd
# src module available here: https://github.com/vcerqueira/tsa4climate/tree/main/src
from src.log import LogTransformation
# a sample here: https://github.com/vcerqueira/tsa4climate/tree/main/content/part_2/assets
assets = 'path_to_data_directory'
DATE_TIME_COLS = ['month', 'day', 'calendar_year', 'hour']
# we'll focus on the data collected at particular station called smf1
STATION = 'smf1'
COLUMNS_PER_FILE = \
{'incoming_solar_final.csv': DATE_TIME_COLS + [f'{STATION}_sin_w/m2'],
'wind_dir_raw.csv': DATE_TIME_COLS + [f'{STATION}_wd_deg'],
'snow_depth_final.csv': DATE_TIME_COLS + [f'{STATION}_sd_mm'],
'wind_speed_final.csv': DATE_TIME_COLS + [f'{STATION}_ws_m/s'],
'dewpoint_final.csv': DATE_TIME_COLS + [f'{STATION}_dpt_C'],
'precipitation_final.csv': DATE_TIME_COLS + [f'{STATION}_ppt_mm'],
'vapor_pressure.csv': DATE_TIME_COLS + [f'{STATION}_vp_Pa'],
'relative_humidity_final.csv': DATE_TIME_COLS + [f'{STATION}_rh'],
'air_temp_final.csv': DATE_TIME_COLS + [f'{STATION}_ta_C'],
}
data_series = {}
for file in COLUMNS_PER_FILE:
file_data = pd.read_csv(f'{assets}/{file}')
var_df = file_data[COLUMNS_PER_FILE[file]]
var_df['datetime'] = \
pd.to_datetime([f'{year}/{month}/{day} {hour}:00'
for year, month, day, hour in zip(var_df['calendar_year'],
var_df['month'],
var_df['day'],
var_df['hour'])])
var_df = var_df.drop(DATE_TIME_COLS, axis=1)
var_df = var_df.set_index('datetime')
series = var_df.iloc[:, 0].sort_index()
data_series[file] = series
mv_series = pd.concat(data_series, axis=1)
mv_series.columns = [re.sub('_final.csv|_raw.csv|.csv', '', x) for x in mv_series.columns]
mv_series.columns = [re.sub('_', ' ', x) for x in mv_series.columns]
mv_series.columns = [x.title() for x in mv_series.columns]
mv_series = mv_series.astype(float)
这段代码生成了以下数据集：
多变量时间序列样本
探索性数据分析
对数尺度的多变量时间序列图。为了可视化，系列被重新采样为每日频率。这是通过取每日平均值来完成的。图片由作者提供。
系列图表明有强烈的年度季节性。辐射水平在夏季达到峰值，其他变量也显示出类似的模式。除了季节性波动外，时间序列的水平在时间上是稳定的。
我们还可以单独可视化太阳辐射变量：
每日总太阳辐射。图片由作者提供。
除了明显的季节性，我们还可以发现一些围绕系列水平的下降峰值。这些情况需要及时预测，以便高效利用备用能源系统。
我们还可以分析每对变量之间的相关性：
显示成对相关性的热图。图片由作者提供。
太阳辐射与一些变量相关。例如，气温、相对湿度（负相关）或风速。
我们已经探讨了如何使用单变量时间序列在上一篇文章中建立预测模型。然而，相关性热图表明，将这些变量纳入模型可能会很有价值。
我们该如何做到这一点？
自回归分布滞后建模入门
自回归分布滞后（ARDL）是一种用于多变量时间序列的建模技术。
ARDL 是一种识别多个变量随时间关系的有用方法。它通过将自回归技术扩展到多变量数据来工作。系列中给定变量的未来值是基于其滞后值和其他变量的滞后值来建模的。
在这种情况下，我们希望基于多个因素的滞后预测太阳辐射，如气温或蒸汽压。
将数据转换为 ARDL 格式
应用 ARDL 方法涉及将时间序列转换为表格格式。这是通过应用时间延迟嵌入到每个变量，然后将结果连接成一个矩阵来完成的。可以使用以下函数来实现：
import pandas as pd
def mts_to_tabular(data: pd.DataFrame,
n_lags: int,
horizon: int,
return_Xy: bool = False,
drop_na: bool = True):
"""
Time delay embedding with multivariate time series
Time series for supervised learning
:param data: multivariate time series as pd.DataFrame
:param n_lags: number of past values to used as explanatory variables
:param horizon: how many values to forecast
:param return_Xy: whether to return the lags split from future observations
:return: pd.DataFrame with reconstructed time series
"""
# applying time delay embedding to each variable
data_list = [time_delay_embedding(data[col], n_lags, horizon)
for col in data]
# concatenating the results in a single dataframe
df = pd.concat(data_list, axis=1)
if drop_na:
df = df.dropna()
if not return_Xy:
return df
is_future = df.columns.str.contains('\+')
X = df.iloc[:, ~is_future]
Y = df.iloc[:, is_future]
if Y.shape[1] == 1:
Y = Y.iloc[:, 0]
return X, Y
这个函数应用于数据如下：
from sklearn.model_selection import train_test_split
# target variable
TARGET = 'Solar Irradiance'
# number of lags for each variable
N_LAGS = 24
# forecasting horizon for solar irradiance
HORIZON = 48
# leaving the last 30% of observations for testing
train, test = train_test_split(mv_series, test_size=0.3, shuffle=False)
# transforming the time series into a tabular format
X_train, Y_train_all = mts_to_tabular(train, N_LAGS, HORIZON, return_Xy=True)
X_test, Y_test_all = mts_to_tabular(train, N_LAGS, HORIZON, return_Xy=True)
# subsetting the target variable
target_columns = Y_train_all.columns.str.contains(TARGET)
Y_train = Y_train_all.iloc[:, target_columns]
Y_test = Y_test_all.iloc[:, target_columns]
我们将预测视野设置为 48 小时。提前预测多个步骤对有效整合多种能源来源到电网中非常有价值。
很难事先确定应该包含多少个滞后。因此，每个变量的值设为 24。这导致总共有 216 个基于滞后的特征。
建立预测模型
在构建模型之前，我们基于日期和时间提取了 8 个额外特征。这些包括年份中的天数或小时等数据，这些数据对建模季节性很有用。
我们通过特征选择减少了解释变量的数量。首先，我们应用相关性过滤器，用于移除与任何其他解释变量相关性大于 95%的特征。然后，我们还基于随机森林的重要性分数应用递归特征消除（RFE）。在特征工程之后，我们使用随机森林训练模型。
我们利用 sklearn 的Pipeline和RandomSearchCV来优化不同步骤的参数：
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import FunctionTransformer
from sklearn.feature_selection import RFE
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import RandomizedSearchCV
from sktime.transformations.series.date import DateTimeFeatures
from src.holdout import Holdout
# including datetime information to model seasonality
hourly_feats = DateTimeFeatures(ts_freq='H',
keep_original_columns=True,
feature_scope='efficient')
# building a pipeline
pipeline = Pipeline([
# feature extraction based on datetime
('extraction', hourly_feats),
# removing correlated explanatory variables
('correlation_filter', FunctionTransformer(func=correlation_filter)),
# applying feature selection based on recursive feature elimination
('select', RFE(estimator=RandomForestRegressor(max_depth=5), step=3)),
# building a random forest model for forecasting
('model', RandomForestRegressor())]
)
# parameter grid for optimization
param_grid = {
'extraction': ['passthrough', hourly_feats],
'select__n_features_to_select': np.linspace(start=.1, stop=1, num=10),
'model__n_estimators': [100, 200]
}
# optimizing the pipeline with random search
model = RandomizedSearchCV(estimator=pipeline,
param_distributions=param_grid,
scoring='neg_mean_squared_error',
n_iter=25,
n_jobs=5,
refit=True,
verbose=2,
cv=Holdout(n=X_train.shape[0]),
random_state=123)
# running random search
model.fit(X_train, Y_train)
# checking the selected model
model.best_estimator_
# Pipeline(steps=[('extraction',
# DateTimeFeatures(feature_scope='efficient', ts_freq='H')),
# ('correlation_filter',
# FunctionTransformer(func=<function correlation_filter at 0x28cccfb50>)),
# ('select',
# RFE(estimator=RandomForestRegressor(max_depth=5),
# n_features_to_select=0.9, step=3)),
# ('model', RandomForestRegressor(n_estimators=200))])
评估模型
我们使用随机搜索选择了一个模型，并结合了a validation split。现在，我们可以评估其在测试集上的预测性能。
# getting forecasts for the test set
forecasts = model.predict(X_test)
forecasts = pd.DataFrame(forecasts, columns=Y_test.columns)
选定的模型仅保留了原始 224 个解释变量中的 65 个。以下是前 20 个特征的重要性：
模型中每个特征的重要性。作者提供的图像。
特征中的一天小时和一年中的一天是排名前四的特征。这一结果突出了数据中季节性效应的强度。除了这些，一些变量的初始滞后期对模型也有用。
主要收获
-
太阳能是一种相关的清洁能源，依赖于太阳辐射；
-
预测太阳辐射是有效整合太阳能到电网中的一个重要方面；
-
太阳辐射依赖于许多变量。这些变量难以建模或可能无法轻易获得；
-
当这些变量可用时，它们表示一个多变量时间序列。这种类型的数据可以用 ARDL 技术建模；
-
你可以通过特征选择过程来估计每个变量的滞后期数量。
感谢阅读，下次故事见！
参考文献
[1] 来自美国爱达荷州西南部四个以西部刺柏为主的实验流域的天气、雪和径流数据。 (许可证：美国公共领域)
[2] Rolnick, David, 等. “利用机器学习应对气候变化。” ACM 计算调查 (CSUR) 55.2 (2022): 1–96。
气候变化中的时间序列：使用深度学习进行精准农业
如何使用时间序列分析和预测来应对气候变化
·发布于 数据科学前沿 ·阅读时间 6 分钟·2023 年 5 月 29 日
–
这是 气候变化中的时间序列 系列的第六部分。文章列表：
精准农业
精准农业旨在改进农业管理。以优化生产，同时节省资源并减少环境影响。
有多种技术致力于提高农业的可持续性。例如：
-
智能灌溉系统：使用传感器优化灌溉过程；
-
精准种植系统：优化种植过程，如种子间距；
-
作物产量预测，帮助农民决定每个季节种植什么作物。
露点温度
在灌溉系统的情况下，过少的灌溉会增加植物的压力并降低作物产量。但过多的灌溉会导致过多的湿气，从而滋生害虫。一个最佳的灌溉过程可以节省大量水资源并保护其他资源。
影响灌溉的一个指标是露点温度。露点温度指示空气中的水分含量。这反过来影响灌溉过程。因此，预测露点温度有助于支持水资源规划。
露点温度的相关性延伸到水文学、气候学和农业等领域。例如，如果预测显示出过多湿气的可能性，可以主动采取害虫控制措施。预测还可以用来预见霜冻，这对植物有害。在这种情况下，农民可以采取预防措施来保护作物。
预测露点温度在能源管理方面也很重要。在高露点温度下，人们倾向于使用空调系统来应对高相对湿度水平。预测这些条件可以用于预测能源需求的增加。这有助于提高电网的效率。
实操：使用深度学习进行露点温度的时空预测
在本文的其余部分，我们将预测多个地点的露点温度。你将学习如何使用深度学习构建时空预测模型。
本教程的完整代码可以在 Github 上找到：
时空预测基础
时空数据是多变量时间序列的一个特例。这些数据集涉及在多个地点观察一个变量，例如露点温度。
这种类型的数据包含了时间依赖性和空间依赖性。特定地点收集的数据点与其滞后值以及附近地点的当前和过去滞后值相关。建模这两种依赖性对于每个地点获得更好的预测可能是重要的。
时空预测通常使用诸如 VAR（向量自回归）或 STAR（时空自回归）等技术。我们将使用 VAR 方法结合深度神经网络。
数据集
我们将使用由美国农业部收集的实际数据集。更多细节请参考文献[1]。该数据集包含了露点温度的信息，这一变量在 6 个附近的站点中被记录。
下载数据后，你可以使用以下代码进行读取：
import pandas as pd
DATE_TIME_COLS = ['month', 'day', 'calendar_year', 'hour', 'water_year']
# reading the data set
data = pd.read_csv(filepath)
# parsing the datetime column
data['datetime'] = \
pd.to_datetime([f'{year}/{month}/{day} {hour}:00'
for year, month, day, hour in zip(data['calendar_year'],
data['month'],
data['day'],
data['hour'])])
data = data.drop(DATE_TIME_COLS, axis=1).set_index('datetime')
data.columns = data.columns.str.replace('_dpt_C', '')
数据的样貌如下：
不同地点的时间序列似乎是相关的。
VAR — 准备数据以进行监督学习
我们将使用 VAR 方法来准备数据以训练深度神经网络。VAR 方法旨在捕捉不同变量之间的时间依赖性。在这种情况下，变量代表在 6 个位置收集的露点温度。
我们可以通过使用滑动窗口将每个变量转换为矩阵格式，然后将结果合并来实现这一点。你可以查看a previous article获取更多关于这一过程的细节。
from sklearn.model_selection import train_test_split
from src.tde import transform_mv_series
N_LAGS, HORIZON = 12, 12
# number of stations
N_STATIONS = data.shape[1]
# leaving last 20% of observations for testing
train, test = train_test_split(data, test_size=0.2, shuffle=False)
# computing the average of each series in the training set
mean_by_location = train.mean()
# mean-scaling: dividing each series by its mean value
train_scaled = train / mean_by_location
test_scaled = test / mean_by_location
# transforming the data for supervised learning
X_train, Y_train = transform_mv_series(train_scaled, n_lags=N_LAGS, horizon=HORIZON)
X_test, Y_test = transform_mv_series(test_scaled, n_lags=N_LAGS, horizon=HORIZON)
然后，我们基于 keras 构建一个堆叠 LSTM 模型。LSTM（长短期记忆）是一种特殊的递归神经网络，可以捕捉时间依赖性。
from keras.models import Sequential
from keras.layers import Dense, Dropout, LSTM, RepeatVector, TimeDistributed
model = Sequential()
model.add(LSTM(64, activation='relu', input_shape=(N_LAGS, N_STATIONS)))
model.add(Dropout(.2))
model.add(RepeatVector(HORIZON))
model.add(LSTM(32, activation='relu', return_sequences=True))
model.add(Dropout(.2))
model.add(TimeDistributed(Dense(N_STATIONS)))
model.compile(optimizer='adam', loss='mse')
定义和编译模型后，我们可以如下训练它：
from keras.callbacks import ModelCheckpoint
# model checkpoint for saving the best model during training
model_checkpoint = ModelCheckpoint(
filepath='best_model_weights.h5',
save_weights_only=True,
monitor='val_loss',
mode='min',
save_best_only=True)
# fitting the model
history = model.fit(X_train, Y_train,
epochs=25,
validation_split=0.2,
callbacks=[model_checkpoint])
训练后，我们可以加载由模型检查点回调保存的最佳权重：
# The best model weights are loaded into the model after training
model.load_weights('best_model_weights.h5')
# predictions on the test set
preds = model.predict_on_batch(X_test)
扩展模型
我们可以通过不同方式扩展这个模型。例如，包含其他气象信息作为解释变量。像露点温度这样的气象数据受到各种因素的影响。它们的包含可能是更好预测性能的关键。
测试其他神经网络配置也可能是有益的。我们应用了堆叠 LSTM，但其他方法也显示出有前景的预测性能。例如 N-BEATS、DeepAR 或 ES-RNN。
我们还可以包括空间信息，例如地理坐标。这样，模型可以改善不同位置之间空间依赖性的建模。
关键要点
-
精准农业旨在优化农业过程。这种优化可以减少农业对环境的影响；
-
露点温度是水文学中的一个关键指标。预测这个变量对各种活动都是有益的，例如灌溉；
-
气象变量通常在几个站点中捕获，这导致了一个时空数据集；
-
深度学习可以用于构建一个时空预测模型。我们使用包含来自六个不同位置的露点温度的数据集训练了一个 LSTM 神经网络；
-
这个模型可以通过包括额外的解释变量、空间信息或更改网络架构来扩展。
感谢阅读，下一个故事见！
参考文献
[1] 来自美国爱达荷州西南部四个以西部刺柏为主的实验流域的天气、降雪和径流数据。（许可证：美国公共领域）
[2] Arikan, Bugrayhan Bickici 等。“北达科他州的露点时间序列预测。” 知识基础工程与科学 2.2 (2021): 24–34。
[3] Liakos, Konstantinos G., 等人。“农业中的机器学习：综述。” 传感器 18.8 (2018): 2674.
时间序列预测：深度学习与统计学——谁能赢？
原文：
towardsdatascience.com/time-series-forecasting-deep-learning-vs-statistics-who-wins-c568389d02df
关于最终难题的全面指南
·发表于Towards Data Science ·14 分钟阅读·2023 年 4 月 5 日
–
使用稳定扩散创建[1]
近年来，深度学习在 NLP 领域取得了显著进展。
时间序列，本质上也具有顺序性，提出了一个问题：如果我们将预训练变压器的全部能力应用于时间序列预测会发生什么？
然而，一些论文，如[2]和[3]，对深度学习模型进行了详细审查。这些论文并没有呈现出完整的画面。即使在 NLP 案例中，有些人将 GPT 模型的突破归因于“更多的数据和计算能力”而不是“更好的 ML 研究”。
本文旨在澄清混淆并提供客观观点，使用来自学术界和行业的可靠数据和来源。具体内容包括：
-
深度学习和统计模型的优缺点。
-
何时使用统计模型，何时使用深度学习。
-
如何处理预测案例。
-
如何通过选择适合你情况和数据集的最佳模型来节省时间和金钱。
让我们深入探讨。
我刚刚发布了AI Horizon Forecast**，**这是一个专注于时间序列和创新 AI 研究的新闻简报。订阅这里以拓宽视野！
[## AutoGluon-TimeSeries：创建强大的集成预测 - 完整教程
亚马逊的时间序列预测框架应有尽有。
aihorizonforecast.substack.com](https://aihorizonforecast.substack.com/p/autogluon-timeseries-creating-powerful?source=post_page-----c568389d02df--------------------------------)
Makridakis 等人的论文[4]
我们不能讨论不同预测模型的进展，而不考虑从**Makridakis 竞赛（M 竞赛）**中获得的见解。
Makridakis 竞赛是一系列大规模挑战，展示了时间序列预测的最新进展。
最近，Makridakis 等人发表了一篇新论文，内容如下：
-
总结了前五次 M 竞赛中的预测现状。
-
提供了各种统计、机器学习和深度学习预测模型的广泛基准。
注： 我们将在本文后面讨论论文的局限性。
基准设置
传统上，Makridakis 及其同事会发布总结最后一次 M 竞赛结果的论文。
然而，作者首次在实验中引入了深度学习模型。为什么？
与 NLP 不同，直到 2018–2019 年，首批 DL 预测模型才足够成熟以挑战传统预测模型。实际上，在 2018 年的 M4 竞赛中，ML/DL 模型排名最后。
图 1： 2018 年 Makridakis 等人检查的八种统计和十种 ML 预测方法的预测准确率（sMAPE）。所有 ML 方法都排在最后。
现在，让我们来看看在新论文中使用的 DL/ML 模型：
-
多层感知机（MLP）： 我们熟悉的前馈网络。
-
WaveNet： 一种自回归神经网络，结合了卷积层（2016 年）。
-
Transformer： 原始 Transformer，介绍于 2017 年。
-
DeepAR： 亚马逊首个成功的自回归网络，结合了 LSTMs（2017 年）
注： 那项研究的深度学习模型已不再是 SOTA（最先进的技术）（更多内容请稍后）。另外，MLP 被认为是 ML 模型而不是“深度”模型。
基准的统计模型是ARIMA和ETS（指数平滑）—— 这些是广为人知且经过实践验证的模型。
此外：
-
ML/DL 模型首先通过超参数调优进行了微调。
-
统计模型是以逐系列的方式训练的。相反，DL 模型是全球性的（一个模型训练于数据集的所有时间序列）。因此，它们利用了跨学习。
-
作者使用了集成方法：从深度学习模型创建了一个Ensemble-DL模型，另一个Ensemble-S由统计模型组成。集成方法是预测的中位数。
-
Ensemble-DL由 200 个模型组成，每个类别有 50 个模型：DeepAR、Transformer、WaveNet 和 MLP。
-
研究使用了 M3 数据集：首先，作者测试了 1,045 个时间序列，然后是完整数据集（3,003 个序列）。
-
作者使用 MASE（均值绝对尺度误差）和 SMAPE（均值绝对百分比误差）来测量预测准确性。这些误差度量常用于预测。
接下来，我们提供了从基准测试中获得的结果和结论的总结。
1. 深度学习模型更好
作者总结道，平均而言，DL 模型优于统计模型。结果显示在 图 2 中：
图 2： 所有模型的平均排名及 95% 置信区间，使用 sMAPE 进行排名。
Ensemble-DL 模型显然优于 Ensemble-S。此外，DeepAR 的结果与 Ensemble-S 非常相似。
有趣的是，图 2 显示虽然 Ensemble-DL 优于 Ensemble-S，但只有 DeepAR 击败了单独的统计模型。这是为什么呢？
我们将在文章后面回答这个问题。
2. 但是，深度学习模型是昂贵的
深度学习模型需要大量的训练时间（和资金）。这是预期之中的。结果显示在 图 3 中：
图 3： SMAPE 与计算时间的关系。ln(CT) 为零对应大约 1 分钟的计算时间，而 ln(CT) 为 2、4、6、8 和 10 分别对应大约 7 分钟、1 小时、7 小时、2 天和 15 天的计算时间。
计算差异是显著的。
因此，降低 10% 的预测误差需要额外大约 15 天的计算时间 (Ensemble-DL)。虽然这个数字看起来很庞大，但有一些因素需要考虑：
-
作者没有说明他们使用了什么类型的硬件。
-
他们也没有提到是否使用了任何并行化或训练优化。
-
如果在集成中使用较少的模型，Ensemble-DL 的计算时间可以显著减少。这在 图 4 中显示了：
图 4： Ensemble-DL 模型中模型数量与 SMAPE 的关系。
我之前提到过，Ensemble-DL 模型是 200 个 DL 模型的集成。
图 4 显示，75 个模型可以以只有 200 个模型三分之一的计算成本实现相当的准确性。如果使用更聪明的集成方法，这个数字可以进一步减少。
最后，当前论文并未探讨深度学习模型的迁移学习能力。我们将在后面讨论这一点。
3. 集成就是你需要的一切
集成的威力是不容置疑的 (图 2，图 3)。
Ensemble-DL 和 Ensemble-SL 都是表现最好的模型。其理念是，每个单独的模型在捕捉不同的时间动态方面都表现出色。结合它们的预测可以识别复杂的模式和进行准确的外推。
4. 短期预测与长期预测
作者调查了模型在短期和长期预测能力上的差异。
确实如此。
图 5展示了每个模型在每个预测期的准确性。例如，第 1 列显示了一步预测误差。类似地，第 18 列显示了第 18 步预测的误差。
图 5： 1045 个系列中每个模型的 sMAPE 误差——数值越低越好（点击这里查看完整图像）
这里有三个关键观察点：
-
首先，长期预测的准确性低于短期预测（这并不意外）。
-
在前 4 个预测期内，统计模型胜出。 在此之后，深度学习模型开始变得更好，Ensemble-DL 胜出。
-
具体来说，在第一个预测期内，Ensemble-S的准确率高出 8.1%。然而，在最后一个预测期，Ensemble-DL的准确率高出 8.5%。
如果你考虑这个问题，就会明白：
-
统计模型是自回归的。随着预测期的增加，误差会累积。
-
相对而言，深度学习模型是多输出模型。因此，它们的预测误差分布在整个预测序列中。
-
唯一的深度学习自回归模型是 DeepAR。这就是为什么 DeepAR 在前几个预测期表现非常好，而其他深度学习模型表现相对较差的原因。
5. 深度学习模型在更多数据下是否有所改善？
在之前的实验中，作者只使用了 M3 数据集中的 1,045 个时间序列。
接下来，作者使用完整的数据集（3,003 个系列）重新进行了实验。他们还分析了每个预测期的预测损失。结果如图 6所示：
图 6： 3003 个系列中每个模型的 sMAPE 误差——数值越低越好（点击这里查看完整图像）
现在，Ensemble-DL和Ensemble-S之间的差距缩小了。统计模型在第一个预测期与深度学习模型相当，但之后，Ensemble-DL 超越了它们。
让我们进一步分析Ensemble-DL和Ensemble-S之间的差异：
图 7： Ensemble-DL 相对于 Ensemble-S 的百分比改善（点击这里查看完整图像）
随着预测步数的增加，深度学习模型优于统计集合模型。
6. 趋势和季节性分析
最后，作者调查了统计模型和深度学习模型如何处理重要的时间序列特征，如趋势和季节性。
为此，作者使用了[5]中的方法论。具体来说，他们拟合了一个多元线性回归模型，将 sMAPE 误差与 5 个关键时间序列特征相关联：预测能力（错误的随机性）、趋势、季节性、线性和稳定性（决定数据正态性的最佳 Box-Cox 参数变换）。结果如图 8所示：
图 8： 不同指标的线性回归系数 0。数值越低越好。
我们观察到：
-
深度学习模型在嘈杂、趋势化和非线性数据上表现更佳。
-
统计模型更适用于季节性和低方差的数据，以及线性关系的数据。
这些见解是非常宝贵的。
因此，在选择适合您的用例的模型之前，进行广泛的探索性数据分析（EDA）并了解数据的性质是至关重要的。
研究的局限性
这篇论文无疑是对当前时间序列预测领域状态的最佳研究之一，但它也存在一些局限性。让我们来看看这些局限性：
缺乏 ML 算法：树 / Boosted Trees
Boosted Trees 模型家族在时间序列预测问题中占有重要地位。
最受欢迎的模型有 XGBoost、LightGBM 和 CatBoost。此外，LightGBM 赢得了 M5 竞赛。
这些模型在表格型数据中表现出色。事实上，到今天为止，Boosted Trees 仍然是表格数据的最佳选择。然而，本文使用的 M3 数据集非常简单，因为它主要包含单变量序列。
在未来的研究中，将 Boosted Trees 添加到数据集中，特别是对于更复杂的数据集，将是一个很好的主意。
选择 M3 作为基准数据集
IJF 期刊主编 Rob Hyndman 教授说：“自 2000 年以来，M3 数据集一直用于测试预测方法；新提出的方法必须超越 M3 才能在 IJF 上发表。”
然而，按现代标准，M3 数据集被认为是小型和简单的，因此不具备现代预测应用和实际场景的代表性。
当然，数据集的选择并不会减少研究的价值。然而，进行一次使用更大数据集的未来基准测试可能会提供有价值的见解。
深度学习模型还不是最先进的（SOTA）。
现在，是时候处理眼前的主要问题了。
该研究中的深度学习模型远未达到最先进水平。
该研究将亚马逊的 DeepAR 识别为在理论预测准确性方面最好的深度学习模型。因此，DeepAR 是唯一一个能够在单独层面上超越统计模型的模型。然而，DeepAR 模型现在已经超过 6 年了。
亚马逊随后发布了其改进版的 DeepAR，称为 Deep GPVAR。实际上，Deep GPVAR 也已过时——亚马逊最新的深度预测模型是 MQTransformer，它于 2020 年发布。
此外，其他强大的模型，如时序融合变换器（TFT）和N-BEATS（最近被 N-HITS 超越），在深度学习集成中也没有被使用。
因此，研究中使用的深度学习模型至少落后于当前技术水平的两代。不容置疑，当前一代深度预测模型将产生更好的结果。
预测并不是一切
准确性在预测中至关重要，但并不是唯一重要的因素。其他关键领域包括：
-
不确定性量化
-
预测可解释性
-
零-shot 学习 / 元学习
-
政策转变隔离
说到零-shot 学习，它是 AI 中最有前途的领域之一。
零-shot 学习是模型在没有专门训练的情况下，正确估计未见数据的能力。这种学习方法更好地反映了人类的认知。
所有的深度学习模型，包括 GPT 模型，都是基于这一原则的。
利用这一原则的首批广受好评的预测模型是N-BEATS / N-HITS。这些模型可以在庞大的时间序列数据集上进行训练，并在完全新的数据上进行预测，其准确性与模型专门训练过的数据相当。
零-shot 学习只是元学习的一个特定实例。自此以来，在时间序列上的元学习取得了进一步的进展。以M6 竞赛为例，其目标是检验数据科学预测与计量经济学是否能像传奇投资者（如沃伦·巴菲特）一样击败市场。获胜解决方案是一种新型架构，其中使用了神经网络和元学习。
不幸的是，本研究并未探索深度学习模型在零-shot 学习设置中的竞争优势。
Nixtla 研究
Nixtla，一家在时间序列预测领域有前景的初创公司，最近发布了对 Makridakis 等人论文[4]的基准后续研究。
具体来说，Nixtla 团队添加了 2 个额外的模型：复杂指数平滑和动态优化 Theta。
这些模型的加入缩小了统计模型和深度学习模型之间的差距。此外，Nixtla 团队正确地指出了这两类模型在成本和资源需求上的显著差异。
确实，许多数据科学家被深度学习的过度炒作所误导，缺乏解决预测问题的正确方法。
我们将在下一部分进一步讨论这个问题。但在此之前，我们需要解决深度学习面临的批评。
深度学习受到批评
过去十年中，深度学习的发展是惊人的。至今尚无减缓的迹象。
然而，每一个威胁到现状的革命性突破通常都会遭遇怀疑和批评。以 GPT-4 为例：这一新发展在下一个十年威胁到 20%的美国工作岗位[6]。
深度学习和变换器在自然语言处理领域的主导地位是不容否认的。然而，面试中人们却提出类似这样的问题：
自然语言处理（NLP）的进步是归功于更好的研究，还是仅仅因为数据的增加和计算能力的提升？
在时间序列预测中，情况更为糟糕。要理解这一点的原因，你必须首先了解传统上是如何处理预测问题的。
在机器学习/深度学习广泛应用之前，预测完全是关于为数据集制定正确的转换。这包括使时间序列平稳，去除趋势和季节性，考虑波动性，并使用如 box-cox 变换等技术。所有这些方法都需要手动干预，并且对数学和时间序列有深刻理解。
随着机器学习的出现，时间序列算法变得更加自动化。你可以在几乎没有预处理的情况下直接应用这些算法（尽管额外的预处理和特征工程总是有帮助的）。如今，许多改进工作主要集中在超参数调优上。
因此，使用高级数学和统计学的人难以理解机器学习/深度学习算法能够超越传统统计模型的事实。而有趣的是，研究人员对一些深度学习概念真正有效的原因和机制却一无所知。
最近文献中的时间序列预测
就我所知，现有文献缺乏足够的证据来说明各种预测模型的优缺点。以下两篇论文是最相关的：
变换器模型对于时间序列预测是否有效？
一篇有趣的论文[2]展示了某些预测变换器模型的弱点。该论文举例说明了现代变换器模型中使用的位置信息编码方案如何未能捕捉时间序列的时间动态。这一点确实是正确的——自注意力机制是不变的。然而，论文未提及那些已经有效解决此问题的变换器模型。
例如，谷歌的 Temporal Fusion Transformer（TFT）使用编码器-解码器 LSTM 层来创建时间感知和上下文感知的嵌入。此外，TFT 使用了为时间序列问题适配的新颖注意力机制，以捕捉时间动态并提供可解释性。
类似地，亚马逊的 MQTransformer 使用其新颖的位置编码方案（上下文相关的季节性编码）和注意力机制（反馈感知注意力）。
我们真的需要深度学习模型来进行时间序列预测吗？
这篇论文 [3] 也很有趣，因为它比较了 统计、提升树、机器学习 和 深度学习 类别的各种预测方法。
不幸的是，它未能达到其标题所述，因为在 12 个模型中表现最好的还是谷歌的 TFT，这是一种纯粹的深度学习模型。论文提到：
…… 表 5 中的结果强调了配置了滚动预测的 GBRT 的竞争力，但也显示了显著强大的基于变换器的模型，如 TFT [12]，确实超越了提升回归树的表现。
总的来说，阅读复杂的预测论文和模型时要小心，特别是关于出版来源的部分。国际预测期刊 (IJF) 就是一个专注于预测的信誉良好的期刊的例子。
如何处理预测问题
这并不简单。每个数据集都是独特的，每个项目的目标也各不相同，使得预测变得具有挑战性。
然而，本文提供了一些可能对大多数方法有益的通用建议。
正如你从本文中了解到的，深度学习模型是预测项目中的一种新兴趋势，但它们仍处于早期阶段。尽管它们具有潜力，但也可能存在陷阱。
不建议立即将深度学习模型作为项目的首选。根据 Makridakis 等人和 Nixtla 的研究，最好从统计模型开始。3–4 个统计模型的集成可能比你预期的更强大。此外，尝试使用提升树，特别是如果你有表格数据的话。对于小型数据集（数千条数据），这些方法可能已足够。
深度学习模型可能提供额外的 3–10% 准确度提升。然而，训练这些模型可能耗时且昂贵。对于一些领域，如金融和零售，这额外的准确度提升可能更具价值，并使使用深度学习模型成为合理选择。更准确的产品销售预测或 ETF 的收盘价可能会转化为数千美元的增量收入。
另一方面，像 N-BEATS 和 N-HITS 这样的深度学习模型具有迁移学习能力。
如果构建了足够大的时间序列数据集，并且一个愿意的实体对这 2 个模型进行预训练并共享其参数，我们可以直接使用这些模型并实现顶尖的预测准确性（或先对我们的数据集进行少量微调）。
结束语
时间序列预测是数据科学的一个关键领域。
但与其他领域相比，它也被严重低估。Makridakis 等人[4]的论文为未来提供了一些有价值的见解，但仍有很多工作和研究需要完成。
此外，深度学习模型在预测中的应用仍未被广泛探索。
例如，多模态架构在深度学习中随处可见。这些架构利用多个领域来学习特定任务。例如，CLIP（由DALLE-2和Stable Diffusion使用）结合了自然语言处理和计算机视觉。
基准 M3 数据集仅包含 3,003 个时间序列，每个序列最多 500 个观测值。相比之下，成功的 Deep GPVAR 预测模型包含平均 44K 个参数。相比之下，Facebook 的 LLaMA 语言 Transformer 模型的最小版本拥有 70 亿个参数，并在 1 万亿个令牌上进行训练。
因此，关于原始问题，没有明确的答案说明哪个模型是最好的，因为每个模型都有其优缺点。
相反，本文旨在提供所有必要的信息，帮助您选择最适合您项目或案例的模型。
感谢阅读！
[## AutoGluon-TimeSeries：创建强大的集成预测——完整教程
亚马逊的时间序列预测框架具备了一切。
aihorizonforecast.substack.com
参考文献
[1] 从 Stable Diffusion 创建，文本提示为“一个蓝青色时间序列抽象，闪亮的数字画，概念艺术”
[2] Ailing Zeng 等人 变压器在时间序列预测中是否有效？（2022 年 8 月）
[3] Shereen Elsayed 等人 我们真的需要深度学习模型来进行时间序列预测吗？（2021 年 10 月）
[4] Makridakis 等人 统计学、机器学习和深度学习预测方法：比较与前进方向 （2022 年 8 月）
[5] Kang 等人 可视化预测算法性能使用时间序列实例空间 （《国际预测学杂志》，2017 年）
[6] Eloundou 等人 GPTs 即 GPTs：大型语言模型对劳动市场潜在影响的早期观察（2023 年 3 月）
文章中使用的所有图像均来自[4]。 M3 数据集以及所有 M 数据集“可以在不需进一步许可的情况下自由使用”。
基于深度学习的时间序列预测（LSTM-RNN）在 PyTorch 中的应用
原文：
towardsdatascience.com/time-series-forecasting-with-deep-learning-in-pytorch-lstm-rnn-1ba339885f0c
一个关于使用 PyTorch 进行单变量时间序列深度学习预测的深入教程
·发布于 Towards Data Science ·阅读时间 12 分钟·2023 年 2 月 9 日
–
介绍
信不信由你，人类在不断被动地预测事物——即使是最微小或看似琐碎的事情。过马路时，我们预测汽车会在哪里以安全地过马路，或者我们尝试预测球在我们试图接住它时会在哪里。我们不需要知道汽车的确切速度或影响球的风向才能完成这些任务——这些能力对我们来说或多或少是自然而然的。这些能力通过一系列事件得到调整，通过多年的经验和实践使我们能够应对我们生活中不可预测的现实。而当需要主动预测大规模现象，如天气或一年后的经济表现时，我们往往因为需要考虑的因素过多而失败。
这就是计算的力量所在——弥补我们无法将即使是最看似随机的事件与未来事件联系起来的不足。众所周知，计算机在执行特定任务时非常擅长——我们可以利用这一点来预测未来。
什么是“时间序列”？
时间序列是指在一段时间内发生的任何可量化的指标或事件。尽管这听起来很琐碎，但几乎任何东西都可以被认为是时间序列。比如一个月内每小时的平均心率，或者一年内每日的股票收盘价，或者某城市一年内每周的交通事故数量。记录这些信息在任何统一的时间段内被认为是时间序列。聪明的人会注意到，这些例子中都有一个频率（每天、每周、每小时等）和一个时间长度（一个月、一年、一天等）来描述事件发生的周期。
对于时间序列，指标是在我们观察指标的时间长度内以统一频率记录的。换句话说，每个记录之间的时间应该是相同的。
在本教程中，我们将探讨如何使用时间序列中的过去数据来预测未来可能发生的情况。
目标
算法的目标是能够接受一系列值，并预测序列中的下一个值。最简单的方法是使用自回归模型，然而，这已经被其他作者广泛讨论过，因此我们将专注于一种更深入的学习方法，即使用递归神经网络。我已将实现笔记本链接在此。本教程中使用的数据集曾在 Kaggle 比赛中使用，可以在这里找到。
数据准备
让我们看看一个示例时间序列。下面的图表显示了 2013 年至 2018 年间油价的一些数据。
作者提供的图像
这只是一个在日期轴上绘制单一数值序列的图表。下表显示了该时间序列的前 10 个条目。从日期列来看，很明显我们有每天频率的价格数据。
date dcoilwtico
2013-01-01 NaN
2013-01-02 93.14
2013-01-03 92.97
2013-01-04 93.12
2013-01-07 93.20
2013-01-08 93.21
2013-01-09 93.08
2013-01-10 93.81
2013-01-11 93.60
2013-01-14 94.27
许多机器学习模型在标准化数据上表现更好。标准化数据的方式是将数据转换为每列的均值为 0，标准差为 1。下面的代码提供了一种使用scikit-learn库进行标准化的方法。
from sklearn.preprocessing import StandardScaler
# Fit scalers
scalers = {}
for x in df.columns:
scalers[x] = StandardScaler().fit(df[x].values.reshape(-1, 1))
# Transform data via scalers
norm_df = df.copy()
for i, key in enumerate(scalers.keys()):
norm = scalers[key].transform(norm_df.iloc[:, i].values.reshape(-1, 1))
norm_df.iloc[:, i] = norm
我们还希望确保我们的数据具有统一的频率——在这个例子中，我们在这 5 年中每天都有油价数据，所以这很好地解决了问题。如果你的数据中不是这种情况，Pandas 提供了几种不同的方法来重新采样数据以适应统一的频率。
序列化
一旦实现了这一点，我们将使用时间序列生成固定长度的片段或序列。在记录这些序列时，我们还会记录序列之后发生的值。例如：假设我们有一个序列：[1, 2, 3, 4, 5, 6]。
通过选择序列长度为 3，我们可以生成以下序列及其相关目标：
[Sequence]: Target
[1, 2, 3] → 4
[2, 3, 4] → 5
[3, 4, 5] → 6
另一种看待这个问题的方法是，我们定义了回顾多少步以预测下一个值。我们将这个值称为 训练窗口，预测的值数目称为 预测窗口。在这个例子中，这些分别是 3 和 1。下面的函数详细说明了如何实现这一点。
# Defining a function that creates sequences and targets as shown above
def generate_sequences(df: pd.DataFrame, tw: int, pw: int, target_columns, drop_targets=False):
'''
df: Pandas DataFrame of the univariate time-series
tw: Training Window - Integer defining how many steps to look back
pw: Prediction Window - Integer defining how many steps forward to predict
returns: dictionary of sequences and targets for all sequences
'''
data = dict() # Store results into a dictionary
L = len(df)
for i in range(L-tw):
# Option to drop target from dataframe
if drop_targets:
df.drop(target_columns, axis=1, inplace=True)
# Get current sequence
sequence = df[i:i+tw].values
# Get values right after the current sequence
target = df[i+tw:i+tw+pw][target_columns].values
data[i] = {'sequence': sequence, 'target': target}
return data
PyTorch 要求我们以以下方式将数据存储在 Dataset 类中：
class SequenceDataset(Dataset):
def __init__(self, df):
self.data = df
def __getitem__(self, idx):
sample = self.data[idx]
return torch.Tensor(sample['sequence']), torch.Tensor(sample['target'])
def __len__(self):
return len(self.data)
然后，我们可以使用 PyTorch 的 DataLoader 来迭代数据。使用 DataLoader 的好处是它内部处理批次和洗牌，因此我们不必担心自己实现这些功能。
训练批次在以下代码执行后就准备好了：
# Here we are defining properties for our model
BATCH_SIZE = 16 # Training batch size
split = 0.8 # Train/Test Split ratio
sequences = generate_sequences(norm_df.dcoilwtico.to_frame(), sequence_len, nout, 'dcoilwtico')
dataset = SequenceDataset(sequences)
# Split the data according to our split ratio and load each subset into a
# separate DataLoader object
train_len = int(len(dataset)*split)
lens = [train_len, len(dataset)-train_len]
train_ds, test_ds = random_split(dataset, lens)
trainloader = DataLoader(train_ds, batch_size=BATCH_SIZE, shuffle=True, drop_last=True)
testloader = DataLoader(test_ds, batch_size=BATCH_SIZE, shuffle=True, drop_last=True)
在每次迭代中，DataLoader 将输出 16（批次大小）个序列及其相关目标，我们将这些序列传递给模型。
模型架构
以下类在 PyTorch 中定义了这一架构。我们将使用一个 LSTM 层，后面跟着一些密集层用于模型的回归部分，并且它们之间有丢弃层。模型将为每个训练输入输出一个单一值。
class LSTMForecaster(nn.Module):
def __init__(self, n_features, n_hidden, n_outputs, sequence_len, n_lstm_layers=1, n_deep_layers=10, use_cuda=False, dropout=0.2):
'''
n_features: number of input features (1 for univariate forecasting)
n_hidden: number of neurons in each hidden layer
n_outputs: number of outputs to predict for each training example
n_deep_layers: number of hidden dense layers after the lstm layer
sequence_len: number of steps to look back at for prediction
dropout: float (0 < dropout < 1) dropout ratio between dense layers
'''
super().__init__()
self.n_lstm_layers = n_lstm_layers
self.nhid = n_hidden
self.use_cuda = use_cuda # set option for device selection
# LSTM Layer
self.lstm = nn.LSTM(n_features,
n_hidden,
num_layers=n_lstm_layers,
batch_first=True) # As we have transformed our data in this way
# first dense after lstm
self.fc1 = nn.Linear(n_hidden * sequence_len, n_hidden)
# Dropout layer
self.dropout = nn.Dropout(p=dropout)
# Create fully connected layers (n_hidden x n_deep_layers)
dnn_layers = []
for i in range(n_deep_layers):
# Last layer (n_hidden x n_outputs)
if i == n_deep_layers - 1:
dnn_layers.append(nn.ReLU())
dnn_layers.append(nn.Linear(nhid, n_outputs))
# All other layers (n_hidden x n_hidden) with dropout option
else:
dnn_layers.append(nn.ReLU())
dnn_layers.append(nn.Linear(nhid, nhid))
if dropout:
dnn_layers.append(nn.Dropout(p=dropout))
# compile DNN layers
self.dnn = nn.Sequential(*dnn_layers)
def forward(self, x):
# Initialize hidden state
hidden_state = torch.zeros(self.n_lstm_layers, x.shape[0], self.nhid)
cell_state = torch.zeros(self.n_lstm_layers, x.shape[0], self.nhid)
# move hidden state to device
if self.use_cuda:
hidden_state = hidden_state.to(device)
cell_state = cell_state.to(device)
self.hidden = (hidden_state, cell_state)
# Forward Pass
x, h = self.lstm(x, self.hidden) # LSTM
x = self.dropout(x.contiguous().view(x.shape[0], -1)) # Flatten lstm out
x = self.fc1(x) # First Dense
return self.dnn(x) # Pass forward through fully connected DNN.
这个类是一个即插即用的 Python 类，我构建它是为了能够动态构建任何大小的神经网络（此类型），基于我们选择的参数——因此请随意调整参数 n_hidden 和 n_deep_players，以添加或删除模型中的参数。更多参数意味着更多模型复杂性和更长的训练时间，所以一定要参考你的使用场景，以确定对数据最合适的参数设置。
作为一个任意选择，我们创建一个具有 5 层全连接层的长短期记忆（LSTM）模型，每层 50 个神经元，最终每个训练样本在每个批次中以单一输出值结束。在这里，sequence_len 指的是训练窗口，而 nout 定义了预测的步数；将 sequence_len 设置为 180 和 nout 设置为 1，意味着模型将回顾过去 180 天（半年），以预测明天会发生什么。
nhid = 50 # Number of nodes in the hidden layer
n_dnn_layers = 5 # Number of hidden fully connected layers
nout = 1 # Prediction Window
sequence_len = 180 # Training Window
# Number of features (since this is a univariate timeseries we'll set
# this to 1 -- multivariate analysis is coming in the future)
ninp = 1
# Device selection (CPU | GPU)
USE_CUDA = torch.cuda.is_available()
device = 'cuda' if USE_CUDA else 'cpu'
# Initialize the model
model = LSTMForecaster(ninp, nhid, nout, sequence_len, n_deep_layers=n_dnn_layers, use_cuda=USE_CUDA).to(device)
模型训练
定义了模型后，我们可以选择损失函数和优化器，设置学习率和训练轮数，并开始训练循环。由于这是一个回归问题（即我们尝试预测一个连续值），一个安全的选择是均方误差作为损失函数。这提供了一种稳健的方法来计算实际值与模型预测值之间的误差。计算公式如下：
图片摘自 Google。
优化器对象存储和计算了反向传播所需的所有梯度。
# Set learning rate and number of epochs to train over
lr = 4e-4
n_epochs = 20
# Initialize the loss function and optimizer
criterion = nn.MSELoss().to(device)
optimizer = torch.optim.AdamW(model.parameters(), lr=lr)
这是训练循环。在每次训练迭代中，我们将计算之前创建的训练集和验证集上的损失：
# Lists to store training and validation losses
t_losses, v_losses = [], []
# Loop over epochs
for epoch in range(n_epochs):
train_loss, valid_loss = 0.0, 0.0
# train step
model.train()
# Loop over train dataset
for x, y in trainloader:
optimizer.zero_grad()
# move inputs to device
x = x.to(device)
y = y.squeeze().to(device)
# Forward Pass
preds = model(x).squeeze()
loss = criterion(preds, y) # compute batch loss
train_loss += loss.item()
loss.backward()
optimizer.step()
epoch_loss = train_loss / len(trainloader)
t_losses.append(epoch_loss)
# validation step
model.eval()
# Loop over validation dataset
for x, y in testloader:
with torch.no_grad():
x, y = x.to(device), y.squeeze().to(device)
preds = model(x).squeeze()
error = criterion(preds, y)
valid_loss += error.item()
valid_loss = valid_loss / len(testloader)
v_losses.append(valid_loss)
print(f'{epoch} - train: {epoch_loss}, valid: {valid_loss}')
plot_losses(t_losses, v_losses)
训练循环的示例输出，显示每个时代的训练和验证损失。
现在模型已经训练完成，我们可以评估我们的预测。
推断
在这里，我们将简单地调用我们训练好的模型来预测未打乱的数据，并查看预测与真实观察值的差异。
def make_predictions_from_dataloader(model, unshuffled_dataloader):
model.eval()
predictions, actuals = [], []
for x, y in unshuffled_dataloader:
with torch.no_grad():
p = model(x)
predictions.append(p)
actuals.append(y.squeeze())
predictions = torch.cat(predictions).numpy()
actuals = torch.cat(actuals).numpy()
return predictions.squeeze(), actuals
历史上的归一化预测与实际油价。图片由作者提供。
初次尝试，我们的预测看起来还不错！而且我们验证损失与训练损失一样低，说明我们没有过拟合模型，因此模型可以被认为具有良好的泛化能力——这对任何预测系统都很重要。
有了对该时间段油价的相对合理估计，让我们看看能否用它来预测未来的情况。
预测
如果我们将历史定义为预测时刻之前的系列，算法则很简单：
-
从历史记录中获取最新有效序列（训练窗口长度）。
-
将最新的序列输入模型，并预测下一个值。
-
将预测值附加到历史记录中。
-
从第 1 步重复进行任意次数的迭代。
一个警告是，根据训练模型时选择的参数，预测得越远，模型越容易受到自身偏差的影响，开始预测均值。因此，如果不必要，我们不希望总是预测过远，因为这会降低预测的准确性。
这在下面的函数中实现：
def one_step_forecast(model, history):
'''
model: PyTorch model object
history: a sequence of values representing the latest values of the time
series, requirement -> len(history.shape) == 2
outputs a single value which is the prediction of the next value in the
sequence.
'''
model.cpu()
model.eval()
with torch.no_grad():
pre = torch.Tensor(history).unsqueeze(0)
pred = self.model(pre)
return pred.detach().numpy().reshape(-1)
def n_step_forecast(data: pd.DataFrame, target: str, tw: int, n: int, forecast_from: int=None, plot=False):
'''
n: integer defining how many steps to forecast
forecast_from: integer defining which index to forecast from. None if
you want to forecast from the end.
plot: True if you want to output a plot of the forecast, False if not.
'''
history = data[target].copy().to_frame()
# Create initial sequence input based on where in the series to forecast
# from.
if forecast_from:
pre = list(history[forecast_from - tw : forecast_from][target].values)
else:
pre = list(history[self.target])[-tw:]
# Call one_step_forecast n times and append prediction to history
for i, step in enumerate(range(n)):
pre_ = np.array(pre[-tw:]).reshape(-1, 1)
forecast = self.one_step_forecast(pre_).squeeze()
pre.append(forecast)
# The rest of this is just to add the forecast to the correct time of
# the history series
res = history.copy()
ls = [np.nan for i in range(len(history))]
# Note: I have not handled the edge case where the start index + n is
# before the end of the dataset and crosses past it.
if forecast_from:
ls[forecast_from : forecast_from + n] = list(np.array(pre[-n:]))
res['forecast'] = ls
res.columns = ['actual', 'forecast']
else:
fc = ls + list(np.array(pre[-n:]))
ls = ls + [np.nan for i in range(len(pre[-n:]))]
ls[:len(history)] = history[self.target].values
res = pd.DataFrame([ls, fc], index=['actual', 'forecast']).T
return res
让我们尝试几个案例。
从系列中间的不同位置进行预测，以便将预测与实际发生的情况进行比较。由于我们编码的预测器可以从任何地方进行预测，并且可以预测任意合理的步数。红线表示预测。请注意，图中的 y 轴显示的是归一化价格。
从 2013 年第三季度预测 200 天。图片由作者提供。
从 2014/15 年年末预测 200 天。图片由作者提供。
从 2016 年第一季度预测 200 天。图片由作者提供。
从数据的最后一天预测 200 天。图片由作者提供。
这只是我们尝试的第一个模型配置！更多地实验架构和实现将使您的模型训练得更好，并更准确地进行预测。
结论
就这样！一个可以预测单变量时间序列中接下来会发生什么的模型。考虑到这种模型可以应用的各种方式和场景，真的很酷。是的，这篇文章仅处理了单变量时间序列，其中只有一个值序列。然而，也有方法可以使用多个测量不同事物的序列来进行预测。这被称为多变量时间序列预测，它主要只需要对模型架构进行一些调整，这些我将在未来的文章中进行介绍。
这种预测模型的真正魔力在于模型的 LSTM 层，它如何处理和记忆序列作为神经网络的递归层。有关不同类型神经网络的更多信息，我强烈推荐3blue1brown 的视频。他有一个很棒的系列，详细介绍了这些算法如何在内部工作，非常直观。
感谢阅读，确保查看我的其他文章！
参考资料：
时间序列数据 — www.kaggle.com/competitions/store-sales-time-series-forecasting/data?select=oil.csv
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
