用于异常检测的具有缺失值的时间序列的 STL 分解

Posted 2023-03-12

【中文标题】用于异常检测的具有缺失值的时间序列的 STL 分解

【英文标题】：STL decomposition of time series with missing values for anomaly detection

【发布时间】：2012-08-17 00:09:57

【问题描述】：

我正在尝试检测时间序列气候数据中的异常值，其中一些观测值缺失。在网上搜索我发现了许多可用的方法。其中，stl 分解似乎很有吸引力，因为它去除了趋势和季节性成分并研究了其余部分。阅读STL: A Seasonal-Trend Decomposition Procedure Based on Loess，stl 在确定分配可变性的设置方面似乎很灵活，不受异常值的影响，并且尽管有缺失值也可以应用。但是，尝试在 R 中应用它，经过四年的观察并根据 http://stat.ethz.ch/R-manual/R-patched/library/stats/html/stl.html 定义所有参数，我遇到了错误：

时间序列包含内部 NAs

当na.action = na.omit，并且

序列不是周期性的或少于两个周期

当na.action = na.exclude。

我已经仔细检查了频率是否正确定义。我在博客中看到了相关问题，但没有找到任何可以解决此问题的建议。是否可以在缺少值的系列中应用 stl？我非常不愿意插入它们，因为我不想引入（并因此检测......）伪影。出于同样的原因，我不知道改用 ARIMA 方法有多么可取（如果缺失值仍然是一个问题）。

如果您知道在缺失值的系列中应用 stl 的方法，或者您认为我的选择在方法上不合理，或者您有更好的建议，请分享。我是该领域的新手，被大量（看似......）相关信息所淹没。

【问题讨论】：

你有多少缺失值？您可以尝试 zoo 包中的 na.approx 函数来插入丢失的数据。

一些额外的想法：显然，现在包 stlplus 可以处理 ts 分解中的缺失值。 Rbeast 可能是另一种选择；它是贝叶斯算法，不像stl只分解时间序列，Rbeast同时进行时间序列分解和变化点检测，允许缺失值。这是一个例子：library(Rbeast); co2[ sample(1:length(co2), 200) ]=NA; plot(beast(co2)).

【参考方案1】：

stl开头我们find

x <- na.action(as.ts(x))

不久之后

period <- frequency(x)
if (period < 2 || n <= 2 * period) 
    stop("series is not periodic or has less than two periods")

也就是说，stl 期望 x 在 na.action(as.ts(x)) 之后成为 ts 对象（否则 period == 1）。让我们先检查na.omit 和na.exclude。

很明显，在getAnywhere("na.omit.ts") 的末尾我们发现

if (any(is.na(object))) 
    stop("time series contains internal NAs")

这很简单，什么也做不了（na.omit 不排除 NAs 来自 ts 对象）。现在getAnywhere("na.exclude.default") 排除NA 观察，但返回类exclude 的对象：

    attr(omit, "class") <- "exclude"

这是一种不同的情况。如上所述，stl 期望 na.action(as.ts(x)) 为 ts，但 na.exclude(as.ts(x)) 属于 exclude 类。

因此，如果对NAs 排除感到满意，那么例如

nottem[3] <- NA
frequency(nottem)
# [1] 12
na.new <- function(x) ts(na.exclude(x), frequency = 12)
stl(nottem, na.action = na.new, s.window = "per")

有效。一般来说，stl 不适用于 NA 值（即使用 na.action = na.pass），它在 Fortran 中崩溃更深（参见完整源代码 here）：

z <- .Fortran(C_stl, ...

na.new 的替代品并不令人愉快：

na.contaguous - 查找时间序列对象中最长的连续非缺失值。 na.approx、na.locf 来自 zoo 或其他一些插值函数。 不确定这个，但是可以找到另一个用于 Python here 的 Fortran 实现。可以使用 Python 进行一些修改后从源代码安装 R，以防此模块确实允许缺失值。

正如我们在paper 中看到的那样，在调用stl 之前，没有一些简单的缺失值程序（比如在一开始就近似它们）可以应用于时间序列。因此，考虑到original implementation 相当冗长这一事实，我会考虑一些其他替代方案，而不是全新的实现。

更新：当NAs 可以是na.approx 和zoo 时，在许多方面都非常理想的选择，所以让我们检查它的性能，即将stl 的结果与完整的比较当有一定数量的NAs，使用na.approx 时的数据集和结果。我使用MAPE 作为准确性的衡量标准，但仅用于趋势，因为季节性分量和余数过零并且会扭曲结果。 NAs 的位置是随机选择的。

library(zoo)
library(plyr)
library(reshape)
library(ggplot2)
mape <- function(f, x) colMeans(abs(1 - f / x) * 100)

stlCheck <- function(data, p = 3, ...)
  set.seed(20130201)
  pos <- lapply(3^(0:p), function(x) sample(1:length(data), x))
  datasetsNA <- lapply(pos, function(x) data[x] <- NA; data)
  original <- data.frame(stl(data, ...)$time.series, stringsAsFactors = FALSE)
  original$id <- "Original"
  datasetsNA <- lapply(datasetsNA, function(x) 
    data.frame(stl(x, na.action = na.approx, ...)$time.series, 
               id = paste(sum(is.na(x)), "NAs"), 
               stringsAsFactors = FALSE))
  stlAll <- rbind.fill(c(list(original), datasetsNA))
  stlAll$Date <- time(data)
  stlAll <- melt(stlAll, id.var = c("id", "Date"))
  results <- data.frame(trend = sapply(lapply(datasetsNA, '[', i = "trend"), mape, original[, "trend"]))
  results$id <- paste(3^(0:p), "NAs")
  results <- melt(results, id.var = "id")
  results$x <- min(stlAll$Date) + diff(range(stlAll$Date)) / 4
  results$y <- min(original[, "trend"]) + diff(range(original[, "trend"])) / (4 * p) * (0:p)
  results$value <- round(results$value, 2)
  ggplot(stlAll, aes(x = Date, y = value, colour = id, group = id)) + geom_line() + 
    facet_wrap(~ variable, scales = "free_y") + theme_bw() +
    theme(legend.title = element_blank(), strip.background = element_rect(fill = "white")) + 
    labs(x = NULL, y = NULL) + scale_colour_brewer(palette = "Set1") +
    lapply(unique(results$id), function(z)
      geom_text(data = results, colour = "black", size = 3,
                aes(x = x, y = y, label = paste0("MAPE (", id, "): ", value, "%"))))

nottem，240 次观察

stlCheck(nottem, s.window = 4, t.window = 50, t.jump = 1)

co2，468 个观察结果

stlCheck(log(co2), s.window = 21)

mdeaths，72 个观察结果

stlCheck(mdeaths, s.window = "per")

在视觉上，我们确实在案例 1 和案例 3 中看到了一些趋势差异。但考虑到样本量 (72)，这些差异在案例 1 中非常小，在案例 3 中也令人满意。

【讨论】：

对，看来stl 在 R 中的实现是基于问题中链接的论文第 4 节（第 20 页）中讨论的 Fortran 实现，它避开了设计目标 4（第 5 页） ) 有利于速度（见第 21 页底部）。我希望有一个允许缺失值的解决方案（S 显然存在一个实现）。尽管我怀疑在大多数情况下，在 stl 之前通过 na.approx 传递输入之间的差异与处理缺失值的 stl 实现并没有太大差异，这确实是我的后备选项。

@user1935457，我更新了我的答案。显然，R 中 stl 的实现缺少一些部分，但正如我在回答中提到的，可能不值得重写它，特别是如果它真的不会像你说的那样有太大差异。只是考虑了一些具有完整数据和na.approx 的测试，我将用结果更新我的答案..

【参考方案2】：

意识到这是一个老问题，但我想我会更新，因为 R 中有一个更新的 stl 包，称为 stlplus。 Here is its homepage on github。您可以使用 install.packages("stlplus") 从 CRAN 安装它，或者使用 devtools::install_github("hafen/stlplus") 直接从 github 安装它。