Произведение столбцов по строкам с динамическим вводом столбцов — операция векторизации

#r

#r

Вопрос:

Я хотел бы векторизовать следующий код для более эффективной обработки. Мне нужно взять произведение столбцов по строкам (т.Е. rowProds), но количество столбцов, из которых я хотел бы получить произведение, должно быть функцией другого ввода.

Если возможно, я бы предпочел, чтобы это было сделано с использованием Base R, но я открыт и ценю любые предложения.

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

 # Generate some data

mat <- data.frame(X = 1:5)
for (i in 1:5) {
  set.seed(i)
  mat[1   i] <- runif(5)
}

# Via a for loop

for (i in 1:nrow(mat)) {  
  mat$calc[i] <- prod(mat[match(mat$X[i], mat$X), 2:(i   1)])
}
mat

# Via a function with mapply

rowprodfun <- function(X) {  
  myprod <- prod(mat[match(X, mat$X), 2:(X   1)])
  return(myprod)
}

mat$calc <- mapply(rowprodfun, mat$X)
mat

mat$calc
# [1] 0.265508663 0.261370165 0.126427355 0.013874517 0.009758232
  

Оба описанных выше метода приводят к одному и тому же столбцу «calc». Мне просто нужен более эффективный способ создания этого столбца.

Ответ №1:

Одним из вариантов было бы преобразовать элементы верхнего треугольника как NA , а затем использовать rowProds из matrixStats

 library(matrixStats)
rowProds(as.matrix(mat[-1] * NA^upper.tri(mat[-1])), na.rm = TRUE)
#[1] 0.265508663 0.261370165 0.126427355 0.013874517 0.009758232
  

Комментарии:

1.Я согласен, что это еще одно возможное решение. Однако я не уверен, что достигается эффективность. library(microbenchmark)' 'library(matrixStats) Y <- microbenchmark( for (i in 1:nrow(mat)) {prod(mat[match(mat$X[i], mat$X), 2:(i 1)])}, mapply(rowprodfun, mat$X), rowProds(as.matrix(mat[-1] * NA ^ upper.tri(mat[-1])), na.rm = T), apply(mat[-1] * NA ^ upper.tri(mat[-1]), 1, prod, na.rm = T), # Not solutions, just benchmark reference points apply(mat[, 2:ncol(mat)], 1, prod), rowSums(mat[, 2:ncol(mat)]) )

2. @rwteg Извините, я не проверял тест. Думал, что rowProds это должно быть быстрее, могут быть преобразования upper.tri etc, замедляющие процесс. Для эффективности по строкам может быть лучше переписать ваш for цикл в Rcpp

3. Похоже upper.tri , что это замедляет вычисление, поскольку, когда я просто использую rowProds(as.matrix(mat[-1])) скорость, она соответствует скорости rowSums. Я могу поиграть с Rcpp . Спасибо за предложение. PS извините за форматирование, я все еще разбираюсь в этом!

4. @rwteg Если это так, есть одно небольшое изменение

5. На самом деле это похоже на умножение data.frames, которое замедляет это. Есть мысли о лучшем решении?

Ответ №2:

Использование upper.tri , предложенное @akrun, было очень полезным. Заключительной частью было преобразование data.frame в матрицу as.matrix перед выполнением поэлементного умножения.

rowProds(as.matrix(mat[-1]) * NA ^ upper.tri(mat[-1]), na.rm = T) результат наиболее эффективного вычисления.

apply(as.matrix(mat[-1]) * NA ^ upper.tri(mat[-1]), 1, prod, na.rm = T) было почти так же эффективно, если пытаться выполнить в базе R.

 library(microbenchmark)
library(matrixStats)
library(ggplot2)

Y <- microbenchmark(
  for.loop = for (i in 1:nrow(mat)) {prod(mat[match(mat$X[i], mat$X), 2:(i   1)])},
  mapply.fun = mapply(rowprodfun, mat$X),
  rowProds = rowProds(as.matrix(mat[-1] * NA ^ upper.tri(mat[-1])), na.rm = T),
  rowProds.matrix = rowProds(as.matrix(mat[-1]) * NA ^ upper.tri(mat[-1]), na.rm = T),
  apply = apply(mat[-1] * NA ^ upper.tri(mat[-1]), 1, prod, na.rm = T),
  apply.matrix = apply(as.matrix(mat[-1]) * NA ^ upper.tri(mat[-1]), 1, prod, na.rm = T)
)

> Y
Unit: microseconds
            expr      min        lq      mean    median        uq       max neval
        for.loop 4094.869 4305.5590 5682.2124 4479.8125 5193.8190 50361.025   100
      mapply.fun  542.962  577.6995 1036.9821  599.2220  658.1245 32426.296   100
        rowProds  518.419  553.9120  654.2657  597.5225  637.1690  2434.267   100
 rowProds.matrix   99.304  116.1065  144.9313  128.0010  153.8650   516.909   100
           apply  547.493  580.1540  686.2317  628.2955  703.0565  1215.812   100
    apply.matrix  117.051  136.6845  158.3808  144.9920  156.5075   339.068   100

  

benchmarkautoplot