Эффективное выполнение критерия распределения по строке

Question

У меня есть матрица, в которой каждая строка является образцом из распределения. Я хочу выполнить скользящее сравнение распределений с использованием ks.test и сохранить статистику теста в каждом случае. Самый простой способ осуществить это концептуально это с петлей:

set.seed(1942) 
mt <- rbind(rnorm(5), rnorm(5), rnorm(5), rnorm(5)) 

results <- matrix(as.numeric(rep(NA, nrow(mt)))) 

for (i in 2 : nrow(mt)) { 

    results[i] <- ks.test(x = mt[i - 1, ], y = mt[i, ])$statistic 

} 

Однако мои реальные данные имеют ~ 400 столбцов и строки ~ 300,000 для одного примера, и у меня есть много примеров. Поэтому я бы хотел, чтобы это было быстро. Тест Колмогорова-Смирнова - это не все, что математически сложно, поэтому, если ответ «реализует его в Rcpp», я неохотно соглашусь с этим, но я был бы несколько удивлен - это уже очень быстро вычислить на одном пара в R.

методы Я пытался, но не смог получить работу: dplyr с помощью rowwise/do/lag, zoo с помощью rollapply (что я использую для создания дистрибутивов), а также заполнение data.table в цикле (редактирование: этот работает, но он все еще медленный).

Answer 1

Быстрая и грязная реализация в Rcpp

// [[Rcpp::depends(RcppArmadillo)]] 
#include <RcppArmadillo.h> 

double KS(arma::colvec x, arma::colvec y) { 
    int n = x.n_rows; 
    arma::colvec w = join_cols(x, y); 
    arma::uvec z = arma::sort_index(w); 
    w.fill(-1); w.elem(find(z <= n-1)).ones(); 
    return max(abs(cumsum(w)))/n; 
} 
// [[Rcpp::export]] 
Rcpp::NumericVector K_S(arma::mat mt) { 
    int n = mt.n_cols; 
    Rcpp::NumericVector results(n); 
    for (int i=1; i<n;i++) { 
    arma::colvec x=mt.col(i-1); 
    arma::colvec y=mt.col(i); 
    results[i] = KS(x, y); 
    } 
    return results; 
} 

для матрицы размеров (400, 30000), он завершает под 1s.

system.time(K_S(t(mt)))[3] 
#elapsed 
# 0.98 

И результат кажется точным.

set.seed(1942) 
mt <- matrix(rnorm(400*30000), nrow=30000) 
results <- rep(0, nrow(mt)) 
for (i in 2 : nrow(mt)) { 
    results[i] <- ks.test(x = mt[i - 1, ], y = mt[i, ])$statistic 
} 
result <- K_S(t(mt)) 
all.equal(result, results) 
#[1] TRUE 

Answer 2

Одним из источников ускорения является запись меньшей версии ks.test, которая делает меньше. ks.test2 ниже является более ограничительным, чем ks.test. Предполагается, например, что у вас нет отсутствующих значений и что вы всегда хотите, чтобы статистика была связана с двухсторонним тестом.

ks.test2 <- function(x, y){ 

    n.x <- length(x) 
    n.y <- length(y) 
    w <- c(x, y) 
    z <- cumsum(ifelse(order(w) <= n.x, 1/n.x, -1/n.y)) 

    max(abs(z)) 

} 

Убедитесь, что выход соответствует ks.test.

set.seed(999) 
x <- rnorm(400) 
y <- rnorm(400) 

ks.test(x, y)$statistic 

    D 
0.045 

ks.test2(x, y) 

[1] 0.045 

Теперь определить экономию от меньшей функций:

library(microbenchmark) 

microbenchmark(
    ks.test(x, y), 
    ks.test2(x, y) 
) 

Unit: microseconds 
      expr  min  lq  mean median  uq  max neval cld 
    ks.test(x, y) 1030.238 1070.303 1347.3296 1227.207 1313.8490 6338.918 100 b 
ks.test2(x, y) 709.719 730.048 832.9532 833.861 888.5305 1281.284 100 a 

Answer 3

Вот dplyr решения, которое получает тот же результат, как и ваш цикл. У меня есть сомнения, если это на самом деле быстрее, чем цикл, но, возможно, он может служить первым шагом к решению.

require(dplyr) 
mt %>% 
    as.data.frame %>% 
    mutate_each(funs(lag)) %>% 
    cbind(mt) %>% 
    slice(-1) %>% 
    rowwise %>% 
    do({ 
    x = unlist(.) 
    n <- length(x) 
    data.frame(ks = ks.test(head(x, n/2), tail(x, n/2))$statistic) 
    }) %>% 
    unlist %>% 
    c(NA, .) %>% 
    matrix 

Answer 4

я смог вычислить парный Крускала-Уоллиса статистики с использованием ks.test() с rollapplyr().

results <- rollapplyr(data = big, 
         width = 2, 
         FUN = function(x) ks.test(x[1, ], x[2, ])$statistic, 
         by.column = FALSE) 

Это получает ожидаемый результат, но медленный набор данных для вашего размера. Медленное медленное замедление. Это может быть связано с тем, что ks.test() вычисляет намного больше, чем просто статистику на каждой итерации; он также получает значение p и выполняет большую проверку ошибок.

В самом деле, если мы моделируем большой набор данных, как так:

big <- NULL 
for (i in 1:400) { 
    big <- cbind(big, rnorm(300000)) 
} 

rollapplyr() решения занимает много времени; Я прекратил выполнение примерно через 2 часа, и в этот момент он вычислил почти все (но не все) результаты.

Похоже, что хотя rollapplyr(), скорее всего, быстрее, чем цикл for, это вряд ли будет лучшим общим решением с точки зрения производительности.