计算R中的移动新近加权平均值答案

【问题标题】：Calculate moving recency-weighted mean in R计算R中的移动新近加权平均值
【发布时间】：2012-09-05 15:28:45
【问题描述】：

在给定时间（天）和完成位置（位置）的情况下，我想计算一匹马参加的一系列比赛的移动新近加权平均完成位置。这样的统计数据在handicapping 中很有用。

目前，我正在使用“循环内循环”方法。有没有更快或更优雅的 R 语言方法来解决这个问题？

#
# Test data
#

day <- c(0, 6, 10, 17, 21, 26, 29, 31, 34, 38, 41, 47, 48, 51, 61)
pos <- c(3, 5, 6, 1, 1, 3, 4, 1, 2, 2, 2, 6, 4, 5, 6)
testdata <- data.frame(id = 1, day = day, pos = pos, wt.pos = NA)

#
# No weight is given to observations earlier than cutoff
#

cutoff <- 30

#
# Rolling recency-weighted mean (wt.pos)
#

for(i in 2:nrow(testdata)) {
  wt <- numeric(i-1)
  for(j in 1:(i-1))
    wt[j] <- max(0, cutoff - day[i] + day[j] + 1)
    if (sum(wt) > 0)
      testdata$wt.pos[i] <- sum(pos[1:j] * wt) / sum(wt)
}

> testdata

   id day pos   wt.pos
1   1   0   3       NA
2   1   6   5 3.000000
3   1  10   6 4.125000
4   1  17   1 4.931034
5   1  21   1 3.520548
6   1  26   3 2.632911
7   1  29   4 2.652174
8   1  31   1 2.954128
9   1  34   2 2.436975
10  1  38   2 2.226891
11  1  41   2 2.119048
12  1  47   6 2.137615
13  1  48   4 3.030534
14  1  51   5 3.303704
15  1  61   6 4.075000

【问题讨论】：

标签： r

【解决方案1】：

我会去

# Calculate `wt` for all values of `i` in one go
wt <- lapply(2:nrow(testdata), function(i)
    pmax(0, cutoff - day[i] + day[1:(i-1)] + 1))

# Fill in the column
testdata$wt.pos[-1] <- mapply(
    function(i, w) if(sum(w) > 0) sum(pos[1:i]*w)/sum(w) else NA,
    1:(nrow(testdata)-1), wt)

请注意，通过同时计算 max 的所有值的 max 的第二个参数，我们对计算进行了矢量化处理，这将速度提高了多个数量级。

我发现没有简单的方法来矢量化外循环和 if 情况（除了用 C 重写它似乎有点过分），但 lapply、mapply 和类似的仍然比 for 快循环。

【讨论】：

您在笔记中将numeric 与as.numeric 混淆了
当应用于包含 5,000 个“testdata”副本的数据帧时，此解决方案比使用循环内循环方法的解决方案快约 20%。然而，使用 mapply 和保留外部“for”循环之间似乎几乎没有任何区别。
真的吗？这很奇怪，但话说回来，R 在计算速度方面有点不可预测，所以最好的方法是总是像你一样验证它。然而，向量化总是比循环更好，因为您在底层 C 或 fortran 代码中完成所有工作。也许像我一样计算wt 的批处理，但保留外部for 循环是最好的解决方案。

【解决方案2】：

此版本演示如何计算 1 个或多个变量（例如，完成位置、速度等级等）和 1 个或多个主题（马）的移动新近加权平均值。

library(plyr)

day <- c(0, 6, 10, 17, 21, 26, 29, 31, 34, 38, 41, 47, 48, 51, 61)
pos <- c(3, 5, 6, 1, 1, 3, 4, 1, 2, 2, 2, 6, 4, 5, 6)
dis <- 100 + 0.5 * (pos - 1)
testdata1 <- data.frame(id = 1, day = day, pos = pos, dis = dis)
day <- c(0, 4, 7, 14, 22, 23, 31, 38, 42, 47, 52, 59, 68, 69, 79)
pos <- c(1, 3, 2, 6, 4, 5, 2, 1, 4, 5, 2, 1, 5, 5, 2)
dis <- 100 + 0.5 * (pos - 1)
testdata2 <- data.frame(id = 2, day = day, pos = pos, dis = dis)
testdata <- rbind(testdata1, testdata2)

# Moving recency-weighted mean
rollmean <- function(day, obs, cutoff = 90) {
  obs <- as.matrix(obs)
  wt <- lapply(2:nrow(obs), function(i)
    pmax(0, cutoff - day[i] + day[1:(i-1)] + 1))
  wt.obs <- lapply(1:(nrow(obs)-1), FUN =
    function(i)
      if(sum(wt[[i]]) > 0) {
        apply(obs[1:i, , drop = F] * wt[[i]], 2, sum) / sum(wt[[i]])
      } else {
        rep(NA, ncol(obs))
      }
  )
  answer <- rbind(rep(NA, ncol(obs)), do.call(rbind, wt.obs))
  if (!is.null(dimnames(answer)))
    dimnames(answer)[[2]] <- paste("wt", dimnames(answer)[[2]], sep = ".")
  return(answer)
}

x <- dlply(testdata, .(id), .fun =
  function(DF) rollmean(DF$day, DF[, c("pos", "dis"), drop = F])
)
y <- do.call(rbind, x)

【讨论】：