R runmed 運行中位數 – 穩健散點圖平滑

說明

計算奇數跨度的運行中位數。這是可能的“最穩健”散點圖平滑。為了提高效率(和曆史原因)，您可以使用兩種不同的算法之一來給出相同的結果。

用法

runmed(x, k, endrule = c("median", "keep", "constant"),
       algorithm = NULL,
       na.action = c("+Big_alternate", "-Big_alternate", "na.omit", "fail"),
       print.level = 0)

參數

x	數值向量，要平滑的 ‘dependent’ 變量。
k	中值窗口的整數寬度；一定是奇數。 Turlach 的默認值是 k <- 1 + 2 * min((n-1)%/% 2, ceiling(0.1*n)) 。使用 k = 3 進行 ‘minimal’ 穩健平滑，消除孤立的異常值。
endrule	字符串，指示應如何處理(數據)開頭和結尾的值。可以縮寫。可能的值為： "keep" 在兩端保留第一個和最後一個 k_2 值，其中 k_2 是 half-bandwidth k2 = k %/% 2 ，即 y[j] = x[j] 對應 j \in \{1,\ldots,k_2; n-k_2+1,\ldots,n\} ； "constant" 將 median(y[1:k2]) 複製到第一個值，並類似地複製最後一個值，使平滑末端保持不變； "median" 默認情況下，通過使用隨後較小帶寬的對稱中位數來平滑末端，但對於應用 Tukey 強大的 end-point 規則的第一個和最後一個值，請參閱 smoothEnds 。
algorithm	字符串(部分匹配 "Turlach" 或 "Stuetzle" )或默認 NULL ，指定應應用哪種算法。默認選擇取決於n = length(x) 和k，其中"Turlach" 將用於更大的問題。
na.action	確定 x 中 NA 或 NaN 情況下的行為的字符串，(部分匹配)其中之一 "+Big_alternate" 在這裏，x 中的所有 NA 首先替換為交替的 \pm B，其中 B 是 “Big” 編號(使用 2B < M* ，其中 M*= .Machine $ double.xmax )。替換值為 “from left” (+B, -B, +B, \ldots) ，即以 "+" 開頭。 "-Big_alternate" 與 "+Big_alternate" 幾乎相同，隻是以 -B ( "-Big..." ) 開頭。 "na.omit" 結果與 runmed(x[!is.na(x)], k, ..) 相同。 "fail" x 中存在 NA 將引發錯誤。
print.level	整數，表示算法的詳細程度；普通用戶很少應該改變。

細節

除了最終值之外，結果 y = runmed(x, k) 僅包含 y[j] = median(x[(j-k2):(j+k2)]) ( k = 2*k2+1 )，計算效率非常高。

這兩種算法在內部完全不同：

"Turlach"

是由 Berwin Turlach 實現的 Härdle-Steiger 算法(參見參考文獻)。使用樹算法，確保性能O(n \log k)，其中n = length(x)是漸近最優的。

"Stuetzle"

是(較舊的)Stuetzle-Friedman 實現，當一個觀察進入和一個離開平滑窗口時，它利用中值更新。雖然這與 O(n \times k) 一樣執行，漸近速度較慢，但對於小型 k 或 n 來說，速度要快得多。

請注意，兩種算法(和 smoothEnds() 實用程序)現在 “work” 也適用於 x 包含非有限條目( \pm Inf 、 NaN 和 NA )：

"Turlach"

…………

"Stuetzle"

目前隻需應用底層數學庫即可工作(‘庫’) 非有限數的算術；將來這可能會發生變化。

目前 long vectors 僅支持 algorithm = "Stuetzle" 。

值

與 x 長度相同的平滑值向量，其中 attr ibute k 包含 (‘oddified’) k 。

例子

require(graphics)

utils::example(nhtemp)
myNHT <- as.vector(nhtemp)
myNHT[20] <- 2 * nhtemp[20]
plot(myNHT, type = "b", ylim = c(48, 60), main = "Running Medians Example")
lines(runmed(myNHT, 7), col = "red")

## special: multiple y values for one x
plot(cars, main = "'cars' data and runmed(dist, 3)")
lines(cars, col = "light gray", type = "c")
with(cars, lines(speed, runmed(dist, k = 3), col = 2))


## nice quadratic with a few outliers
y <- ys <- (-20:20)^2
y [c(1,10,21,41)] <- c(150, 30, 400, 450)
all(y == runmed(y, 1)) # 1-neighbourhood <==> interpolation
plot(y) ## lines(y, lwd = .1, col = "light gray")
lines(lowess(seq(y), y, f = 0.3), col = "brown")
lines(runmed(y, 7), lwd = 2, col = "blue")
lines(runmed(y, 11), lwd = 2, col = "red")

## Lowess is not robust
y <- ys ; y[21] <- 6666 ; x <- seq(y)
col <- c("black", "brown","blue")
plot(y, col = col[1])
lines(lowess(x, y, f = 0.3), col = col[2])


lines(runmed(y, 7),      lwd = 2, col = col[3])
legend(length(y),max(y), c("data", "lowess(y, f = 0.3)", "runmed(y, 7)"),
       xjust = 1, col = col, lty = c(0, 1, 1), pch = c(1,NA,NA))

## An example with initial NA's - used to fail badly (notably for "Turlach"):
x15 <- c(rep(NA, 4), c(9, 9, 4, 22, 6, 1, 7, 5, 2, 8, 3))
rS15 <- cbind(Sk.3 = runmed(x15, k = 3, algorithm="S"),
              Sk.7 = runmed(x15, k = 7, algorithm="S"),
              Sk.11= runmed(x15, k =11, algorithm="S"))
rT15 <- cbind(Tk.3 = runmed(x15, k = 3, algorithm="T", print.level=1),
              Tk.7 = runmed(x15, k = 7, algorithm="T", print.level=1),
              Tk.9 = runmed(x15, k = 9, algorithm="T", print.level=1),
              Tk.11= runmed(x15, k =11, algorithm="T", print.level=1))
cbind(x15, rS15, rT15) # result for k=11  maybe a bit surprising ..
Tv <- rT15[-(1:3),]
stopifnot(3 <= Tv, Tv <= 9, 5 <= Tv[1:10,])
matplot(y = cbind(x15, rT15), type = "b", ylim = c(1,9), pch=1:5, xlab = NA,
        main = "runmed(x15, k, algo = \"Turlach\")")
mtext(paste("x15 <-", deparse(x15)))
points(x15, cex=2)
legend("bottomleft", legend=c("data", paste("k = ", c(3,7,9,11))),
       bty="n", col=1:5, lty=1:5, pch=1:5)

作者

Martin Maechler maechler@stat.math.ethz.ch, based on Fortran code from Werner Stuetzle and S-PLUS and C code from Berwin Turlach.

參考

Härdle, W. and Steiger, W. (1995) Algorithm AS 296: Optimal median smoothing, Applied Statistics 44, 258-264. doi:10.2307/2986349.

Jerome H. Friedman and Werner Stuetzle (1982) Smoothing of Scatterplots; Report, Dep. Statistics, Stanford U., Project Orion 003.

也可以看看

smoothEnds 實現 Tukey 的端點規則，默認從 runmed(*, endrule = "median") 調用。 smooth 的複合平滑器使用運行中位數 3。