ケィオスの時系列解析メモランダム

時系列解析,生体情報学,数学・物理などの解説です.

【Rの処理時間】計算時間短縮のための工夫 ― meanよりもsumが速い ―

プログラムを作って何か処理するとき,処理時間が短い方が良いに決まってます.ということで,今回はRの処理時間についてのお話です.

 処理時間を無駄に長くしないための一般論として,私が説明できることは,「forループは使わずにベクトルのまま計算した方が良い」,「ifをできるだけ避ける」くらいです.私は,Rの内部の仕組みについて知らないので,Rに特有の挙動には精通していません.

 今回は「標本平均の計算」を例として,実験的に処理時間を検証してみます.

サンプルデータの生成

 以下では,xにデータが入っているとします.xの具体例として,長さN,平均3,標準偏差1の正規乱数を用意するRスクリプトは以下です.

# 数値実験用の乱数生成
# データの長さ
N <- 1e8 # 10^8と同じ
# 平均3,標準偏差1の正規乱数
x <- rnorm(N,3,1)

3通りの方法で標本平均を計算

 比較する標本平均の計算方法は以下の3通りです.
1. コマンドmeanを使う

mean(x)

2. コマンドsumで和を計算して,xの長さで割る

sum(x)/length(x)

3. forループを使って,プログラムっぽく計算

 s <- 0
 n <- length(x)
 for(i in 1:n){
  s <- s+x[i]
 }
 s/n

計算時間を比較

 標本平均を計算する3通りの方法の処理時間を比べてみます.

 使ったRスクリプトは以下です.

# 数値実験用の乱数生成
# データの長さ
N <- 1e8
# 平均3,標準偏差1の正規乱数
x <- rnorm(N,3)
##########################
# 結果比較用の変数を準備
mu <- c() # 標本平均の計算結果
cal.time <- c() # 処理時間
##########################
# 3通りの方法で標本平均を計算
##########################
# 1. コマンドmeanを使用
cal.time[1] <- system.time(mu[1] <- mean(x))[[1]]
# 2. コマンドsumで和を計算して,xの長さで割る
cal.time[2] <- system.time(mu[2] <- sum(x)/length(x))[[1]]
# 3. forループを使って,プログラムっぽく計算 
# forループを使用して関数fを定義
f <- function(x){
 s <- 0
 n <- length(x)
 for(i in 1:n){
  s <- s+x[i]
 }
 s/n
}
cal.time[3] <- system.time(mu[3] <- f(x))[[1]]
########
# 結果の比較
data.frame(method=c("mean","sum/N","for loop"),time=cal.time,mu)

 このRスクリプトを実行すると,処理時間time [s]と標本平均の計算結果muの表が出力されます.私の環境では,以下のようになりました.

    method time       mu
1     mean 0.14 2.999914
2    sum/N 0.08 2.999914
3 for loop 2.11 2.999914

ここで,methodは,それぞれ,上で説明した3通りの方法に対応します.

 当然,標本平均の計算結果muは,どの方法も同じです.しかし,処理時間time [s]は違います.

 意外なのは,meanよりも,sumで計算した後に標本数で割った方が処理時間が短かいことだと思います.半分くらいの時間になるので,

 Rでは,meanよりも,sumの方が速い

ということです.

欠損値NAがあると,処理時間は長くなる

 現実のデータでは,欠損値があることが普通です.データに欠損値NAが含まれると処理時間は長くなります.

 xに欠損値を一つ入れて,上と同じ数値実験をしてみます.ここでは,meansumに,

na.rm=TRUE

のオプションを設定して,欠損値を除きます.

 使ったスクリプトは以下です.

# 数値実験用の乱数生成
# データの長さ
N <- 1e8
# 平均3,標準偏差1の正規乱数
x <- rnorm(N,3)
# NAを1つ挿入
x[1] <- NA
##########################
# 結果比較用の変数を準備
mu <- c() # 標本平均の計算結果
cal.time <- c() # 処理時間
##########################
# 3通りの方法で標本平均を計算
##########################
# 1. コマンドmeanを使用
cal.time[1] <- system.time(mu[1] <- mean(x,na.rm=TRUE))[[1]]
# 2. コマンドsumで和を計算して,欠損値を除いたxの長さで割る
cal.time[2] <- system.time(mu[2] <- sum(x,na.rm=TRUE)/(length(x)-sum(is.na(x))))[[1]]
# 3. forループを使って,プログラムっぽく計算 
# forループを使用して関数fを定義
f <- function(x){
 s <- 0
 n <- 0
 for(tmp in x){
  if(!is.na(tmp)){
   s <- s+tmp
   n <- n + 1
  }
 }
 s/n
}
cal.time[3] <- system.time(mu[3] <- f(x))[[1]]
########
# 結果の比較
data.frame(method=c("mean","sum/N","for loop"),time=cal.time,mu)

 このRスクリプトを実行すると,処理時間time [s]と標本平均の計算結果muの表が出力されます.私の環境では,以下のようになりました.

    method  time       mu
1     mean  0.66 2.999827
2    sum/N  0.33 2.999827
3 for loop 12.66 2.999827

ここで,methodは,それぞれ,上で説明した3通りの方法に対応します.

 どの計算方法も,欠損値を気にする必要がなかったケースと比べて計算時間が長くなっています.

 注目すべき点は,今回も,sumを使った標本平均の計算が,最も処理時間が短いということです.meanの半分です.

 おまけで注意しておきますが,否定を使った処理!is.na(x)は,is.na(x)よりも処理時間がかかります.つまり,欠損値がある場合は,

sum(x,na.rm=TRUE)/(sum(!is.na(x)))

ではなく,

sum(x,na.rm=TRUE)/(length(x)-sum(is.na(x)))

とした方が,速いです.

まとめ

 私は18~22歳の間,大学に行かずに,工事現場などで毎日働いていました.その間,学問的なことには全く興味を失っていたので,勉強はしませんでした.その結果,23歳にして中学生レベルの知識しかもたない自分に気がつきました.何の能力もない自分に気づいた後は,勉強して新しいことを知る喜びを感じられたので,誰かに押しつけられるのではなく,自分自身でいろんなことに興味をもって勉強できるようになりました.今となっては人生には無駄な時間はないと感じます.

 人生には無駄な時間はないといっても,Rの処理時間は短い方が良いということが今回のメッセージです.