处理时间序列中的间隔-开发者社区

原文：towardsdatascience.com/handling-gaps-in-time-series-dc47ae883990

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/25bd6f39a6c87567493d47435a112a45.png

图片由 Willian Justen de Vasconcellos 在 Unsplash 提供

时间是物理学和自然界中最定义明确的连续体。因此，在时间序列数据集中——一系列按时间顺序排列的观察值——连续性的重要性也就不足为奇了。

单独这个概念就推动了本文背后的动机。现实世界的数据集可能由于各种原因存在缺失值，例如传感器故障、数据摄取失败，或者在某些时间段内信息的缺失。然而，这并不改变您特征数据生成过程的本质。

因此，理解这些中断的原因，并在时间序列数据集中分析和处理它们，对于任何后续任务都是至关重要的。

本文的目标

在对您的时序数据进行全面探索性分析之后，您可能会发现存在相当数量的缺失值。通过理解您数据的本质，您应该能够区分代表缺失的间隔和表示实际中断的间隔，将其视为间歇序列。

本文将重点关注第一种场景——缺失值分析和评估插补结果的方法。虽然执行插补的实际技术有很多 [1][2]，但我将详细阐述评估这些结果质量的方法，以便您能够评估在您的数据集中何时何地应用它们。

在我们开始之前，值得提一下，插补可能并不总是最佳选择。记住，在时间序列问题中，你的目标通常是构建数据生成过程的数学表示（是什么导致数据发生），但这个概念，根据定义，是一个我们无法完全达到的柏拉图理想[3]。缺失值所做的是模糊这个理想，但不良的插补可能会让您离它更远。

考虑到这一点，我们将处理一个假设场景，我们的目标是预测空气中细颗粒物（也称为 PM 2.5）的水平，它是空气污染和空气质量指数的主要贡献者之一。那么，让我们看看我们的数据集。

数据集描述

数据包含从 2016 年 1 月 1 日到 2022 年 7 月 3 日，在加拿大不列颠哥伦比亚省温哥华市四个监测站，细颗粒物 PM 2.5（直径为 2.5 微米及以下的细颗粒物）的平均一小时测量值，单位为µg/m3（每立方米微克）。

该数据集是原始数据集的子集，由不列颠哥伦比亚省数据目录 [9] 发布，每个站点包含 57,014 条 PM 2.5 记录，有 4,958 个缺失值。它可在开放政府许可 – 不列颠哥伦比亚省 [10] 下获取。

importpandasaspd master_df=pd.read_csv("data/2016_2021_master_df.csv")master_df=pd.melt(master_df,value_vars=["Vancouver_Clark_Drive_PM25","Vancouver_International_Airport_#2_PM25","North_Vancouver_Mahon_Park_PM25","North_Vancouver_Second_Narrows_PM25"],ignore_index=False)display(master_df.head())

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/4be5d65043a4328f15c0868496828e53.png

主数据帧

图书馆和依赖项

本文的代码示例可以在任何 Python 3.8+环境中重现。

我们将使用Pandas和Numpy进行数据处理，Matplotlib和Seaborn进行可视化，以及SciPy进行一些度量。我们还将使用Darts，这是一个具有用户友好实现的众多预测和异常检测模型的库。

我将利用我在之前的帖子中介绍的一些自定义可视化函数。您将在我叫用名为tshelpers的模块时看到它们。

数据预处理

在预处理过程中，我们希望识别缺失值并量化缺失和连续性。为此，我们将首先引入一个标志isMissing。

# Dummy variable to track missing samples and missing sample lengthmaster_df["isMissing"]=np.where(master_df["PM 2.5"].isnull(),1,0)master_df.head()

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/505ba5a199043527dff9251bcac848b1.png

isMissing 标志

下面的代码片段将四个站点分离成四个独立的数据帧。

# Isolating stations on independent data framesdatasets={}forstationinstations:datasets[station]=master_df[master_df["Station"]==station]datasets[station]=datasets[station][["PM 2.5","isMissing"]]datasets[station].reset_index(inplace=True)

Darts 有一系列有用的方法和函数来帮助量化定位缺口和完整序列，例如[TimeSeries.gap](https://unit8co.github.io/darts/generated_api/darts.timeseries.html#darts.timeseries.TimeSeries.gaps)和[darts.utils.missing_values.extract_subseries](https://unit8co.github.io/darts/generated_api/darts.utils.missing_values.html?highlight=extract+subseries#darts.utils.missing_values.extract_subseries)。我强烈建议你在项目中已经使用这个库的情况下查看这些内容，但由于 Darts 是一个相当重的依赖项，我将在下面定义一个简单的实现，你可以复制并满足你分析缺失性的需求。

以下函数计算缺失和非缺失区间的长度，以及它们的起始和结束时间戳。

defcount_sequences(data,time,value):# Indexer for nonmissing value on datais_not_nan=~data[value].isna()# Auxiliary indexer to group sequences# diff and cumsum aggregate nonmissing sequencesgroup_idx=is_not_nan.diff().cumsum().fillna(0)# fillna(0) resolves cumsum's 'NaN' on idx = 0# Not-missing counternot_nan_counts=is_not_nan.groupby(group_idx).sum()# Instantiate sequence lengths DataFrame and retrieve position indices and timesequences_df=pd.Series(np.arange(len(data))).groupby(group_idx).agg(['min','max'])sequences_df['seq_start_time']=sequences_df['min'].map(data[time])sequences_df['seq_end_time']=sequences_df['max'].map(data[time])sequences_df['not_nan_count']=not_nan_counts sequences_df['nan_count']=(sequences_df['max']-sequences_df['min'])-(sequences_df['not_nan_count']-1)# Assert sum of sequence lengths == total series lengthassertsum(sequences_df[['not_nan_count','nan_count']].sum())==data.shape[0]# Tidy upsequences_df.rename(columns={'min':'seq_start_idx','max':'seq_end_idx'},inplace=True)sequences_df.sort_values('seq_start_idx',inplace=True)sequences_df.reset_index(drop=True,inplace=True)returnsequences_df

孤立与连续缺失值

从论文“在条件监测数据集中缺失值的三个步骤插补”中很好地阐述的概念，我们可以根据缺失段的长度和相邻变量的完整性来分类缺失区间[5]。

我们将借鉴他们的原始定义，并将其推广到单变量序列，其中我们将缺失段分类为连续或孤立。因此，为了本文的目的，我们定义如下：

孤立缺失值 (IMV)：一个单独的观测值或长度等于或小于 T 的观测序列，其中某个变量的数据缺失。
连续缺失值 (CMV)：一个长度大于 T 的观测序列，其中某个变量的数据缺失。

在这里，T 是用于将缺失序列分类为孤立（IMV）或连续（CMV）的阈值。这个定义可以推广到单变量或多变量场景，并决定了是否将完整段（缺失样本之间的部分）视为独立的。

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/88a0091e5b1b0aba004bde3e7359378a.png

孤立（IMV）与连续缺失值（CMV），阈值 T = 2。图片由作者提供。

设置 IMV 阈值 T 为 2 小时后，我们现在可以使用数据集中的count_sequences()函数来测量缺失段的长度，并将它们分类为单个或连续缺失。

IMV_THRESHOLD=2# T (in hours)metadata={}forstationinstations:metadata[station]=count_sequences(datasets[station],time='DATE_PST',value='PM 2.5')# Continuous Missing Values (CMV) flagmetadata[station]['isCMV']=(metadata[station]['nan_count']>IMV_THRESHOLD).astype(int)# Isolated Missing Values (IMV) flagmetadata[station]['isIMV']=(metadata[station]['nan_count'].isin([iforiinrange(1,IMV_THRESHOLD+1)])).astype(int)metadata['North_Vancouver_Second_Narrows_PM25'].head()

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/a7f58f2fbf88eec05184ee1a53f8579f.png

缺失值的元数据

使用此元数据，我们现在可以在原始数据框中创建一个新的变量seqMissing，以存储缺失序列的长度，并重新组合我们的主数据框。

# Missing values metadataforstationinstations:fori,rowinmetadata[station].iterrows():seq_range=range(row['seq_start_idx'],row['seq_end_idx']+1)datasets[station].loc[seq_range,'seqMissing']=row['nan_count']datasets[station]['seqMissing']=datasets[station]['seqMissing'].astype(int)# Recomposing Master DataFramemaster_df=pd.concat(datasets).reset_index(level=0)# Renaming station columnmaster_df["Station"]=master_df["level_0"]master_df=master_df[["DATE_PST","Station","PM 2.5","isMissing","seqMissing"]]# Redefining DATE_PST indexmaster_df.set_index("DATE_PST",inplace=True)master_df[master_df["isMissing"]==1].head()

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/2e1209cccd6811868f70960677868d59.png

附加到主数据框的缺失元数据

在这个阶段，你可能想知道如何定义 IMV 阈值 T。

你可以从查看数据集的粒度开始，当它用于下游数据科学或分析任务时。如果你知道你将把数据汇总到六小时周期，你可以定义T=6并将等于或低于该值的任何内容视为单个或短序列。我们将在稍后讨论这样做的原因。

可视化缺失值

对于我们的第一次可视化，我们将查看数据集中单个与连续缺失值的分布。我们希望根据每个变量的总缺失值来量化它们。

以下代码将汇总（求和）IMV 和 CMV 间隙的总长度，并在饼图中绘制结果（是的，饼图——如果你有超过 5 个变量，请尝试使用条形图）。

# Continuous Missing Values &amp; Isolated Missing Values counttot_missing=[]stations_cmv_imv=[]stations_cmv_imv_label=[]# Count of CMV and IMV intervals and percentage out of total missingforstationinstations:# Total missing valuestot_missing_current=datasets[station]['isMissing'].sum()tot_missing.append(tot_missing_current)# Total Continuous Missing Valuetot_cmv=metadata[station].loc[metadata[station]['isCMV']==1,'nan_count'].sum()stations_cmv_imv.append(tot_cmv)stations_cmv_imv_label.append(f"{tot_cmv/tot_missing_current*100:.0f}%nCMV")# Total Isolated Missing Valuetot_imv=tot_missing_current-tot_cmv stations_cmv_imv.append(tot_imv)stations_cmv_imv_label.append(f"{tot_imv/tot_missing_current*100:.0f}%nIMV")# Stacked pie charts with IMV and CMS percentagesfig,ax=plt.subplots(figsize=(9,9))ax.axis("equal")width=0.3# Outer pie chart (total missing values per station)pie,texts1=ax.pie(tot_missing,radius=1,labels=tot_missing,labeldistance=0.85,textprops={"fontsize":14,"weight":"bold"})fortintexts1:t.set_horizontalalignment("center")plt.legend(pie,stations,loc=(0,-0.05))plt.setp(pie,width=width,edgecolor="white")# Inner pie chart (percentage of IMV and CMV out of total missing)pie2,texts2=ax.pie(stations_cmv_imv,radius=1-width,labels=stations_cmv_imv_label,labeldistance=0.78,textprops={"weight":"bold"})fortintexts2:t.set_horizontalalignment("center")plt.setp(pie2,width=width,edgecolor="white")plt.title("Missing Values per Station",fontsize=14,weight="bold",y=0.96)plt.suptitle("Isolated Missing Values (IMV) and Continuous Missing Values (CMV)",fontsize=12,y=0.85,)plt.show()

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/42dc195ee45fcb8dd0ad6bbb795c41af.png

是的，一个饼图

我们可以看到，我们大部分的缺失值都是 3 小时或更长的序列的一部分。图表还显示，北温哥华第二海峡的缺失值数量最多，并且绝大多数来自长序列。

以下可视化将展示 CMV 序列长度的分布情况。我们将绘制序列长度计数的直方图和补充箱线图。

# Distribution of CMV Sequence Lengthforstationinstations:_,(ax_hist,ax_box)=plt.subplots(2,sharex=True,gridspec_kw={"height_ratios":(.9,.10)},figsize=(12,4))plot_title=station[:-5].replace('_',' ')cmv_counts=metadata[station].loc[metadata[station]['isCMV']==1,'nan_count'].to_list()sns.histplot(x=cmv_counts,kde=True,bins=40,ax=ax_hist)sns.boxplot(x=cmv_counts,width=0.7,color='0.5',linecolor='#003366',flierprops={'marker':'.'},ax=ax_box)ax_hist.set(title=plot_title,ylabel=f'Count of Sequences')ax_hist.yaxis.set_major_locator(MaxNLocator(integer=True))plt.suptitle('Continuous Missing Value (CMV) | Sequence Length Distribution',fontweight='bold',fontsize=12)ax_box.set(yticks=[],xlabel='Missing Values Sequence Length')sns.despine(ax=ax_hist)sns.despine(ax=ax_box,left=True)

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/b1d9f8e6be2b117201dbaca1db7b3e67.png

CMV 分布：温哥华克拉克大道

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/6806702a43094b7be5bb57efa8db2a24.png

CMV 分布：温哥华国际机场

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/ffb195dfc07bc3c1b66cfcca8cdabe00.png

CMV 分布：北温哥华马洪公园

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/213ffa1d6557ebcf0a19fca1df6b4d23.png

CMV 分布：北温哥华第二海峡

我们可以看到，大部分长序列的长度大约在 3 到 24 小时之间，有些异常值超出了这个范围。这也清楚地解释了为什么北温哥华第二海峡和温哥华克拉克大道车站是缺失值最多的地方：它们都有 350 小时或更长的序列。

这些结果暗示，我们可能可以从处理缺失值的三种不同方法中受益。一种用于 IMV 或短序列的插补方法，另一种（或相同的）用于 3 到 24 小时长的序列，可能将超过这个长度的间隙视为中断者，在这种情况下，考虑中间的完整段作为独立的部分可能更好，因为我们可能缺乏一个准确插补它们的好模型。

为实验而选择的子集间隔

为了评估短序列插补的质量，我们将寻找可以人工创建缺失值以进行插补的完整区间，并将结果与真实数据进行比较。

首先，我们想要了解我们的数据集时间跨度内缺失值的位置以及它们的分布情况。我们将使用热力图来完成这项工作。

fromtshelpers.plotimportplot_missing plot_missing(master_df.pivot(columns="Station",values="PM 2.5"))

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/c8afcb5f43ad1c2eee96f28ea792a7cf.png

缺失值热力图

它们分布得相当分散。以下代码片段将遍历每个站点并捕获没有缺失值的月份。这些月份将被分配为我们的实验子集，我们将人为地创建间隔并进行插补以评估插补方法。

# Looking for months without missing valuesmonths_complete={}forstationinstations:months_complete[station]=[]station_subset=master_df[master_df["Station"]==station]foryearinrange(2016,2022):ifyear==2022:range_max=7# Limiting upper range for 2022 as data goes up to Julyelse:range_max=13formonthinrange(1,range_max):ifmonth==12:subset=station_subset[datetime(year,month,1,1):datetime(year,month,31,23)]else:subset=station_subset[datetime(year,month,1,1):datetime(year,month+1,1)]ifsubset['PM 2.5'].isna().sum()==0:months_complete[station].append((month,year))# Experimental subsetsexp_subsets={}forstationinmonths_complete:exp_subsets[station]={}formonthsinmonths_complete[station]:# Iterating on list of tuples (month, year) | (4, 2016)month=months[0]year=months[1]ifmonth==12:next_month=1next_year=year+1else:next_month=month+1next_year=year exp_subsets[station][f'{month}-{year}']=master_df[master_df["Station"]==station][datetime(year,month,1):datetime(next_year,next_month,1)]forstationinexp_subsets:print(f'{station}:n{exp_subsets[station].keys()}n')

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/fdeb7e21a26c0f33efbb3673acb3d167.png

我们总共有 20 个实验子集，每个子集大约有 720 个完整的观测值（24 小时的 30 天）。

人工缺失数据

为了程序化创建人工缺失序列，我们将定义一个函数，该函数接受开始日期和结束日期、所需缺失间隔的长度以及创建它们的位置（相对于序列的开始或结束）的填充。

defcreate_missing(data,time=None,value=None,start=None,end=None,missing_length=1,missing_index="end",padding=24,verbose=True,):# Missing data subsetsubset_missing=subset.copy()subset_array=subset_missing[value].to_numpy()ifmissing_index=="end":missing_start=padding+missing_length missing_end=padding subset_array[-missing_start:-missing_end]=np.NaNelifmissing_index=="start":missing_start=padding missing_end=padding+missing_length subset_array[missing_start:missing_end]=np.NaN subset_missing[value]=subset_arrayreturnsubset,subset_missing

可以通过以下方式生成一个长度为 6 小时的人工缺失间隔，并在序列末尾填充 24 小时。

fromtshelpers.plotimportplot_compare station_subset=master_df[master_df["Station"]=="North_Vancouver_Mahon_Park_PM25"]subset,subset_missing=create_missing(data=station_subset,value="PM 2.5",start=datetime(2019,9,1),end=datetime(2019,9,10),missing_length=6,padding=24,missing_index="end",)plot_compare(data=subset,# Subset to plotdata_missing=subset_missing,# Subset with missing values to plotvalue="PM 2.5",# Variable to plot from subsetvalue_missing="PM 2.5",# Variable to plot from subset_missingmissing_only=True,# Plots only the section with missing values from subset_missingstart=datetime(2019,9,1),end=datetime(2019,9,10),fill=True,data_label="Artificial missing values",# Label for subsetdata_missing_label="Original data",# Label for subset_missing)

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/3ce707f7905698cbb92538325112103c.png

人工缺失值示例

这个概念将指导我们接下来的实验方法。

短序列插补

对于短序列插补，你通常希望避免从复杂的方法开始。这样做有几个原因：

短序列按定义是分散的，并且通常以大量出现。任何用于处理这些序列的插补方法都将高度受益于轻量级和操作高效。
在可用数据的样本量有限的情况下，简单的插补方法仍然可以提供合理的估计。
简单的插补方法也易于解释，因为它们通常来自你数据的基本统计属性。
采用从简单基线开始的概念。如果需要，只在复杂度尺度上进一步发展。

在我们的例子中，我们将比较线性插值与三次插值在缺失值插补 IMVs 方面的性能。简单来说：

线性插值将在两个已知点之间拟合一条直线来估计中间的缺失值。

三次插值将在一系列点之间拟合一个三次多项式函数来估计中间的缺失值。

下面的图展示了使用线性插值进行的插补。

fromtshelpers.metricsimportrmse_score,mae_score# Imputation with linear interpolationsubset_linear=subset.copy()subset_linear["PM 2.5"]=(subset_missing["PM 2.5"].interpolate(method="linear").tolist())plot_compare(subset_linear,subset_missing,...,plot_sup_title=f"MAE:{mae_score(subset,subset_linear,value='PM 2.5',verbose=False):.4f}|RMSE:{rmse_score(subset,subset_linear,value='PM 2.5',verbose=False):.4f}",)

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/0f182032c181cfdf3922a3f41ef39a63.png

线性插值示例

你可以看到我们引入了两个指标来评估插补的准确性：

平均绝对误差（MAE）：它定义为生成或预测值与实际值之间绝对差的平均值。

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/ad619ca3c9091c17c3c4186e6ef5b432.png

均方根误差（RMSE）：它定义为生成或预测值与实际值之间平方差的平均值的平方根。由于其二次性质，它对大误差给予更多的权重。

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/2511cb8103fb0806e899b8c1a3d7ef41.png

上面的公式显示了这两个指标作为 n 个观测值样本 yi 和 n 个相应的生成或预测值ŷ的函数[6][7]。

尽管有许多适合不同目的的准确度指标，但我们将坚持使用这两个指标作为我们的示例。我们鼓励你寻找更适合你数据集特性的指标。

在我们继续之前，关于 Pandas 三次插值实现的补充说明。

Pandas 的默认三次样条插值方法传递给 SciPy 模块中的[UnivariateSpline](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.UnivariateSpline.html)函数，并且没有考虑周围点的第一导数来拟合曲线。下面的图表展示了 Pandas 对三次样条插值的实现。

# Cubic interpolation - Pandas implementationsubset_spline=subset.copy()subset_spline["PM 2.5"]=(subset_missing["PM 2.5"].interpolate(method="spline",order=3).tolist())plot_compare(subset_spline,subset_missing,...)

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/14514f3e5fafdb310e5d602db0c42b38.png

使用 Pandas 进行三次样条插值

为了解决这个问题，我们将使用 SciPy 的[splrep](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.splrep.html)和[splev](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.splev.html#scipy.interpolate.splev)函数的组合。现在生成的结果将拟合一个曲线，该曲线评估相邻观测值的导数。

subset_spline=subset.copy()# Indexing missing valuesmissing_idx=np.where(subset_missing["PM 2.5"].isna())[0].tolist()# Indexing boundary values at t-2, t-1, t+1, and t+2boundary_idx=[]foridxinrange(min(missing_idx)-2,max(missing_idx)+3):ifidxnotinmissing_idx:boundary_idx.append(idx)boundary=subset_missing["PM 2.5"][boundary_idx].tolist()# Fitting cubic splinecubic_spline=interpolate.splrep(x=boundary_idx,y=boundary,k=3)imputed=interpolate.splev(missing_idx,cubic_spline)subset_spline["PM 2.5"][missing_idx]=imputed plot_compare(subset_spline,subset_missing,...)

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/4a3e37a7961dfca260fdfe128f53a5e1.png

使用 SciPy splrep 和 splev 进行三次样条插值

实现增加了复杂性，但结果现在符合预期。

实验方法

为了准确测量插补方法的有效性，我们将遍历：

a. 20 个实验子集。

b. 缺失序列的长度（对于 IMVs 为 1 或 2 小时）。

c. 缺失序列在子集中的位置。

对于 b.和 c.，我们将使用我们的create_missing()函数逐步增加缺失序列的填充。这将有效地将缺失样本从右向左滑动，同时覆盖整个实验集。下面的图像展示了这个过程。

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/11ed3d6ec5a4613540616dbd5e681237.png

人工生成缺失值的滑动窗口

在做这件事的时候，我们将存储每个迭代的 MAE 和 RMSE 指标，以便我们可以分析结果。为了简洁，下面的代码片段包含了这个实现的简化版本。在文章的末尾，你可以找到一个链接到包含原始代码的 jupyter 笔记本。

fromdarts.metricsimportmae,rmse imputation_results={}forstation_subsetinexp_subsets[station]:# Generating backwards missing intervalspaddings=[iforiinrange(2,len(exp_subsets[station][station_subset]))]forpaddinginpaddings:# Creating missing interval with length = 1 and length = 2forlengthin[1,2]:subset,subset_missing=create_missing(data=exp_subsets[station][station_subset],padding=padding,missing_length=length)### Cubic spline[...]# Fitting cubic splinecubic_spline=interpolate.splrep(x=boundary_idx,y=boundary,k=3)imputed=interpolate.splev(missing_idx,cubic_spline)subset_spline["PM 2.5"][missing_idx]=imputed# Evaluating cubic splineimputation_results[imputation_step]["RMSE_Cubic"].append(rmse(subset,subset_spline))imputation_results[imputation_step]["MAE_Cubic"].append(mae(subset,subset_spline))### Linear interpolationsubset_linear=subset.copy()subset_linear["PM 2.5"]=subset_missing["PM 2.5"].interpolate(method="linear",inplace=False).tolist()# Evaluating linear interpolationimputation_results[imputation_step]["RMSE_Linear"].append(rmse(subset,subset_linear))imputation_results[imputation_step]["MAE_Linear"].append(mae(subset,subset_linear))imputation_results[imputation_step]["Missing_Length"].append(length)

如果我们现在绘制结果，我们可以通过查看我们的准确率指标来比较两种插补方法。对于下面的图表，为了改善可视化，移除了超出三个标准差之外的异常值。

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/58fceda179cf956ecc4e964cb02ae2ca.png

立方和线性插值的插补指标

线性插值显示了更好的结果，并且正如预期的那样，对于单个观测值（1 小时）的插补，误差指标更小。为了具有统计学意义地得出结论，你可能会想要进行统计测试，例如配对样本的 Wilcoxon 符号秩检验（其中样本不是从正态分布中抽取的，如上面的结果所示）。

# Linear interpolation imputation of IMVsmaster_df.loc[master_df["seqMissing"]<=2,"PM 2.5"]=master_df.loc[master_df["seqMissing"]<=2,"PM 2.5"].interpolate()

通过使用线性插值来插补 IMVs，我们现在只剩下 4,175 个缺失值（从 4,958 个减少）。

长序列的插补

在我们开始之前，我们可以列出与长序列插补相关的一些挑战。

不确定性增加：缺失值的序列越长，插补值的认知不确定性就越大。使用预测模型来填补中断的时间序列中的长间隔并不罕见，统计模型往往严重依赖于最近的观测值（或滞后），特别是在单变量问题中。
对数据分布的影响：插补长序列的缺失值可以显著改变数据的分布，这对后续分析和解释尤其有害。
结果评估：对于长序列，测量插补方法的有效性具有挑战性。传统的准确率指标可能无法完全捕捉数据中时间动态和趋势的影响。

简而言之，不仅在长序列中插补缺失值更困难，而且评估结果的有效性也变得更加具有挑战性。

你可能已经注意到，现在是时候采用更复杂的插补算法了。虽然模型选择不是本文的重点，但你可以在以下选项中探索作为你的基准：

kNN：k-最近邻算法可以用来使用具有该特征值的 k 个最近邻的数据来插补缺失值。
MICE：缩写为通过链式方程的多变量插补，它是一个迭代算法，通过一系列预测模型来插补缺失数据。
FFT：可以使用快速傅里叶变换来通过从时间序列中提取潜在的周期性模式来进行预测。
SMA：简称为简单移动平均，它是一种算法，通过计算固定有限数量的先前可用值来生成后续缺失值的一组。

如果需要，可以深入研究更复杂的架构。以下是一些有趣的资源，如果您想了解更多关于最近提出的用于多元时间序列缺失值插补的内容：

填补空白：剑桥微软研究人员的 EDDI 多元时间序列缺失值插补
DeepMVI：多维时间序列缺失值插补 by IIT Bombay 研究人员在孟买，印度
M-RNN：多方向循环神经网络 by Yoon et al. 2017 用于缺失值估计
SAITS：基于自注意力的时间序列插补 by Du et al. 2023 (GitHub 仓库这里)。

这次，然而，我们将把我们的重点转移到问题的另一端。让我们介绍两种额外的评估我们的插补结果的方法。

插补对后续分析的影响。
插补前后分布的分析。

对后续分析的影响。

比较两种或多种插补方法的一个实用方法是评估它们对下游分析或建模任务的影响。

在我们的例子中，我们将假设我们的目标是预测温哥华市未来 PM 2.5 的浓度水平，为此，我们需要向模型提供符合最低完整性和连续性标准的数据。为了简化，我们将使用一个简单的 XGBoost 模型。

让我们先定义一个函数，将 Darts 的TimeSeries.plot方法包装起来，并添加一些额外的功能，以便在预测期间轻松查看我们的误差指标。

fromdartsimportTimeSeriesfromdarts.metricsimportmae,rmse# Auxiliar plotting functiondefdart_plot(train,val,plot_title="Time Series Train/Validation Sets",pred=None,pred_label="predicted"):plt.figure(figsize=(12,2.5))train.plot(label="training")val.plot(label="validation",linewidth=1,color="darkblue")ifpred:pred.plot(label=pred_label,linewidth=1.5,color="darkorange")plt.title(f"MAE:{mae(pred,val):.3f}| RMSE:{rmse(pred,val):.3f}",fontsize=10)plt.suptitle(plot_title,fontweight="bold",fontsize=11,y=1.04)else:plt.title(plot_title,fontweight="bold",fontsize=11)plt.legend()plt.xticks(fontsize=10)plt.yticks(fontsize=10)plt.ylabel("PM 2.5",fontsize=11)plt.xlabel("")sns.despine()plt.show()

使用TimeSeries.split_after()方法，我们将我们的实验子集分为训练集和验证集。

exp_subset=master_df[master_df['Station']=='Vancouver_Clark_Drive_PM25']exp_subset=exp_subset.loc[datetime(2019,8,8):datetime(2020,9,6),['PM 2.5']]validation_cutoff=pd.Timestamp("2019-11-29 19:00:00")exp_subset_series=TimeSeries.from_dataframe(exp_subset)train,val=exp_subset_series.split_after(validation_cutoff)

我们还需要定义模型将使用哪些滞后（或过去观测值）来进行预测。在我们的例子中，由于我们的数据具有季节性，我们将使用过去 48 小时内的 6 个滞后的倍数。

max_lag=48step=6[-iforiinrange(step,max_lag+step,step)]

这为我们提供了一组[-6, -12, -18, -24, -30, -36, -42, -48]滞后，用于预测。

我们将通过定义参数likelihood='poisson'来使用 XGBoost 模型的概率版本，以绘制预测区间。以下是 24 个时间步长（以小时为单位）的样本预测，从模型中采样 200 次。

fromdarts.modelsimportXGBModel model=XGBModel(lags=lags,output_chunk_length=2,likelihood='poisson',random_state=42)model.fit(train)pred=model.predict(n=24,num_samples=200)dart_plot(train[-48:],val[:24],pred=pred)

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/ecd0b97d2b39a5b36174f58765ad971e.png

基于 XGBoost 的概率预测示例

在我们的例子中，我们需要预测一个相对较高的 PM 2.5 指数期（数据集中所有四个站点的平均 PM 2.5 测量值为 6.2 µg/m3）。模型预计 11 月 30 日下午指数会增加，但不确定性很高，预测区间一直上升至 80 µg/m3。

将我们的重点重新回到原始问题上，我们想看看训练集中插补的缺失值对我们预测期的影响。为此，我们将定义一个函数，评估验证集中有多少点落在 95%预测区间内。

defevaluate_prediction_interval(pred:TimeSeries,val:TimeSeries,lower_quantile:float=0.05,upper_quantile:float=0.95)->np.ndarray:# Lower and Upper prediction interval serieslower_bound=pred.quantile_timeseries(lower_quantile).all_values().flatten()upper_bound=pred.quantile_timeseries(upper_quantile).all_values().flatten()# Validation setval_flatten=val.all_values().flatten()# Accuracy vectoraccuracy_vector=np.logical_and(val_flatten>=lower_bound,val_flatten<=upper_bound)returnaccuracy_vector

这个函数返回一个布尔数组，我们可以用它来确定预测区间的覆盖率。

coverage=sum(evaluate_prediction_interval(pred,val[:24]))/len(pred)print(f'Prediction Interval coverage:{coverage*100:.2f}%')

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/9adad1c8f8687caa936d9d977b988376.png

该指标告诉我们，验证集的 83.33%位于我们的预测区间内。在接下来的例子中，我们将从我们的训练集中移除 12 小时，看看两种插补方法（线性插补和 FFT 重组）如何影响这些结果。

_,train_missing=create_missing(data=train,value="PM 2.5",missing_length=12,padding=12,missing_index="end",)# Linear interpolationtrain_linear=train_missing.copy()train_linear["PM 2.5"]=(train_missing["PM 2.5"].interpolate(method="linear").tolist())plot_compare(train_linear[-72:],train[-72:])# FFT imputationtrain=TimeSeries.from_dataframe(train)fft_model=FFT(nr_freqs_to_keep=None)fft_model.fit(train[:(-24)])# -(missing_length + padding)fft_pred=fft_model.predict((12))# predicting n = missing_length stepstrain_fft=train_missing.pd_dataframe().copy()fft_pred_df=fft_pred.pd_dataframe().copy()train_fft.loc[train_fft["PM 2.5"].isna()]=fft_pred_df[["PM 2.5"]]plot_compare(train_fft[-72:],train[-72:])

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/890fcbf08d697749b07a145513702cfa.png

线性插补的一个 12 小时长序列示例

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/40019cd70a4e455a402be3e6b6851c35.png

FFT 插补法的一个 12 小时长序列示例

如我们所见，FFT 插补法无法捕捉到在人工生成的缺失区间内 PM 2.5 的峰值。因此，让我们看看这两种插补方法如何影响我们的 XGBoost 模型的预测。

model=XGBModel(lags=lags,output_chunk_length=2,likelihood='poisson',random_state=42)model.fit(train_linear)pred=model.predict(n=24,num_samples=200)dart_plot(train_linear[-48:],val[:24],pred=pred)

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/22579025c2b9863eafe40aba02317101.png

使用线性插补对训练集进行插补的预测结果

coverage=sum(evaluate_prediction_interval(pred,val[:24]))/len(pred)print(f'Prediction Interval coverage:{coverage*100:.2f}%')

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/2340f8578bef0cd909241892a10f9388.png

model=XGBModel(lags=lags,output_chunk_length=2,likelihood='poisson',random_state=42)model.fit(train_fft)pred=model.predict(n=24,num_samples=200)dart_plot(train_fft[-48:],val[:24],pred=pred)

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/daa12f5a626a72ac602c647f656fa398.png

使用 FFT 插补法对训练集进行插补的预测结果

coverage=sum(evaluate_prediction_interval(pred,val[:24]))/len(pred)print(f'Prediction Interval coverage:{coverage*100:.2f}%')

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/bc428e8d6620a05b7a378b4b18af5b47.png

预测区间的覆盖率清楚地表明，FFT 插补法的性能不佳正在影响我们的预测。

对数据分布的影响

另一种评估大间隔拉伸插补是否在合理值范围内的方法是查看插补前后的数据分布。

评估这一点的办法之一是简单地绘制原始训练集的直方图，并将其与插补后的数据进行比较。然而，这并不能给我们一个定量的指标来确定新值如何适应原始数据的分布。

为了解决这个问题，我们可以借鉴两个概率分布之间的相对熵测量的概念，也称为库尔巴克-莱布勒（KL）散度。它是一种统计距离，用于衡量一个概率分布 Q 与另一个参考 P 分布的差异[8]。它通过以下公式计算：

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/db4bc399191022293b4f24e2498a0c38.png

我们可以使用 KL 散度来比较插补前后训练集的分布。那么，让我们看看它对我们上一节中的例子有什么启示。

# KL Divergencefromscipy.specialimportrel_entr KL_Div_linear=sum(rel_entr(train[-48:].values(),train_linear[-48:].values()))KL_Div_fft=sum(rel_entr(train[-48:].values(),train_fft_imputed[-48:].values()))fig,axes=plt.subplots(1,2,figsize=(10,5),sharey=True)sns.histplot(train[-48:].pd_dataframe(),x="PM 2.5",label="OriginalnData",kde=True,color="navy",bins=10,ax=axes[0])sns.histplot(train_linear[-48:].pd_dataframe(),x="PM 2.5",label="LinearnInterpolation",kde=True,color="orange",bins=10,ax=axes[0])sns.histplot(train[-48:].pd_dataframe(),x="PM 2.5",label="OriginalnData",kde=True,color="navy",bins=10,ax=axes[1])sns.histplot(train_fft_imputed[-48:].pd_dataframe(),x="PM 2.5",label="FFTnImputation",kde=True,color="orange",bins=10,ax=axes[1])plt.suptitle(f"PM 2.5 Distribution",fontsize=13,fontweight="bold",y=0.95)axes[0].set(title=f"KL Divergence Linear:{KL_Div_linear[0]:.3f}")axes[0].legend()axes[1].set(title=f"KL Divergence FFT:{KL_Div_fft[0]:.3f}")axes[1].legend()plt.xlabel("PM 2.5")plt.tight_layout()sns.despine()plt.show()

https://github.com/OpenDocCN/towardsdatascience-blog-zh-2024/raw/master/docs/img/712808920e2eedf6ab5e5b5e1ece6297.png

原始数据与线性及 FFT 插补的分布比较

更高的 KL 散度意味着分布与原始训练数据之间的差异更大。因此，从上面的结果中，我们可以清楚地看到它如何惩罚 FFT 模型的糟糕插补。

结论

面对一个存在缺失值的时序数据集时，很容易一头扎进模型选择中。正如我们在最新章节中看到的，可能性众多，选项的数量随着我的写作而增长。

但与插补算法同样重要的是，以下这一点是根本性的：

彻底分析你的缺失值。你可能会经常从将你的插补方法分成不同大小的区间中受益的独立步骤中受益。
定义一个合适的评估策略。评估插补性能不应仅依赖于准确度指标。你总是可以从查看数据的分布及其对下游任务的影响中受益。

你可以在这里找到包含本文中展示的实现细节的 jupyter notebook。

喜欢这个故事吗？

你可以在 Medium 上关注我，了解更多关于数据科学、机器学习、可视化和数据分析的文章。

你还可以在 LinkedIn 和 X 上找到我，我在那里分享这些内容的简短版本。

每当 Erich Silva 发布时，都会收到电子邮件

参考文献

[1] Fang, C., & Wang, C. (2020). 时序数据插补：关于深度学习方法的综述。ArXiv。/abs/2011.11347

[2] Anna Richter, Jyotirmaya Ijaradar, Ulf Wetzker 等人. 多变量时序插补：关于可用方法的综述，重点关注混合 GANs。TechRxiv.2022 年 11 月 21 日。

[3] ReneBt (stats.stackexchange.com/users/195935/renebt)，数据生成过程（DGP）实际上是什么意思？URL（版本：2020-08-28）：stats.stackexchange.com/q/451230

[4] 服务，公民部。“不列颠哥伦比亚省数据目录。”不列颠哥伦比亚省。不列颠哥伦比亚省，2022 年 2 月 2 日。www2.gov.bc.ca/gov/content/data/bc-data-catalogue

[5] 刘辉，王宇，陈伟。(2020)。条件监测数据集中缺失值的三角插补。IET Generation, Transmission & Distribution，14(16)。doi.org/https://doi.org/10.1049/iet-gtd.2019.1446

[6] 维基媒体基金会。(2023 年 11 月 25 日)。平均绝对误差。维基百科。en.wikipedia.org/wiki/Mean_absolute_error

[7] 维基媒体基金会。(2024 年 1 月 8 日)。均方根偏差。维基百科。en.wikipedia.org/wiki/Root-mean-square_deviation

[8] 维基媒体基金会。(2024 年 1 月 10 日)。库尔布克-莱布勒散度。维基百科。en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence

[9] 加拿大不列颠哥伦比亚省公民服务部。“不列颠哥伦比亚省数据目录。”不列颠哥伦比亚省，2022 年 2 月 2 日。www2.gov.bc.ca/gov/content/data/bc-data-catalogue

[10] Engagement, G. C. 和 P. (2024 年 1 月 30 日)。不列颠哥伦比亚省开放政府许可。不列颠哥伦比亚省。www2.gov.bc.ca/gov/content/data/open-data/open-government-licence-bc

除非另有说明，所有图像均由作者提供。

处理时间序列中的间隔

目录