COLOR(magenta){SIZE(20){Rコードの最適化例：トランプ}} [[Rコード最適化のコツと実例集]] に戻る~
The R Book 中の例を最適化する。~
トランプを良く切って、10枚取る。最初の5枚の絵柄か。残り5枚の絵柄が同じになる
確率をモンテカルロ法で求める。1:13がダイヤ, 14:26がクラブ, 27:39がハート, 40:52がスペードとする
## The R Book 中のオリジナルコード(141ページのプログラム)
orig <- function (rep1) { # repl は実験回数
(本を御覧下さい)
}
## 匿名子Aさんの例。判定条件を整理し、10枚取り出す操作を sample 関数で記述
test1 <- function (rep1) {
check2 <- function(x) { # 5枚の絵柄が同じかどうか判定する補助関数
any(all(x <= 13), all((14 <= x) & (x <= 26)),
all((27 <= x) & (x <= 39)), all(x >= 40))
}
pass <- 0; try <- 0
while(try < rep1) {
X <- sample(1:52, 10) # 52枚から無作為に10枚取る
if (check2(X[1:5]) | check2(X[6:10])) pass <- pass+1
try <- try+1
}
pass/try
}
## pass <- pass + (check2(X[1:5]) | check2(X[6:10])) とすれば(見かけ上) if 文追放可能
## test1 の論理判断部分を改善 -> 論理判断部分の最適化が重要という例
test1b <- function (rep1) {
aux <- function(x) {
x <- (x-1) %/% 13 # 絵柄情報 0,1,2,3 に変更
sum(x==x[1]) } # これが 5 なら 5 枚の絵柄が一致
pass <- 0; try <- 0
while(try < rep1) {
X <- sample(1:52, 10) # 52枚から無作為に10枚取る
if (aux(X[1:5])==5 | aux(X[6:10])==5) pass <- pass+1
try <- try+1
}
pass/try
}
## 抽出結果を一旦ベクトルに格納し、全体として操作。判定条件を一回の == 演算に置き換え
test2 <- function (repl) {
x <- numeric(10*repl) # repl 回の抽出結果を一括して納めるベクトルを用意
# 後の剰余演算のために 1:52 ではなく 0:51 としておく
for (i in 1:repl) x[(10*i-9):(10*i)] <- sample(0:51, 10) # 抽出結果を全部並べる
# 5枚ずつ行列に並べ、絵柄情報(0,1,2,3)だけにする
y <- (matrix(x, nrow=5)) %/% 13
z <- y - matrix(y[1,], nrow=5, ncol=2*repl, byrow=TRUE)
# 各5枚の絵柄から、最初の絵柄を引く(これが 0,0,0,0,0 なら絵柄一致!)
w <- matrix(colSums(abs(z)), nrow=2) # 各列毎の絶対値総和を求め、それを列数2の行列に直す
# もし w[1,i]*w[2,i]=0 ならi番目の実験成功
pass <- sum(w[1,]*w[2,]==0) # 成功した回数を数える(ここで始めて論理判断)
pass/repl
}
## 結局絵柄だけが問題と喝破し、絵柄 0,1,2,3 の抽出に置き換え
test3 <- function (repl) {
x <- numeric(10*repl) # repl 回の抽出結果を一括して納めるベクトルを用意
a <- rep(0:3,13) # 絵柄 0,1,2,3 が13通り
for (i in 1:repl) x[(10*i-9):(10*i)] <- sample(a, 10) # 抽出結果を全部並べる
y <- matrix(x, nrow=5)
z <- y - matrix(y[1,], nrow=5, ncol=2*repl, byrow=TRUE)
w <- matrix(colSums(abs(z)), nrow=2)
pass <- sum(w[1,]*w[2,]==0)
pass/repl
}
## ベクトルに代入する代わりに、行列に代入 (これだけでも違う)
test4 <- function (repl) {
x <- matrix(0, nrow = 10, ncol=repl) # 全ての抽出結果を納める 10 x repl 行列
a <- rep(0:3,13)
for (i in 1:repl) x[,i] <- sample(a, 10) # 毎回の結果を一列に納める
y <- matrix(x, nrow=5)
z <- y - matrix(y[1,], nrow=5, ncol=2*repl, byrow=TRUE)
w <- matrix(colSums(abs(z)), nrow=2)
pass <- sum(w[1,]*w[2,]==0)
pass/repl
}
## 一度に10枚だけ取り出すのはもったいないので,50枚取り出す
## 引き続く5回の実験は同分布だが独立にならないことを注意 (MCMC法では独立でない標本を使う方が普通 <- 言い訳)
## 標本抽出の繰り返しがボトルネックであることが良く分かる
test5 <- function (repl) { # repl が 5 の倍数であることが前提
x <- matrix(0, nrow = 50, ncol=repl/5) # 全ての抽出結果を納める 50 x (repl/5) 行列
a <- rep(0:3,13)
for (i in 1:(repl/5)) x[,i] <- sample(a, 50) # 50枚の結果を一列に納める
y <- matrix(x, nrow=5)
z <- y - matrix(y[1,], nrow=5, ncol=2*repl, byrow=TRUE)
w <- matrix(colSums(abs(z)), nrow=2)
pass <- sum(w[1,]*w[2,]==0)
pass/repl
}
## 究極のコード。インチキなどと言わないで下さい。そもそもこれを計算したかったのでしょう？
## R ではなく頭を使うことも大事と言う教訓のつもりです。
test6 <- function () {2*4*choose(13, 5)/choose(52, 5)}
## 一万回の実験を10回の繰り返したときの平均時間計測(使用環境で異なる可能性)
37.17 0.00 74.33 0.00 0.00 # orig オリジナル版との速度比
13.85 0.00 27.76 0.00 0.00 # test1 2.7 倍
0.794 0.000 1.585 0.00 0.00 # test1b 46.8 倍
0.556 0.007 1.118 0.00 0.00 # test2 66.9 倍
0.548 0.006 1.099 0.00 0.00 # test3 67.8 倍
0.453 0.013 0.929 0.00 0.00 # test4 82.1 倍
0.176 0.011 0.366 0.00 0.00 # test5 211.2 倍
0 0 0 0 0 # test6 ∞ 倍 (実際は 970496 倍)
## 念のためチェック
c(test1(10000),test1b(10000),test2(10000),test3(10000),test4(10000),test5(10000))
[1] 0.0035 0.0042 0.0027 0.0035 0.0043 0.0039
## コードの簡潔さ＋素直さを含めた最適化という点ではこんなのどうでしょう
# 頻繁に使う定数は変数 a にいれておく
# aux 関数は一回の実験とその結果の判定関数。判定パスなら TRUE、さもなければ FALSE を返す
# aux 関数の引数 i は sapply 関数中で使うために必要なダミー変数
# 最後に sapply 関数でループを隠す
a <- rep(0:3, 13)
aux <- function(i) {x <- sample(a,10); sum(abs(x[2:5]-x[1]))==0 | sum(abs(x[7:10]-x[6]))==0}
sum(sapply(1:repl, aux))/repl # repl 回のシミュレーション
[1] 0.0037
[1] 0.718 0.000 1.431 0.000 0.000 # それなりに早い
![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
