Re: Rで「ガチャとは心の所作」

id:bob3 さんから Rで「ガチャとは心の所作」 #TokyoSciPy #rstatsj - 極めて個人的なメモ

二重ループで効率悪いのですが、きっと id:a_bicky さんや裏RjpWikiの中の人がもっと良い方法を示唆してくれるものと思います。

とかプレッシャーを与えられたのでこんな感じかなぁと書いてみました。

ベースライン

元のコードを自分好みにちょっと変更した次の関数をベースラインとします。

gacha.sim <- function(prob, nop = 10000, iter.max = 10000) {
    nop.result <- numeric(nop)
    items <- seq(along = prob)
    result <- matrix(sample(items, nop * iter.max, replace = TRUE, prob = prob), iter.max)
    for (i in seq(length = nop)) {
        for (j in seq(length = iter.max)) {
            if (setequal(unique(result[1:j, i]), items)) {
                nop.result[i] <- j
                break
            }
        }
    }
    list(result = nop.result, prob = prob)
}

※アイテムコンプできなくても終了しますが・・・

ピュアRで個人的にベストだと思う関数

ベースラインの大きなボトルネックは2つで

  • matrix の生成
  • 終了判定(アイテムコンプ)のための for ルーブ
    です。

というわけで、これらを解消したのが次の関数です。

gacha.sim2 <- function(prob, nop = 10000, iter.max = 10000) {
    nop.result <- numeric(nop)
    items <- seq(along = prob)
    for (i in seq(length = nop)) {
        nop.result[i] <- max(match(items, sample.int(length(items), iter.max, TRUE, prob)))
    }
    list(result = nop.result, prob = prob)
}

matrix の生成に関しては、一度に全ての値を生成するのではなく1人分の値しか生成しないことと、sample.int を使うことがポイントです。
ベースラインの書き方で sample を使うと最終的に次のコードが実行され、.Internal の結果を直接返さない分無駄が生じます。

        x[.Internal(sample(length(x), size, replace, prob))]

今回の場合、sample に渡すのは文字列ベクトルではなく整数のベクトルなので、sample.int、あるいは sample の第1引数に整数ベクトルの長さ(長さ1の数値ベクトル)を使うことで .Internal の結果を直接得ることができます。

終了判定の for ループは、繰り返し回数が増えれば増えるほど速度面で致命的になるので、それを解消するために match を使っています。

ベンチマーク

次のような感じになりました。確率に偏りがある場合にかなり高速化されているのがわかるかと思います。

> library(rbenchmark)
> p <- rep(1, 6); benchmark(gacha.sim(p), gacha.sim2(p), replications = 1)
           test replications elapsed relative user.self sys.self user.child sys.child
1  gacha.sim(p)            1  10.586 1.618902     9.537    1.044          0         0
2 gacha.sim2(p)            1   6.539 1.000000     4.907    1.631          0         0
> p <- c(100, 50, 10, 10, 3, 1); benchmark(gacha.sim(p), gacha.sim2(p), replications = 1)
           test replications elapsed relative user.self sys.self user.child sys.child
1  gacha.sim(p)            1  81.277 12.81972     79.26    1.939          0         0
2 gacha.sim2(p)            1   6.340  1.00000      4.56    1.772          0         0

bob3 さんと同じ関数を使っても2倍ぐらい時間がかかるんですが、どんなハイスペックマシン使ってるんでしょうかね・・・

おまけ

次のような書き方もできますが、アイテムの取得確率に偏りがあると繰り返し回数が多くなって遅くなるので inline とか使わないと効果はなさそうです。

gacha.sim3 <- function(prob, nop = 10000, iter.max = 10000) {
    nop.result <- numeric(nop)
    items <- seq(along = prob)
    for (i in seq(length = nop)) {
        my.items <- c()
        for (j in seq(length = iter.max)) {
            my.items[sample.int(length(items), 1, prob = prob)] <- 1
            if (sum(my.items, na.rm = TRUE) == length(items)) {
                nop.result[i] <- j
                break
            }
        }
    }
    list(result = nop.result, prob = prob)
}

追記(2012/4/16)

R 2.14 から R-core になった compiler を試してみました。あと、バージョンをR 2.14.1 から R 2.15.0 に上げました。
他の環境は同じです。

library(compiler)
cgacha.sim <- cmpfun(gacha.sim)
cgacha.sim2 <- cmpfun(gacha.sim2)
cgacha.sim3 <- cmpfun(gacha.sim3)

以下、ベンチマークの結果です。

> library(rbenchmark)
> # 等確率
> p <- rep(1, 6); benchmark(gacha.sim(p), gacha.sim2(p), gacha.sim3(p), replications = 1)
           test replications elapsed relative user.self sys.self user.child sys.child
1  gacha.sim(p)            1  10.301 4.893587     9.326    0.971          0         0
2 gacha.sim2(p)            1   6.713 3.189074     4.889    1.822          0         0
3 gacha.sim3(p)            1   2.105 1.000000     1.934    0.171          0         0
> # 等確率(コンパイル版)
> p <- rep(1, 6); benchmark(cgacha.sim(p), cgacha.sim2(p), cgacha.sim3(p), replications = 1)
            test replications elapsed relative user.self sys.self user.child sys.child
1  cgacha.sim(p)            1   9.844 6.696599     8.918    0.923          0         0
2 cgacha.sim2(p)            1   6.775 4.608844     4.830    1.944          0         0
3 cgacha.sim3(p)            1   1.470 1.000000     1.293    0.176          0         0
> # 確率に偏りがある場合
> p <- c(100, 50, 10, 10, 3, 1); benchmark(gacha.sim(p), gacha.sim2(p), gacha.sim3(p), replications = 1)
           test replications elapsed relative user.self sys.self user.child sys.child
1  gacha.sim(p)            1  54.373 8.506414    51.633    2.719          0         0
2 gacha.sim2(p)            1   6.392 1.000000     4.384    2.005          0         0
3 gacha.sim3(p)            1  20.533 3.212297    19.686    0.838          0         0
> # 確率に偏りがある場合(コンパイル版)
> p <- c(100, 50, 10, 10, 3, 1); benchmark(cgacha.sim(p), cgacha.sim2(p), cgacha.sim3(p), replications = 1)
            test replications elapsed  relative user.self sys.self user.child sys.child
1  cgacha.sim(p)            1  72.244 10.871934    70.234    1.987          0         0
2 cgacha.sim2(p)            1   6.645  1.000000     4.354    2.289          0         0
3 cgacha.sim3(p)            1  12.826  1.930173    11.780    1.042          0         0

gacha.sim と gacha.sim3 はコンパイルするとかなり速くなると思ったんですが、それほど速くなりませんでしたね・・・。
1回の実行しか見てないんで厳密さには欠けますが、コンパイルしたところで gacha.sim や gacha.sim3 が gacha.sim2 より速くなることはまずなさそうですし、compiler に頼るよりも自分で最適化した方が速くなりそうですね。

広告
一部のデータフレームにおいてfor ループがcolMeansより速い理由 潜在的意味インデキシング(LSI)徹底入門
※このエントリーははてなダイアリーから移行したものです。過去のコメントなどはそちらを参照してください