Модифицированная функция cumsum - PullRequest
1 голос
/ 08 июля 2011

Ниже приведена упрощенная версия сегмента кода, над которым я работаю (для избежания путаницы оставлено множество дополнительных вычислений).Это просто модифицированная форма функции 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")    

Ответы [ 3 ]

1 голос
/ 09 июля 2011

Рассмотрим cumval [5]. Используя j [] для мусора и jk [] для мусора1 и пропуская символы *, его расширение будет:

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) )
1 голос
/ 01 ноября 2011

Этот алгоритм идеально подходит для пакета 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

... так что было бы интересно узнать, является ли этот алгоритм каким-либо образом "типичным" - в этом случае компилятор может быть еще более оптимизирован для подобных ситуаций ...

1 голос
/ 08 июля 2011

Это быстрее, но не дает правильных результатов.Запустите

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 даже не имеют одинаковую длину.

Нужно переписать отношение повторения.Я не вижу способа преобразовать повторение в непериодическую формулу.

...