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。