Какой самый быстрый способ применить t.test к каждому столбцу большой матрицы?

Предположим, у меня есть большая матрица:

M <- matrix(rnorm(1e7),nrow=20)

Далее предположим, что каждый столбец представляет выборку. Скажем, я хотел бы применить t.test() к каждому столбцу, есть ли способ сделать это намного быстрее, чем использование apply()?

apply(M, 2, t.test)

Запуск анализа на моем компьютере занял чуть меньше 2 минут:

> system.time(invisible( apply(M, 2, t.test)))
user  system elapsed 
113.513   0.663 113.519

r apply matrix statistics hypothesis-test
person Alex    schedule 12.07.2012    source источник
comment
apply очень гибкая функция и поэтому включает в себя множество вещей, которые вам не нужны в каждом конкретном случае. Вероятно, кодирование той же логики вручную с помощью цикла for даст некоторое увеличение производительности.   -  person ffriend    schedule 13.07.2012

Ответы (2)


arrow_upward
8
arrow_downward

Если у вас многоядерная машина, использование всех ядер дает некоторые преимущества, например, использование mclapply.

> library(multicore)
> M <- matrix(rnorm(40),nrow=20)
> x1 <- apply(M, 2, t.test)
> x2 <- mclapply(1:dim(M)[2], function(i) t.test(M[,i]))
> all.equal(x1, x2)
[1] "Component 1: Component 9: 1 string mismatch" "Component 2: Component 9: 1 string mismatch"
# str(x1) and str(x2) show that the difference is immaterial

Этот мини-пример показывает, что все идет по плану. Теперь увеличьте:

> M <- matrix(rnorm(1e7), nrow=20)
> system.time(invisible(apply(M, 2, t.test)))
   user  system elapsed 
101.346   0.626 101.859
> system.time(invisible(mclapply(1:dim(M)[2], function(i) t.test(M[,i]))))
  user  system elapsed 
55.049   2.527  43.668

Это использует 8 виртуальных ядер. Ваш пробег может отличаться. Небольшой выигрыш, но он достигается при очень небольших усилиях.

ИЗМЕНИТЬ

Если вас интересует только сама t-статистика, извлечение соответствующего поля ($statistic) немного ускорит работу, особенно в многоядерном случае:

> system.time(invisible(apply(M, 2, function(c) t.test(c)$statistic)))
   user  system elapsed 
 80.920   0.437  82.109 
> system.time(invisible(mclapply(1:dim(M)[2], function(i) t.test(M[,i])$statistic)))
   user  system elapsed 
 21.246   1.367  24.107

Или, что еще быстрее, вычислите значение t напрямую.

my.t.test <- function(c){
  n <- sqrt(length(c))
  mean(c)*n/sd(c)
}

потом

> system.time(invisible(apply(M, 2, function(c) my.t.test(c))))
   user  system elapsed 
 21.371   0.247  21.532 
> system.time(invisible(mclapply(1:dim(M)[2], function(i) my.t.test(M[,i]))))
   user  system elapsed 
144.161   8.658   6.313
person Ryogi    schedule 12.07.2012
comment
Я думаю, что просто вычислю статистику t напрямую, что, как вы показали, намного быстрее. - person Alex; 13.07.2012

arrow_upward
9
arrow_downward

Вы можете добиться большего успеха с помощью функции colttests из пакета genefilter (на Bioconductor).

> library(genefilter)
> M <- matrix(rnorm(40),nrow=20)
> my.t.test <- function(c){
+   n <- sqrt(length(c))
+   mean(c)*n/sd(c)
+ }
> x1 <- apply(M, 2, function(c) my.t.test(c))
> x2 <- colttests(M, gl(1, nrow(M)))[,"statistic"]
> all.equal(x1, x2)
[1] TRUE
> M <- matrix(rnorm(1e7), nrow=20)
> system.time(invisible(apply(M, 2, function(c) my.t.test(c))))
   user  system elapsed 
 27.386   0.004  27.445 
> system.time(invisible(colttests(M, gl(1, nrow(M)))[,"statistic"]))
   user  system elapsed 
  0.412   0.000   0.414

Ссылка: «Одновременное вычисление тысяч тестовых статистических данных в R», SCGN, том 18 (1), 2007 г., http://stat-computing.org/newsletter/issues/scgn-18-1.pdf.

person Heather Turner    schedule 13.07.2012