Ускорьте расчет Tidyverse нескольких квантилей

Question

У меня есть эта замечательная маленькая функция summarise_posterior (приведенная ниже) как часть моего пакета driver ( доступна здесь ).

Функция отличная и супер полезная. Одна проблема заключается в том, что я работаю с большими и большими данными, и это может быть очень медленным. Короче говоря, мой вопрос таков: есть ли способ ускорения этого процесса, сохраняя при этом ключевую гибкость этой функции (см. Примеры в документации).

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

#' Shortcut for summarize variable with quantiles and mean
#'
#' @param data tidy data frame
#' @param var variable name (unquoted) to be summarised
#' @param ... other expressions to pass to summarise
#'
#' @return data.frame
#' @export
#' @details Notation: \code{pX} refers to the \code{X}\% quantile
#' @import dplyr
#' @importFrom stats quantile
#' @importFrom rlang quos quo UQ
#' @examples
#' d <- data.frame("a"=sample(1:10, 50, TRUE),
#'                 "b"=rnorm(50))
#'
#' # Summarize posterior for b over grouping of a and also calcuate
#' # minmum of b (in addition to normal statistics returned)
#' d <- dplyr::group_by(d, a)
#' summarise_posterior(d, b, mean.b = mean(b), min=min(b))
summarise_posterior <- function(data, var, ...){
  qvar <- enquo(var)
  qs <- quos(...)


  data %>%
    summarise(p2.5 = quantile(!!qvar, prob=0.025),
              p25 = quantile(!!qvar, prob=0.25),
              p50 = quantile(!!qvar, prob=0.5),
              mean = mean(!!qvar),
              p75 = quantile(!!qvar, prob=0.75),
              p97.5 = quantile(!!qvar, prob=0.975),
              !!!qs)
}

Опции Rcpp также приветствуются.

Спасибо!

Answer 1

Вот решение, которое использует вложение, чтобы избежать многократного вызова quantile. В любое время, когда вам нужно сохранить вектор результатов внутри summarize, просто оберните его внутри list. После этого вы можете удалить эти результаты, сопоставить их с именами и использовать spread, чтобы поместить их в отдельные столбцы:

summarise_posterior2 <- function(data, var, ...){
  qvar <- ensym(var)
  vq <- c(0.025, 0.25, 0.5, 0.75, 0.975)

  summarise( data, .qq = list(quantile(!!qvar, vq, names=FALSE)),
             .nms = list(str_c("p", vq*100)), mean = mean(!!qvar), ... ) %>%
  unnest %>% spread( .nms, .qq )  
}

Это не дает вам такой же скорости, как у решения @ jay.sf

d <- data.frame("a"=sample(1:10, 5e5, TRUE), "b"=rnorm(5e5))    
microbenchmark::microbenchmark( f1 = summarise_posterior(d, b, mean.b = mean(b), min=min(b)),
                                f2 = summarise_posterior2(d, b, mean.b = mean(b), min=min(b)) )
# Unit: milliseconds
#  expr      min       lq     mean   median       uq      max neval
#    f1 49.06697 50.81422 60.75100 52.43030 54.17242 200.2961   100
#    f2 29.05209 29.66022 32.32508 30.84492 32.56364 138.9579   100

, но он будет корректно работать с group_by и внутри вложенных функций (, тогда как решения на основе substitute сломаются при вложении) .

r1 <- d %>% dplyr::group_by(a) %>% summarise_posterior(b, mean.b = mean(b), min=min(b))
r2 <- d %>% dplyr::group_by(a) %>% summarise_posterior2(b, mean.b = mean(b), min=min(b))
all_equal( r1, r2 )     # TRUE

Если вы профилируете код, вы увидите, где основные зависания

Rprof()
for( i in 1:100 )
  d %>% dplyr::group_by(a) %>% summarise_posterior2(b, mean.b = mean(b), min=min(b))
Rprof(NULL)
summaryRprof()$by.self %>% head
#             self.time self.pct total.time total.pct
# ".Call"          1.84    49.73       3.18     85.95
# "sort.int"       0.94    25.41       1.12     30.27
# "eval"           0.08     2.16       3.64     98.38
# "tryCatch"       0.08     2.16       1.44     38.92
# "anyNA"          0.08     2.16       0.08      2.16
# "structure"      0.04     1.08       0.08      2.16

.Call соответствует в основном бэкэнду C ++ dplyr, тогда как sort.int является рабочим позади quantile(). Решение @ jay.sf значительно ускоряется за счет отделения от dplyr, но также теряет соответствующую гибкость (например, интеграция с group_by). В конечном счете, вам решать, что важнее.

Answer 2

Почему бы не что-то подобное?

summarise_posterior2 <- function(data, x, ...){
  x <- deparse(substitute(x))
  nm <- deparse(substitute(...))
  M <- matrix(unlist(data[, x]), ncol=length(data[, x]))
  qs <- t(sapply(list(...), do.call, list(M)))
  'rownames<-'(cbind(p2.5 = quantile(M, prob=0.025),
        p25 = quantile(M, prob=0.25),
        p50 = quantile(M, prob=0.5),
        mean = mean(M),
        p75 = quantile(M, prob=0.75),
        p97.5 = quantile(M, prob=0.975), qs), NULL
  )
}

> summarise_posterior2(df1, X4, mean=mean, mean=mean, min=min)
     p2.5 p25 p50 mean p75 p97.5 mean mean min
[1,] 28.2  30  32   32  34  35.8   32   32  28

# > summarise_posterior(df1, X4, mean.b = mean(X4), min=min(X4))
#   p2.5 p25 p50 mean p75 p97.5 mean.b min
# 1 28.2  30  32   32  34  35.8     32  28

Работает в шесть раз быстрее:

> microbenchmark::microbenchmark(orig.fun=summarise_posterior(df1, X4, max(X4), min(X4)),
+                                new.fun=summarise_posterior2(df1, X4, max=max, min=min))
Unit: microseconds
     expr      min       lq      mean   median       uq      max neval
 orig.fun 4289.541 4324.490 4514.1634 4362.500 4411.225 8928.316   100
  new.fun  716.071  734.694  802.9949  755.867  778.317 4759.439   100

Данные

df1 <- data.frame(matrix(1:144, 9, 16))