Hatena::ブログ(Diary)

極めて個人的なメモ Twitter

Θ・)ノシ Bob#3のメモ帳です。
1972 | 12 |
2003 | 03 | 04 | 05 | 06 | 11 | 12 |
2004 | 01 | 02 | 03 | 04 | 05 | 06 | 09 | 10 | 11 | 12 |
2005 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 09 | 10 | 11 |
2006 | 03 | 04 | 05 | 06 | 07 | 09 | 10 | 11 | 12 |
2007 | 01 | 02 | 03 | 05 | 06 | 07 | 10 | 11 | 12 |
2008 | 01 | 02 | 04 | 05 | 08 | 09 | 10 | 11 |
2009 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 |
2010 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 11 |
2011 | 05 | 06 | 09 | 10 |
2012 | 02 | 03 | 05 | 06 | 08 |
2013 | 01 | 02 | 03 |
2015 | 02 |

2012-03-18

[]Rで「ガチャとは心の所作」 #TokyoSciPy #rstatsj Rで「ガチャとは心の所作」 #TokyoSciPy #rstatsjを含むブックマーク Rで「ガチャとは心の所作」 #TokyoSciPy #rstatsjのブックマークコメント

Tokyo.SciPy#3id:AntiBayesian さんがLTされていた「ガチャとは心の所作」が非常に興味深かったのでRでもやってみました。

(もっとも、ustがなかったので私は発表資料を読んだだけですが…)


「ガチャとは心の所作」 d:id:AntiBayesian:20120318


Rで書いたらこんな感じになりました。

f:id:bob3:20120318173714p:image:w640


ガチャのシミュレーション関数はこんな感じ。

# probが各アイテムの出現確率
# nopがガチャを回す人数(既定値は10000人)
# iter.maxが一人が回す回数の上限(既定値は10000回)
gacha.sim <- function(prob, nop=10000, iter.max=10000){
  nop.result <- numeric(nop)
  items <- length(prob)
  result <- matrix(sample(1:items, nop*iter.max, replace=TRUE, prob=prob),
                   iter.max)
  for (i in 1:nop) {
    for (j in 1:iter.max) {
      if (length(unique(result[1:j, i]))==items){
        nop.result[i] <- j
        break
        }
      }
    }
  list("result" = nop.result, "prob"=prob)
  }

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


シミュレーションの実行はこんな感じで。

# ここでは4つの出現確率パターンで試してみる。
result1 <- gacha.sim(c(  1,  1,  1,  1, 1, 1))
result2 <- gacha.sim(c(100, 50, 10, 10, 3, 1))
result3 <- gacha.sim(c(  5,  5,  3,  3, 2, 1))
result4 <- gacha.sim(c( 10,  3,  2,  2, 1, 1))

シミュレーションの結果をヒストグラムにしてみましょう。

library(MASS)
windows(12, 9)
par(mfrow=c(2, 2))

x <- result1
truehist(x$result, xlim=c(min(x$result), max(x$result)),
         xlab="コンプまでの回数", ylab="人数", prob=FALSE)
legend(30,600,
       legend=c("出現確率 [1, 1, 1, 1, 1, 1]",
                paste("平均", round(mean(x$result), 1)),
                paste("標準偏差", round(sd(x$result), 1)),
                paste("中央値", median(x$result)),
                paste("最大値", max(x$result))))
par(new=TRUE)
plot(density(x$result), xlim=c(min(x$result), max(x$result)), col="red",
     lwd=2, xlab="", ylab="", axes=FALSE, main="")

x <- result2
truehist(x$result, xlim=c(min(x$result), max(x$result)),
         xlab="コンプまでの回数", ylab="人数", prob=FALSE)
legend(500,800,
       legend=c("出現確率 [100, 50, 10, 10, 3, 1]",
                paste("平均", round(mean(x$result), 1)),
                paste("標準偏差", round(sd(x$result), 1)),
                paste("中央値", median(x$result)),
                paste("最大値", max(x$result))))
par(new=TRUE)
plot(density(x$result), xlim=c(min(x$result), max(x$result)), col="red",
     lwd=2, xlab="", ylab="", axes=FALSE, main="")

x <- result3
truehist(x$result, xlim=c(min(x$result), max(x$result)),
         xlab="コンプまでの回数", ylab="人数", prob=FALSE)
legend(75,600,
       legend=c("出現確率 [5, 5, 3, 3, 2, 1]",
                paste("平均", round(mean(x$result), 1)),
                paste("標準偏差", round(sd(x$result), 1)),
                paste("中央値", median(x$result)),
                paste("最大値", max(x$result))))
par(new=TRUE)
plot(density(x$result), xlim=c(min(x$result), max(x$result)), col="red",
     lwd=2, xlab="", ylab="", axes=FALSE, main="")

x <- result4
truehist(x$result, xlim=c(min(x$result), max(x$result)),
         xlab="コンプまでの回数", ylab="人数", prob=FALSE)
legend(75,1000,
       legend=c("出現確率 [10, 3, 2, 2, 1, 1]",
                paste("平均", round(mean(x$result), 1)),
                paste("標準偏差", round(sd(x$result), 1)),
                paste("中央値", median(x$result)),
                paste("最大値", max(x$result))))
par(new=TRUE)
plot(density(x$result), xlim=c(min(x$result), max(x$result)), col="red",
     lwd=2, xlab="", ylab="", axes=FALSE, main="")

描画ももっとスマートな方法がありそう。

軸の目盛もちゃんと合わせた方がいいよなぁ。


ちなみに実行速度は

system.time(result1 <- gacha.sim(c(  1,  1,  1,  1, 1, 1)))
# ユーザ   システム       経過  
#     4.94       0.33       5.27 

system.time(result2 <- gacha.sim(c(100, 50, 10, 10, 3, 1)))
# ユーザ   システム       経過  
#    31.11       0.20      31.47 

system.time(result3 <- gacha.sim(c(  5,  5,  3,  3, 2, 1)))
# ユーザ   システム       経過  
#    5.99       0.20       6.20 

system.time(result4 <- gacha.sim(c( 10,  3,  2,  2, 1, 1)))
# ユーザ   システム       経過  
#    6.85       0.22       7.07 

ちょっと思ったこととして、我々がガキの頃1回20円で回してたガチャガチャは非復元抽出だったわけだけど、今のソーシャルゲームのガチャは復元抽出なんだよね。

そのあたりが感覚的にわかりにくくて文句言う人が出てきてるのかも。


なお、coupon collector's problemについてはこちらのびいだまブログさんの記事も参考になります。


【追記】

期待通りにビッキーさんがやってくれました!

d:id:a_bicky:20120318:1332083264

こちらに載っているgacha.sim2関数は爆速です!

一応、私の環境で試してみると、

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    6.27 2.125424      5.95     0.29         NA        NA
# 2 gacha.sim2(p)            1    2.95 1.000000      2.95     0.00         NA        NA

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   44.72  17.0038     44.35     0.33         NA        NA
# 2 gacha.sim2(p)            1    2.63   1.0000      2.62     0.00         NA        NA

すげー。

3倍から22倍も早くなってる。

r-de-rr-de-r 2012/04/13 16:20 ご指名があったことに気づきましたので,書いて見ました

r-de-rr-de-r 2012/04/15 22:42 グラフを簡単に描く方法は,明らかと思いますが...
results <- list(
result1 = gacha.sim(c( 1, 1, 1, 1, 1, 1)),
result2 = gacha.sim(c(100, 50, 10, 10, 3, 1)),
result3 = gacha.sim(c( 5, 5, 3, 3, 2, 1)),
result4 = gacha.sim(c( 10, 3, 2, 2, 1, 1))
)
library(MASS)
windows(12, 9)
par(mfrow = c(2, 2))
draw <- function(x) {
truehist(x$result, xlim = c(min(x$result), max(x$result)),
xlab = "コンプまでの回数", ylab = "人数", prob = FALSE)
legend("topright", legend = c(paste("出現確率 [", paste(x$prob,
collapse = ", "), "]"), paste("平均", round(mean(x$result),
1)), paste("標準偏差", round(sd(x$result), 1)), paste("中央値",
median(x$result)), paste("最大値", max(x$result))))
par(new = TRUE)
plot(density(x$result), xlim = c(min(x$result), max(x$result)),
col = "red", lwd = 2, xlab = "", ylab = "", axes = FALSE,
main = "")
}
sapply(results, draw)