Effectuer efficacement un test de distribution par rangée

Question

J'ai une matrice dans laquelle chaque ligne est un échantillon d'une distribution. Je veux faire une comparaison de roulement des distributions en utilisant ks.test et enregistrer la statistique de test dans chaque cas. La façon la plus simple à mettre en œuvre ce plan conceptuel est une boucle:

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 

} 

Cependant, mon vrai données a ~ 400 colonnes et rangées ~ 300.000 pour un seul exemple, et j'ai beaucoup d'exemples. Donc j'aimerais que ça soit rapide. Le test de Kolmogorov-Smirnov n'est pas compliqué du point de vue mathématique, donc si la réponse est "implémentez-le dans Rcpp", je l'accepterai à contrecoeur, mais je serais plutôt surpris - il est déjà très rapide de calculer sur un seul paire R.

Méthodes J'ai essayé, mais ont été incapables de faire fonctionner correctement: dplyr en utilisant rowwise/do/lag, en utilisant zoorollapply (qui est ce que j'utilise pour générer les distributions), et alimenter une data.table dans une boucle (modifier: celui-ci fonctionne, mais c'est encore lent).

Answer 1

Une mise en œuvre rapide et sale dans CRPP

// [[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; 
} 

pour la matrice de la taille (400, 30000), il termine en 1 s.

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

Et le résultat semble être précis.

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

Une des sources d'accélération est d'écrire une version plus petite de ks.test qui en fait moins. ks.test2 ci-dessous est plus restrictif que ks.test. Cela suppose, par exemple, que vous n'avez pas de valeurs manquantes et que vous voulez toujours que la statistique soit associée à un test à deux faces.

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

} 

Vérifiez que la sortie est conforme à 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 

déterminer maintenant les économies de la fonction plus petite:

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

est ici une solution dplyr qui obtient le même résultat que votre boucle. J'ai des doutes si c'est en fait plus rapide que la boucle, mais peut-être cela peut-il servir de premier pas vers une solution.

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

j'ai pu calculer la statistique Kruskal-Wallis avec ks.test() en utilisant par paires rollapplyr().

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

Ceci obtient le résultat attendu, mais il est lent pour un jeu de données de votre taille. Lent lentement lent. C'est peut-être parce que ks.test() calcule beaucoup plus que la statistique à chaque itération; il obtient également la valeur p et fait beaucoup de vérification d'erreur.

En effet, si nous simulons un grand ensemble de données comme ceci:

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

La solution rollapplyr() prend beaucoup de temps; J'ai arrêté l'exécution après environ 2 heures, à quel point il avait calculé presque tous les résultats (mais pas tous).

Il semble que si rollapplyr() est probablement plus rapide qu'une boucle for, il ne sera probablement pas la meilleure solution globale en termes de performances.