我正在尝试检测一系列气候数据中的异常值，其中一些观测值缺失。在网上搜索后，我发现了许多可用的方法。其中，stl分解似乎很有吸引力，因为它可以消除趋势和季节性组件并研究余项。阅读 STL：基于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开始时，我们找到
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")
也就是说，在na.action(as.ts(x))之后，stl期望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 不会从 ts 对象中排除 NAs）。现在，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深处崩溃（请参见完整的源代码此处）：
z <- .Fortran(C_stl, ...
与na.new的替代方法并不令人愉悦:
na.contiguous - 查找时间序列对象中最长连续的非缺失值。na.approx，na.locf来自zoo或其他插值函数。正如我们可以在论文中看到的那样，在调用stl之前，没有一些简单的处理缺失值的过程（例如在一开始近似它们），可以应用于时间序列。鉴于原始实现相当冗长，我会考虑一些其他替代方法而不是全新的实现。
更新:当存在NAs时，一个在许多方面都相当优化的选择可能是zoo中的na.approx，因此让我们检查其性能，即使用na.approx比较具有完整数据集和具有一些NAs的结果的stl。我使用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之间趋势上的一些差异。但在案例1中这些差异相当小，在考虑样本量（72）的情况下，在案例3中也是可以接受的。
stl实现是基于问题链接中讨论的第4节（第20页）中提到的Fortran实现，该实现放弃了设计目标4（第5页），而选择速度（请参见第21页底部）。我希望有一种解决方案可以允许缺失值（S中显然存在这样的实现）。虽然我怀疑，在大多数情况下，通过na.approx将输入传递给stl与处理缺失值的stl实现之间的差异不会太大，这也是我的备选方案。 - user1935457stl缺少一些部分，但正如我在回答中提到的那样，如果真的没有太大区别，重写它可能不值得。只是考虑了一些使用完整数据和na.approx进行测试的情况，我将在我的回答中更新结果。 - Julius Vainora我知道这个问题比较老，但是有一个新的R语言的stl包叫做stlplus。 在GitHub上可以找到它的主页。你可以通过CRAN使用install.packages("stlplus")或者直接从GitHub上使用devtools::install_github("hafen/stlplus")来安装该软件包。
zoo包中的na.approx函数来插值你的缺失数据。 - Fernandostlplus包可以处理时间序列分解中的缺失值。Rbeast可能是另一个选择；它是一种贝叶斯算法，不同于只能分解时间序列的stl，Rbeast可以同时进行时间序列分解和变点检测，并允许存在缺失值。这里有一个例子：library(Rbeast);co2[ sample(1:length(co2), 200) ]=NA;plot(beast(co2))。 - zhaokg