Merlion笔记（四）：添加一个新的预测模型

本文将实现一个预测模型，其预测值正好等于该时间点的前一个观测值。有关更真实的示例，请参阅对 Sarima 的实现。

1 模型配置类

创建新模型的第一步是定义一个配置类，该类继承自 ForecasterConfig:

from merlion.models.forecast.base import ForecasterConfig

class RepeatRecentConfig(ForecasterConfig):
    def __init__(self, max_forecast_steps=None, **kwargs):
        super().__init__(max_forecast_steps=max_forecast_steps, **kwargs)

2 模型类

接下来，定义模型本身，该模型必须继承自 ForecasterBase 基类，并实现所有抽象方法。

from collections import OrderedDict
from typing import List, Tuple

import numpy as np
import pandas as pd

from merlion.models.forecast.base import ForecasterBase
from merlion.utils.time_series import to_pd_datetime


class RepeatRecent(ForecasterBase):
    # RepeatRecent 的配置类是上面定义的 RepeatRecentConfig
    config_class = RepeatRecentConfig

    @property
    def require_even_sampling(self):
        """
        许多预测模型假设输入的时间序列是均匀采样的。这个模型不需要这种假设，因此重写该属性。
        """
        return False

    def __init__(self, config):
        """
        设置模型配置和其他局部变量。在这里，我们将 most_recent_value 初始化为 None。
        """
        super().__init__(config)
        self.most_recent_value = None


    def _train(self, train_data: pd.DataFrame, train_config=None) -> Tuple[pd.DataFrame, None]:
        # 训练模型。在这里，我们只是收集每个单变量的最新观察值。
        # 列表推导式，用来遍历 train_data 的每一列（键值对形式）。对于每一列，k 是列名，v.values[-1] 
        # 是该列的最后一个观测值。最终生成一个列表，其中每个元素是 (列名, 最近观测值) 这样的元组。
        self.most_recent_value = [(k, v.values[-1]) for k, v in train_data.items()]

        # 模型的目标值是每个时间序列的前一个值，即每一行的预测值是上一时间点的实际值。
        # 将一个全 0 的数组与 train_data（去掉最后一行后的数据）拼接起来，形成一个新的数组 pred，这个数组的每一行都是前一个时间点的数值。
        pred = np.concatenate((np.zeros((1, self.dim)), train_data.values[:-1]))
        train_forecast = pd.DataFrame(pred, index=train_data.index, columns=train_data.columns)

        # 这个模型没有误差的概念
        train_stderr = None

        # 返回训练的预测结果和标准误差
        return train_forecast, train_stderr

    def _forecast(self, time_stamps: List[int], time_series_prev: pd.DataFrame = None, return_prev=False) -> Tuple[pd.DataFrame, None]:

        # 如果提供了 time_series_prev，则使用其最近的值。否则，使用从训练数据中存储的最近值
        if time_series_prev is not None:
            most_recent_value = [(k, v.values[-1]) for k, v in time_series_prev.items()]
        else:
            most_recent_value = self.most_recent_value

        # 预测值只是将最近的一个值重复用于每一个未来的时间点。
        i = self.target_seq_index # 目标序列的索引
        datetimes = to_pd_datetime(time_stamps) # 测试序列的时间戳
        name, val = most_recent_value[i]
        forecast = pd.DataFrame([val] * len(datetimes), index=datetimes, columns=[name])

        # 如果需要，给 time_series_prev 的 target_seq_index 预加上“预测”值。
        if return_prev and time_series_prev is not None:
            pred = np.concatenate(([0], time_series_prev.values[:-1, i]))
            prev_forecast = pd.DataFrame(pred, index=time_series_prev.index, columns=[name])
            forecast = pd.concat((prev_forecast, forecast))

        return forecast, None

3 运行模型：一个简单的例子

尝试在一些实际数据上运行这个模型！我们将首先从 M4 数据集中获取时间序列并将其可视化。

import matplotlib.pyplot as plt
import pandas as pd

from merlion.utils import TimeSeries, UnivariateTimeSeries
from ts_datasets.forecast import M4

time_series, metadata = M4(subset="Hourly")[0]

# Visualize the full time series
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
ax.plot(time_series)

# Label the train/test split with a dashed line
ax.axvline(time_series[metadata["trainval"]].index[-1], ls="--", lw=2, c="k")

plt.show()

现在，将数据分成训练和测试部分，并在其上运行我们的预测模型。

train_data = TimeSeries.from_pd(time_series[metadata["trainval"]])
test_data = TimeSeries.from_pd(time_series[~metadata["trainval"]])

# Initialize a model & train it. The dataframe returned & printed
# below is the model's "forecast" on the training data. None is
# the uncertainty estimate.
model = RepeatRecent(RepeatRecentConfig())
model.train(train_data=train_data)

(                        H1
 time                      
 2015-01-07 12:00:00    0.0
 2015-01-07 13:00:00  605.0
 2015-01-07 14:00:00  586.0
 2015-01-07 15:00:00  586.0
 2015-01-07 16:00:00  559.0
 ...                    ...
 2015-02-05 11:00:00  820.0
 2015-02-05 12:00:00  790.0
 2015-02-05 13:00:00  784.0
 2015-02-05 14:00:00  752.0
 2015-02-05 15:00:00  739.0
 
 [700 rows x 1 columns],
 None)

# Let's run our model on the test data now
forecast, err = model.forecast(test_data.to_pd().index)
print("Forecast")
print(forecast)
print()
print("Error")
print(err)

Forecast
                        H1
time                      
2015-02-05 16:00:00  684.0
2015-02-05 17:00:00  684.0
2015-02-05 18:00:00  684.0
2015-02-05 19:00:00  684.0
2015-02-05 20:00:00  684.0
2015-02-05 21:00:00  684.0
2015-02-05 22:00:00  684.0
2015-02-05 23:00:00  684.0
2015-02-06 00:00:00  684.0
2015-02-06 01:00:00  684.0
2015-02-06 02:00:00  684.0
2015-02-06 03:00:00  684.0
2015-02-06 04:00:00  684.0
2015-02-06 05:00:00  684.0
2015-02-06 06:00:00  684.0
2015-02-06 07:00:00  684.0
2015-02-06 08:00:00  684.0
2015-02-06 09:00:00  684.0
2015-02-06 10:00:00  684.0
2015-02-06 11:00:00  684.0
2015-02-06 12:00:00  684.0
2015-02-06 13:00:00  684.0
2015-02-06 14:00:00  684.0
2015-02-06 15:00:00  684.0
2015-02-06 16:00:00  684.0
2015-02-06 17:00:00  684.0
2015-02-06 18:00:00  684.0
2015-02-06 19:00:00  684.0
2015-02-06 20:00:00  684.0
2015-02-06 21:00:00  684.0
2015-02-06 22:00:00  684.0
2015-02-06 23:00:00  684.0
2015-02-07 00:00:00  684.0
2015-02-07 01:00:00  684.0
2015-02-07 02:00:00  684.0
2015-02-07 03:00:00  684.0
2015-02-07 04:00:00  684.0
2015-02-07 05:00:00  684.0
2015-02-07 06:00:00  684.0
2015-02-07 07:00:00  684.0
2015-02-07 08:00:00  684.0
2015-02-07 09:00:00  684.0
2015-02-07 10:00:00  684.0
2015-02-07 11:00:00  684.0
2015-02-07 12:00:00  684.0
2015-02-07 13:00:00  684.0
2015-02-07 14:00:00  684.0
2015-02-07 15:00:00  684.0

Error
None

4 可视化

# Qualitatively, we can see what the forecaster is doing by plotting
print("Forecast w/ ground truth time series")
fig, ax = model.plot_forecast(time_series=test_data,
                              time_series_prev=train_data,
                              plot_time_series_prev=True)
plt.show()

print()
print("Forecast without ground truth time series")
fig, ax = model.plot_forecast(time_stamps=test_data.to_pd().index,
                              time_series_prev=train_data,
                              plot_time_series_prev=True)

Forecast w/ ground truth time series

​
Forecast without ground truth time series

5 定量评估

也可以对模型进行定量评估。计算模型预测结果与真实数据的对称平均百分比误差（sMAPE，symmetric Mean Average Percent Error）。

from merlion.evaluate.forecast import ForecastMetric
smape = ForecastMetric.sMAPE.value(ground_truth=test_data, predict=forecast)
print(f"sMAPE = {smape:.3f}")

sMAPE = 20.166

6 定义一个基于预测器的异常检测器

将一个预测模型转换为异常检测模型是非常简单的。只需要在合适的目录下创建一个新文件，并定义包含一些基本头部的类结构。通过多重继承 ForecastingDetectorBase 类，大部分繁重的工作都可以自动处理。

任何基于预测的异常检测器返回的异常评分，都是基于预测值与真实时间序列值之间的残差。

from merlion.evaluate.anomaly import TSADMetric
from merlion.models.anomaly.forecast_based.base import ForecastingDetectorBase
from merlion.models.anomaly.base import DetectorConfig
from merlion.post_process.threshold import AggregateAlarms
from merlion.transform.normalize import MeanVarNormalize


# 定义一个配置类，该类按顺序继承自 RepeatRecentConfig 和 DetectorConfig
class RepeatRecentDetectorConfig(RepeatRecentConfig, DetectorConfig):
    # 设置一个默认的异常评分后处理规则
    _default_post_rule = AggregateAlarms(alm_threshold=3.0)

    # 默认的数据预处理变换是均值-方差归一化，
    # 这样异常评分大致与 z-score 对齐
    _default_transform = MeanVarNormalize()

# 定义一个模型类，该类按顺序继承自 ForecastingDetectorBase 和 RepeatRecent
class RepeatRecentDetector(ForecastingDetectorBase, RepeatRecent):
    # 我们只需要设置配置类
    config_class = RepeatRecentDetectorConfig

# Train the anomaly detection variant
model2 = RepeatRecentDetector(RepeatRecentDetectorConfig())
model2.train(train_data)

                     anom_score
time                           
2015-01-07 12:00:00   -0.212986
2015-01-07 13:00:00   -0.120839
2015-01-07 14:00:00    0.000000
2015-01-07 15:00:00   -0.171719
2015-01-07 16:00:00   -0.305278
...                         ...
2015-02-05 11:00:00   -0.190799
2015-02-05 12:00:00   -0.038160
2015-02-05 13:00:00   -0.203519
2015-02-05 14:00:00   -0.082679
2015-02-05 15:00:00   -0.349798

[700 rows x 1 columns]

# Obtain the anomaly detection variant's predictions on the test data
model2.get_anomaly_score(test_data)

                     anom_score
time                           
2015-02-05 16:00:00   -0.413397
2015-02-05 17:00:00   -0.756835
2015-02-05 18:00:00   -0.966714
2015-02-05 19:00:00   -1.202032
2015-02-05 20:00:00   -1.291072
2015-02-05 21:00:00   -1.380111
2015-02-05 22:00:00   -1.341952
2015-02-05 23:00:00   -1.246552
2015-02-06 00:00:00   -1.163873
2015-02-06 01:00:00   -0.953994
2015-02-06 02:00:00   -0.686876
2015-02-06 03:00:00   -0.286198
2015-02-06 04:00:00    0.178079
2015-02-06 05:00:00    0.559676
2015-02-06 06:00:00    0.928554
2015-02-06 07:00:00    1.246552
2015-02-06 08:00:00    1.329232
2015-02-06 09:00:00    1.348311
2015-02-06 10:00:00    1.316512
2015-02-06 11:00:00    1.081193
2015-02-06 12:00:00    0.756835
2015-02-06 13:00:00    0.540597
2015-02-06 14:00:00    0.426117
2015-02-06 15:00:00    0.108119
2015-02-06 16:00:00   -0.311638
2015-02-06 17:00:00   -0.712316
2015-02-06 18:00:00   -0.966714
2015-02-06 19:00:00   -1.214752
2015-02-06 20:00:00   -1.316512
2015-02-06 21:00:00   -1.373751
2015-02-06 22:00:00   -1.399191
2015-02-06 23:00:00   -1.316512
2015-02-07 00:00:00   -1.221112
2015-02-07 01:00:00   -1.049393
2015-02-07 02:00:00   -0.737755
2015-02-07 03:00:00   -0.381598
2015-02-07 04:00:00    0.076320
2015-02-07 05:00:00    0.489717
2015-02-07 06:00:00    0.814075
2015-02-07 07:00:00    0.966714
2015-02-07 08:00:00    0.979434
2015-02-07 09:00:00    0.922194
2015-02-07 10:00:00    0.782275
2015-02-07 11:00:00    0.642356
2015-02-07 12:00:00    0.457917
2015-02-07 13:00:00    0.222599
2015-02-07 14:00:00    0.120839
2015-02-07 15:00:00   -0.158999

# Visualize the anomaly detection variant's performance, with filtered anomaly scores
fig, ax = model2.plot_anomaly(test_data, time_series_prev=train_data,
                              filter_scores=True, plot_time_series_prev=False,
                              plot_forecast=True)

​
​