Kiểm tra phân phối hàng khôn ngoan một cách hiệu quả

Question

Tôi có ma trận trong đó mỗi hàng là mẫu từ bản phân phối. Tôi muốn so sánh các bản phân phối bằng cách sử dụng ks.test và lưu thống kê thử nghiệm trong mỗi trường hợp. Cách đơn giản nhất để thực hiện điều này khái niệm là với một vòng lặp:

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 

} 

Tuy nhiên, dữ liệu thực tế của tôi có ~ 400 cột và ~ 300.000 hàng cho một ví dụ duy nhất, và tôi có rất nhiều ví dụ. Vì vậy, tôi muốn điều này được nhanh chóng. Bài kiểm tra Kolmogorov-Smirnov không phức tạp về toán học, vì vậy nếu câu trả lời là "thực hiện nó trong Rcpp", tôi sẽ chấp nhận một cách miễn cưỡng, nhưng tôi sẽ hơi ngạc nhiên - nó đã rất nhanh để tính toán trên một đĩa đơn cặp trong R.

Phương pháp tôi đã cố gắng nhưng đã không thể có được làm việc: dplyr sử dụng rowwise/do/lag, zoo sử dụng rollapply (đó là những gì tôi sử dụng để tạo ra các bản phân phối), và populating một data.table trong một vòng lặp (chỉnh sửa: cái này hoạt động, nhưng nó vẫn còn chậm).

Answer 1

Một thực hiện nhanh chóng và dơ bẩn trong 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; 
} 

cho ma trận kích thước (400, 30000), nó hoàn thành dưới 1s.

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

Và kết quả có vẻ chính xác.

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

Một nguồn tăng tốc là viết phiên bản nhỏ hơn là ks.test. ks.test2 dưới đây hạn chế hơn ks.test. Ví dụ, giả sử rằng bạn không có giá trị bị thiếu và bạn luôn muốn thống kê được kết hợp với kiểm tra hai mặt.

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)) 

} 

Xác minh rằng đầu ra phù hợp với 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 

Bây giờ xác định các khoản tiết kiệm từ các chức năng nhỏ:

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

Dưới đây là một giải pháp dplyr mà được kết quả tương tự như vòng lặp của bạn. Tôi có nghi ngờ của tôi nếu điều này thực sự nhanh hơn vòng lặp, nhưng có lẽ nó có thể phục vụ như là một bước đầu tiên hướng tới một giải pháp.

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

tôi đã có thể tính toán các số liệu thống kê cặp Kruskal-Wallis sử dụng ks.test() với rollapplyr().

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

Điều này đạt được kết quả mong đợi nhưng chậm đối với tập dữ liệu có kích thước của bạn. Chậm chậm chạp. Điều này có thể là do ks.test() tính toán nhiều hơn chỉ số liệu thống kê tại mỗi lần lặp; nó cũng nhận được giá trị p và thực hiện rất nhiều kiểm tra lỗi.

Thật vậy, nếu chúng ta mô phỏng một bộ dữ liệu lớn như vậy:

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

Giải pháp rollapplyr() mất một thời gian dài; Tôi ngừng thực hiện sau khoảng 2 giờ, tại thời điểm đó nó đã tính toán gần như tất cả (nhưng không phải tất cả) kết quả.

Dường như trong khi rollapplyr() có khả năng nhanh hơn vòng lặp for, nó sẽ không có khả năng là giải pháp tổng thể tốt nhất về hiệu suất.