需要帮助对齐python fill_between和数据点

Question

警告！我在这里问这样一个基础问题，真是太尴尬了！！！

大家好，我正在做一个项目，涉及到时间序列的总电子含量数据。我的目标是进行统计分析，找出由于地震引起的TEC异常。

我在参考这篇研究论文中的滑动IQR方法。但是，不知道为什么，我用论文中的公式得不到和论文一样的结果。所以我决定使用滚动平均 + 1.6标准差的方法。

问题是，当我使用ax.fill_between(x1 = index, y1 = ub, y2=lb)这个方法时，我的置信区间带比数据点多了一步。为了更好地理解，请看下面的图。

这是我目前正在做的：

df = DAEJ.copy()
window = 12
hourly = df.resample(rule="h").median()
hourly["ma"] = hourly["TEC"].rolling(window=window).mean()
hourly["hour"] = hourly.index.hour
hourly["std_err"] = hourly["TEC"].rolling(window=window).std()
hourly["ub"] = hourly["ma"] + (1.67* hourly["std_err"])
hourly["lb"] = hourly["ma"] - (1.67* hourly["std_err"])
hourly["sig2"] = hourly["TEC"].rolling(window=window).var()
hourly["kur"] = hourly["TEC"].rolling(window=window).kurt()
hourly["pctChange"] = hourly.TEC.pct_change(12, fill_method="bfill")
hourly = hourly.dropna()
dTEC = hourly[(hourly["TEC"] > hourly["ub"])]



fig, ax = plt.subplots(figsize=(12,4))
hourly["TEC"].plot(ax=ax, title="TEC Anomaly", label="Station: DAEJ")
ax.fill_between(x=hourly.index, y1=hourly["ub"], y2=hourly["lb"], color="red", label="Conf Interval", alpha=.4)
ax.legend()

这是我得到的结果！

从图中可以看到，数据和颜色带没有对齐。我知道计算一个12小时窗口的滚动平均会导致值向后移动12小时，但即使我去掉了前面的值，图形仍然没有对齐。

Answer 1

你需要把你的 rolling 窗口居中，方法是设置 center=True：

window = 12
hourly = df.resample(rule="h").median()
hourly["ma"] = hourly["TEC"].rolling(window=window, center=True).mean()
hourly["hour"] = hourly.index.hour
hourly["std_err"] = hourly["TEC"].rolling(window=window, center=True).std()
hourly["ub"] = hourly["ma"] + (1.67* hourly["std_err"])
hourly["lb"] = hourly["ma"] - (1.67* hourly["std_err"])
hourly["sig2"] = hourly["TEC"].rolling(window=window, center=True).var()
hourly["kur"] = hourly["TEC"].rolling(window=window, center=True).kurt()
hourly["pctChange"] = hourly.TEC.pct_change(12, fill_method="bfill")
hourly = hourly.dropna()
dTEC = hourly[(hourly["TEC"] > hourly["ub"])]

另外，你尝试做的事情应该用 shift 来完成：

# your original code
# ...

# then shift
ax.fill_between(x=hourly.index, y1=hourly["ub"].shift(-window//2),
                y2=hourly["lb"].shift(-window//2),
                color="red", label="Conf Interval", alpha=.4)

输出结果：