кумулятивное применение функции по размеру в R

Question

Я хотел бы ускорить то, что показано в псевдо-коде ниже, до максимально возможного в R (векторизованный или любой метод, который быстрее, чем простой цикл цикла).

визуализации У меня есть 4-мерный массив A (заполненный произвольно 1 раз в качестве примера):

A = array(runif(nx*ny*nz*nt), c(nx,ny,nz,nt))

, и я хочу сделать это для цикла быстрее (заполнить output массив, который имеет более высокое второе измерение в кумулятивных модах от предыдущего значения ... больше как совокупного продукт второго измерения массива входных A:

output = array(1, c(nx, ny+1, nz, nt)) 
for (x in 1:nx) 
{ 
    for (z in 1:nz) 
    { 
    for (t in 1:nt) 
    { 
     for (y in 2:(ny+1)) 
     { 
     output[x,y,z,t] = output[x,y-1,z,t] * (1 - A[x,y-1,z,t]) 
     } 
    } 
    } 
} 

Как я могу сделать это быстрее? используя apply()? или какой-то умный кумулятивный продукт с abind() в конце?

Answer 1

Вы, конечно, можете использовать apply и cumprod, чтобы получить тот же результат (aperm необходимо, так как результаты функции, вызываемой применить конец в первом измерении):

output1 <- aperm(apply(A,c(1,3,4),function(v) cumprod(1-v)),c(2,1,3,4)) 

Сравнивая результат output в различия все очень близко к .Machine$double.eps:

> max(abs(output[,2:11,,]-output1)) 
[1] 1.110223e-16 
> .Machine$double.eps 
[1] 2.220446e-16 

Обратите внимание, что output1 не содержит output[,1,,] , Но эта матрица просто заполнена единицами:

> all(output[,1,,]==1) 
[1] TRUE 

Таким образом output1 может быть легко расширен, если это желательно.

Для nx = ny = nz = nt = 10 этот метод явно лучше:

nx = ny = nz = nt = 10 
A = array(runif(nx*ny*nz*nt), c(nx,ny,nz,nt)) 

f.old <- function(){ 
    output = array(1, c(nx, ny+1, nz, nt)) 
    for (x in 1:nx) 
    { 
    for (z in 1:nz) 
    { 
     for (t in 1:nt) 
     { 
     for (y in 2:(ny+1)) 
     { 
      output[x,y,z,t] = output[x,y-1,z,t] * (1 - A[x,y-1,z,t]) 
     } 
     } 
    } 
    } 
} 

f.new <- function() aperm(apply(A,c(1,3,4),function(v) cumprod(1-v)),c(2,1,3,4)) 

microbenchmark Тогда доходность (на моей машине):

> microbenchmark(f.old(),f.new()) 
Unit: milliseconds 
    expr  min  lq  mean median  uq  max neval 
f.old() 49.553825 53.486576 61.701149 57.710147 62.862921 136.02883 100 
f.new() 2.036781 2.365426 2.988266 2.685126 3.396083 10.88668 100