修改的 cumsum 函数答案

【问题标题】：Modified cumsum function修改的 cumsum 函数
【发布时间】：2011-11-01 05:54:19
【问题描述】：

下面是我正在处理的一段代码的简化版本（为了避免混淆，省略了很多额外的计算）。它只是cumsum 函数的修改形式。我不想重新发明轮子，那么这个功能已经存在了吗？如果不是，什么方案可以提供最好的速度？

#Set up the data   
set.seed(1)   
junk <- rnorm(1000000)   
junk1 <- rnorm(1000000)   
cumval <- numeric(1000000)   

#Initialize the accumulator   
cumval[1] <- 1   

#Perform the modified cumsum
system.time({   
for (i in 2:1000000) cumval[i] <- junk[i] + (junk1[i] * cumval[i-1])       
})   

#Plot the result
plot(cumval, type="l")

【问题讨论】：

您介意解释一下它是如何使用的吗？请注意，junk[1] 和 junk1[1] 永远不会在您的算法中使用...

标签： performance r

【解决方案1】：

这个算法非常适合compiler 包！

#Set up the data   
set.seed(1)   
junk <- rnorm(1000000)   
junk1 <- rnorm(1000000)

# The original code
f <- function(junk, junk1) {
  cumval <- numeric(1000000)
  cumval[1] <- 1
  for (i in 2:1000000) cumval[i] <- junk[i] + (junk1[i] * cumval[i-1])
  cumval
}
system.time( f(junk, junk1) ) # 4.11 secs

# Now try compiling it...
library(compiler)
g <- cmpfun(f)
system.time( g(junk, junk1) ) # 0.98 secs

...所以知道这个算法是否以任何方式“典型”会很有趣 - 在这种情况下，编译器可能会针对这种情况进行更多优化...

【讨论】：

【解决方案2】：

它更快，但没有给出正确的结果。运行这个

set.seed(1)

N <- 10

junk  <- rnorm(N)

junk1 <- rnorm(N)

cumval <- numeric(N)
cumval.1 <- numeric(N)
cumval[1] <- 1

for( i in 2:N ) cumval[i] <- junk[i] + junk1[i]*cumval[i-1]
cumval

cumval.1 <- cumsum( junk[-1] + (junk1[-1] * cumval.1[-N]) ) 

cumval.1

你会发现 cumval 和 cumval.1 的长度甚至都不一样。

需要重写递归关系。我没有看到将循环转换为非循环公式的方法。

【讨论】：

这似乎不是答案，而是@Seth 对答案的评论。但是，您的观点似乎是正确的。

【解决方案3】：

考虑 cumval[5]。将 j[] 用于 junk 并将 jk[] 用于 junk1 并省略 * 符号，其扩展为：

j[5] +jk[5]j[4] + jk[5]jk[4]j[3] + jk[5]jk[4]jk[3]j[2] + jk[5]jk[4]jk[3]jk[2]

模式表明这可能是（接近于？）第 5 项的表达式：

    sum(  j[1:5] * c(1, Reduce("*" , rev(jk[2:5]), accumulate=TRUE) )

【讨论】：

这是修改后的“cumsum”的另一个版本（请注意 Gabor G 的回应），其中涉及“减少”。 r.789695.n4.nabble.com/…我还没有弄清楚为什么会这样，但是有了你的解决方案，有些东西可能会凝固。