时间序列缺失值处理之滑动平均（移动平均）插补模型（内含python代码）

---

前言

玩弄过机器学习的各位想必都或多或少有遇到过数据集中含有缺失值的情况。如果数据集是一个多维的数据（矩阵形式）,那可以考虑KNN或missforest这一类方法，但如果是一个一维向量表示的时间序列呢？这里给大家提供一个简单有效的方法——滑动平均插补（Moving Average Imputation）。你可别瞧不起它，不信你接着看…

---

一、滑动平均插补的基本思路

移动平均插补其实就是用一个滑动窗口在时间序列上滑动处理缺失值，比如下面这幅图

虚线框表示滑动窗口，框中心点为处理点（注意只有缺失的位置才做处理，不缺的直接跳过）。绿色的部分为滑动窗口处理过的部分，红色的部分为未处理的部分。它的思路非常简单，核心就在于滑动窗口的设计。理解了一个滑动窗口内的缺失值的插补思路，其实就等于理解了滑动平均插值。

1.1滑动窗口

• 固定长度的滑动窗口：处理点左右各取K个数，构成大小为2K+1的滑动窗口，对窗口内n个非缺失元素求和后取均值来填补缺失值，这里2≤n≤2K。当n取值过大时，可能出现窗口中非缺失元素的个数＜n的情况，从而导致无法填补缺失值。这时候就要用到自适应长度的滑动窗口。
• 自适应长度的滑动窗口：与固定长度的滑动窗口思想类似，不同的是，如果窗口内的非缺失元素<n，则扩大窗口（左边不够就扩大左边，右边不够就扩大右边），直到窗口内的元素满足n个为止。

1.2窗口的滑动方式

• 单向：往一个方向滑动，一般为正向滑动。
• 双向交替：正反交替滑动。
• 双向平均：做一次正向滑动和一次反向滑动后取两个结果的平均。

通过组合不同的滑动窗口和滑动方式就能得到不同的滑动平均模型。

二、缺失值生成

在实验开始前，我们先人为的生成数据缺失以便考察模型对这些缺失的插补效果。下面是生成了两种类型的缺失：

2.1随机点缺失

往原始序列中随机添加nan来表示缺失，代码如下：

import numpy as np

def random_missing(Time_series, missing_rate):
    mask = np.round(np.random.rand(len(Time_series)) + 0.5 - missing_rate) == 0
    Time_series[mask] = np.nan
    return Time_series

2.2随机段缺失

有的时候数据可能会出现呈现出一段一段的缺失。为了生成这种缺失，可以在随机点缺失的代码上稍微调整。首先把向量数据重塑为矩阵，然后随机挑选一些行并让这些行全部缺失，最后再将这个矩阵展开为向量，就得到了长度不固定的段缺失。
代码如下：

import numpy as np

def random_missing2(Time_series, missing_rate, miss_gap=2):
    if len(Time_series) % miss_gap:
        print('无法重构为指定形状')
        raise Exception
    else:
        length = len(Time_series) // miss_gap
    M = np.round(np.random.rand(length) + 0.5 - missing_rate).reshape(-1, 1) == 0
    mask = np.tile(M, (1, miss_gap)).ravel()
    Time_series[mask] = np.nan
    return Time_series

三、缺失值插补

3.1插补模型

把前面介绍的2种滑动窗口和3种滑动方式进行排列组合得到以下6中滑动平均模型，分别是：

• 固定长度的滑动窗口+单向滑动（MA1）：

def MA1(series, neighbor=1):
    # 给时间序列两端都填充数值（这里填充的是未缺失部分的均值）以方便处理两端
    series = np.pad(series, (neighbor, neighbor), 'constant', constant_values=np.nanmean(series))
    for i in range(neighbor, len(series) - neighbor):
        if np.isnan(series[i]):
            series[i] = np.nanmean(series[i - neighbor:i + neighbor])
    return series[neighbor:len(series) - neighbor]

由于这种方式可能会有部分缺失值无法插补，故不采用这种方式进行实验，但后续的两部分代码会基于上述代码来编写

• 固定长度的滑动窗口+双向交替滑动（MA2）：双向交替滑动解决了单向滑动有时无法处理完所有缺失值的不足。

def MA2(series, neighbor=1):
    while np.isnan(series).any():
        # 正向滑动
        ts_forward = MA1(series, neighbor=1)
        # 反向滑动
        ts_backward = MA1(ts_forward[::-1], neighbor=1)[::-1]
        series = ts_backward
    return series

• 固定长度的滑动窗口+双向平均（MA3）：

def MA3(series, neighbor=1):
    while np.isnan(series).any():
        # 正向滑动
        ts_forward = MA1(series, neighbor=1)
        # 反向滑动
        ts_backward = MA1(series[::-1], neighbor=1)[::-1]
        # 都有值的部分计算均值
        temp_series = (ts_forward + ts_backward) / 2
        # 只有一个值的就用那个值来插补
        half_obs_idx = np.bitwise_xor(np.isnan(ts_forward), np.isnan(ts_backward))
        ts_forward[np.isnan(ts_forward)] = 0
        ts_backward[np.isnan(ts_backward)] = 0
        temp_series[half_obs_idx] = (ts_forward + ts_backward)[half_obs_idx]

        series = temp_series
    return series

• 自适应长度的滑动窗口+单向滑动（MA4）：

def MA4(series, neighbor=1):
    series = np.pad(series, (neighbor, neighbor), 'constant', constant_values=np.nanmean(series))
    for i in range(neighbor, len(series) - neighbor):
        left_neighbor = right_neighbor = neighbor
        if np.isnan(series[i]):
            while np.count_nonzero(~np.isnan(series[i - left_neighbor:i])) < neighbor:
                left_neighbor += 1
            while np.count_nonzero(~np.isnan(series[i:i + right_neighbor])) < neighbor:
                right_neighbor += 1
            series[i] = np.nanmean(series[i - left_neighbor:i + right_neighbor])
    return series[neighbor:len(series) - neighbor]

• 自适应长度的滑动窗口+双向交替：没必要这样，因为自适应的滑动窗口只需单向滑动一次就可以填补完所有的缺失值。
• 自适应长度的滑动窗口+双向平均（MA5）：

def MA5(series, neighbor=1):
    while np.isnan(series).any():
        # 正向滑动
        ts_forward = MA4(series, neighbor=1)
        # 反向滑动
        ts_backward = MA4(series[::-1], neighbor=1)[::-1]
        # 取均值
        temp_series = (ts_forward + ts_backward) / 2

        series = temp_series
    return series

3.2实验结果

先生成个简单的正弦波进行缺失值插补实验，数据如下：

t = np.arange(0, 5, 1/100)
TS_1 = np.sin(2 * np.pi * t)

通过对比MA2-MA5这4个模型的插补值与真实值的均方根误差（RMSE）来评估模型的性能，需要注意的未缺失的部分不参与误差计算哦，代码及结果如下：

import matplotlib.pyplot as plt
np.random.seed(1)	# 设置随机种子以便复现结果
# *****************  两类缺失+缺失率30%  *****************
TS_2 = [random_missing(TS_1.copy(), missing_rate=0.3), random_missing2(TS_1.copy(), missing_rate=0.3, miss_gap=4)]

# *****************  插补结果对比  *****************
plt.figure(figsize=[10, 4])
for ts in range(len(TS_2)):
    maks = np.isnan(TS_2[ts])
    IMP = [MA2(TS_2[ts], neighbor=1), MA3(TS_2[ts], neighbor=1), MA4(TS_2[ts], neighbor=1), MA5(TS_2[ts], neighbor=1)]
    RMSE = []
    for i in IMP:
        RMSE.append(round(np.sqrt(np.sum((i[maks] - TS_1[maks]) ** 2) / len(i[maks])), 3))
    # 创建数据
    x = ['MA2', 'MA3', 'MA4', 'MA5']
    plt.subplot(1, 2, ts+1)
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签SimHei
    plt.rcParams['axes.unicode_minus'] = False
    for a, b in zip(x, RMSE):
        plt.text(a, b, str(b), ha='center', va='bottom')
    plt.bar(x, RMSE, color=['red', 'blue', 'green', 'yellow'])
    plt.title('缺失类型' + str(ts+1) + '误差对比')
    plt.ylabel('RMSE')
plt.show()

感觉还行，下面用两幅图来展示插补效果。这里仅选择误差最小的模型（MA5）所得结果进行展示


怎么说呢？感觉不对劲起来了呢，赶紧换个真实数据集继续搞。这里用的是airline-passengers数据集，我已经事先下好了数据集，当然也可以在线下载（使用被注释的那3行代码），但是需要科学上网。

# import pandas as pd
# df = pd.read_csv('https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv', header=0, index_col=0, parse_dates=True)
# TS_1 = df.values().astype(float).ravel()
TS_1 = np.load('airline-passengers.npy').astype(float).ravel()

---

一样的，先对比RMSE，后看插补效果。

结果不堪入目o(╥﹏╥)o ~~~

总结

最后做一下简单的总结吧，滑动平均仅仅能够处理间隔很小的缺失，缺失间隔稍微大一点立马就拉B胯！！！这里使用的是简单的滑动窗口，当然你还可以进一步尝试加权滑动窗口，详见参考文献，但我觉得应该都大差不差。写这篇文章呢是因为前段时间有个审稿人让我对比这个方法，真就挺无语的，不过弄都弄了咱就贴上来呗。

参考文献

ImputeTS: time series missing value imputation in R