apply関数は遅かった・・・

うーん、こりゃ参ったねｗ

n <- 20
rep <- 100000
d <- matrix(rnorm(n*rep), nrow=n, ncol=rep )

# apply
Rprof()
varA <- apply(d, 2, var)
Rprof(NULL)
print(a <- summaryRprof()$sampling.time)

# matrix operation (1)
Rprof()
varB <- NULL
for(i in 1:1000) {
varB <- cbind(
  varB, 
  (diag(t(d[,((i-1)*100+1):((i)*100)])
   %*% d[,((i-1)*100+1):((i)*100)])
   - (rep(1, n) %*% d[,((i-1)*100+1):((i)*100)])^2 / n)
   / (n - 1)
)}
Rprof(NULL)
print(b <- summaryRprof()$sampling.time)

# matrix operation (2)
Rprof()
varC <- NULL
for(i in 1:1000) {
varC <- cbind(
  varC, 
  (rep(1, n) %*% d[,((i-1)*100+1):((i)*100)]^2 -
  (rep(1, n) %*% d[,((i-1)*100+1):((i)*100)])^2 / n)
  / (n - 1)
)}
Rprof(NULL)
print(c <- summaryRprof()$sampling.time)

max(varA-varB)
max(varB-varC)
b/a
c/b
c/a

> # apply
> Rprof()
> varA <- apply(d, 2, var)
> Rprof(NULL)
> print(a <- summaryRprof()$sampling.time)
[1] 5.02

> # matrix operation (1)
> Rprof()
> varB <- NULL
> for(i in 1:1000) {
+ varB <- cbind(
+   varB, 
+   (diag(t(d[,((i-1)*100+1):((i)*100)])
+    %*% d[,((i-1)*100+1):((i)*100)])
+    - (rep(1, n) %*% d[,((i-1)*100+1):((i)*100)])^2 / n)
+    / (n - 1)
+ )}
> Rprof(NULL)
> print(b <- summaryRprof()$sampling.time)
[1] 0.68

> # matrix operation (2)
> Rprof()
> varC <- NULL
> for(i in 1:1000) {
+ varC <- cbind(
+   varC, 
+   (rep(1, n) %*% d[,((i-1)*100+1):((i)*100)]^2 -
+   (rep(1, n) %*% d[,((i-1)*100+1):((i)*100)])^2 / n)
+   / (n - 1)
+ )}
> Rprof(NULL)
> print(c <- summaryRprof()$sampling.time)
[1] 0.34

> max(varA-varB)
[1] 1.332268e-15
> max(varA-varC)
[1] 1.332268e-15
> b/a
[1] 0.1354582
> c/b
[1] 0.5
> c/a
[1] 0.06772908

applyなら大丈夫だろうと思ってたら、今まで色々と損してました。
シミュレーション組むときは行列演算を駆使して書こう。

~~diag(t(A) %*% A)で損していると思われるのに、applyよりも全然早いです。(1)~~
メモリがオーバーフローしないように、念のためforループで分割しています。

高速 tapply(X, INDEX, var) もどき関数

Triad sou.

apply関数は遅かった・・・

関連記事