r - R中向量的子向量总和

Question

给定一个向量 x长度为 k，我想获得一个 k × k 矩阵 X哪里X[i,j]是 x[i] + ... + x[j] 的总和.我现在的做法是

set.seed(1)
x <- rnorm(10)

X <- matrix(0,10,10)
for(i in 1:10) 
  for(j in 1:10)
    X[i,j] <- sum(x[i:j])

#             [,1]       [,2]       [,3]      [,4]        [,5]       [,6]        [,7]      [,8]      [,9]      [,10]
# [1,]  -0.6264538 -0.4428105 -1.2784391 0.3168417  0.64634948 -0.1741189  0.31331014 1.0516348 1.6274162  1.3220278
# [2,]  -0.4428105  0.1836433 -0.6519853 0.9432955  1.27280329  0.4523349  0.93976395 1.6780887 2.2538700  1.9484816
# [3,]  -1.2784391 -0.6519853 -0.8356286 0.7596522  1.08915996  0.2686916  0.75612063 1.4944453 2.0702267  1.7648383
# [4,]   0.3168417  0.9432955  0.7596522 1.5952808  1.92478857  1.1043202  1.59174924 2.3300739 2.9058553  2.6004669
# [5,]   0.6463495  1.2728033  1.0891600 1.9247886  0.32950777 -0.4909606 -0.00353156 0.7347931 1.3105745  1.0051861
# [6,]  -0.1741189  0.4523349  0.2686916 1.1043202 -0.49096061 -0.8204684 -0.33303933 0.4052854 0.9810667  0.6756783
# [7,]   0.3133101  0.9397640  0.7561206 1.5917492 -0.00353156 -0.3330393  0.48742905 1.2257538 1.8015351  1.4961467
# [8,]   1.0516348  1.6780887  1.4944453 2.3300739  0.73479315  0.4052854  1.22575376 0.7383247 1.3141061  1.0087177
# [9,]   1.6274162  2.2538700  2.0702267 2.9058553  1.31057450  0.9810667  1.80153511 1.3141061 0.5757814  0.2703930
# [10,]  1.3220278  1.9484816  1.7648383 2.6004669  1.00518611  0.6756783  1.49614672 1.0087177 0.2703930 -0.3053884

但我不禁感到必须有一种更优雅的 R 方式(除了将其翻译成 Rcpp)。

Answer 1

这是另一种方法，它似乎比 OP 的 for 循环(因子 ~30)快得多，并且比当前存在的其他答案(因子 >=18)快:

n <- 5
x <- 1:5
z <- lapply(1:n, function(i) cumsum(x[i:n]))
m <- mapply(function(y, l) c(rep(NA, n-l), y), z, lengths(z))
m[upper.tri(m)] <- t(m)[upper.tri(m)]
m

#     [,1] [,2] [,3] [,4] [,5]
#[1,]    1    3    6   10   15
#[2,]    3    2    5    9   14
#[3,]    6    5    3    7   12
#[4,]   10    9    7    4    9
#[5,]   15   14   12    9    5

基准(向下滚动查看结果)

library(microbenchmark)
n <- 100
x <- 1:n

f1 <- function() {
  X <- matrix(0,n,n)
  for(i in 1:n) {
    for(j in 1:n) {
      X[i,j] <- sum(x[i:j])
    }
  }
  X
}

f2 <- function() {
  mySum <- function(i,j) sum(x[i:j])
  outer(1:n, 1:n, Vectorize(mySum))
}

f3 <- function() {
  matrix(apply(expand.grid(1:n, 1:n), 1, function(y) sum(x[y[2]:y[1]])), n, n)
}

f4 <- function() {
  z <- lapply(1:n, function(i) cumsum(x[i:n]))
  m <- mapply(function(y, l) c(rep(NA, n-l), y), z, lengths(z))
  m[upper.tri(m)] <- t(m)[upper.tri(m)]
  m
}

f5 <- function() {
  X <- diag(x)
  for(i in 1:(n-1)) {
    for(j in 1:(n-i)){
      X[j+i,j] <- X[j,j+i] <- X[j+i,j+i] + X[j+i-1,j]
    }  
  }
  X
}

microbenchmark(f1(), f2(), f3(), f4(), f5(), times = 25L, unit = "relative")
#Unit: relative
# expr      min       lq     mean   median       uq      max neval
# f1() 29.90113 29.01193 30.82411 31.15412 32.51668 35.93552    25
# f2() 29.46394 30.93101 31.79682 31.88397 34.52489 28.74846    25
# f3() 56.05807 53.82641 53.63785 55.36704 55.62439 45.94875    25
# f4()  1.00000  1.00000  1.00000  1.00000  1.00000  1.00000    25
# f5() 16.30136 17.46371 18.86259 17.87850 21.19914 23.68106    25

all.equal(f1(), f2())
#[1] TRUE
all.equal(f1(), f3())
#[1] TRUE
all.equal(f1(), f4())
#[1] TRUE
all.equal(f1(), f5())
#[1] TRUE

更新了 Neal Fultz 的编辑功能。