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  # それなりに早い

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Google
WWW を検索 OKADAJP.ORG を検索
Last-modified: 2015-03-01 (日) 01:15:59 (1717d)