はじめ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")
つまり、(そうでない場合) の後にオブジェクトであることがstl期待xされます。まずはチェックしてみましょう。tsna.action(as.ts(x))period == 1na.omitna.exclude
明らかに、getAnywhere("na.omit.ts")私たちが見つけた最後に
if (any(is.na(object)))
stop("time series contains internal NAs")
これは簡単で、何もできません (オブジェクトからna.omit除外されません)。観測を除外するようになりましたが、クラスのオブジェクトを返します:NAstsgetAnywhere("na.exclude.default")NAexclude
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機能せず、Fortran ではより深くクラッシュします (完全なソース コードはこちらを参照してください)。NAna.action = na.pass
z <- .Fortran(C_stl, ...
の代替na.newは楽しいものではありません：
na.contaguous- 時系列オブジェクト内の非欠損値の最長連続ストレッチを見つけます。na.approx、na.locffromzooまたはその他の補間関数。- これについてはよくわかりませんが、別の Fortran 実装が Python 用にここにあります。このモジュールが実際に欠損値を許可する場合は、Python を使用して、いくつかの変更後にソースから R をインストールすることができます。
論文でわかるように、 を呼び出す前に時系列に適用できる欠損値の単純な手順 (最初にそれらを概算するなど) はありませんstl。したがって、元の実装が非常に長いという事実を考慮して、まったく新しい実装以外の代替案について考えます。
更新:がからのものであるNAs可能性がある場合、多くの面で非常に最適な選択であるため、そのパフォーマンスを確認してみましょう。つまり、 の結果を完全なデータセットと比較し、いくつかのを使用して を使用した場合の結果を比較します。私は正確さの尺度としてMAPEを使用していますが、季節成分と残りがゼロを横切り、結果を歪めるため、トレンドのみに使用しています. の位置はランダムに選択されます。na.approxzoostlNAsna.approxNAs
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 では満足のいくものです。