じゃんけんのシミュレーション

人員が異動で減って、一時的に仕事量が倍になっている。

あと２週間の辛抱だ。

お世話になっている驚異のアニヲタ社会復帰への道さんが、卒業記念に本を書いたとのこと。サイトが怪しいが、５００円とちょっとで安い。おすすめである。

購入はここ＝>　[Rで始めた医学・統計学・Bioinformatics]

この本の中で初心者の頃に書いたという、じゃんけんのシミュレーションの

コードが有る。これが相当時間がかかるので改造してみる。

グリコ、パイナツプル、チヨコレイトのあれだ。

ランダムに相手が手を出してくるとき、そのような比率で手を出せば

勝てる可能性が高くなるか。

じゃんけんのシミュレーション

# シミュレーションによって勝率を計算する
library(MCMCpack)  # ディリクレ分布による確率を発生させるのに必要なパッケージ

library(doSNOW)  # 並列計算に必要なパッケージ

library(foreach)  # 並列計算に必要なパッケージ

library(gplots)  # 色付けに必要なパッケージ

library(scatterplot3d)

niter <- 10000  # シミュレーション回数
m <- c(3, 6, 6)  # グー，チョキ、パー

x <- rdirichlet(niter, alpha = rep(1, length(m)))  # Xのじゃんけんの手の確率 alpha=1で一様分布
y <- rdirichlet(niter, alpha = rep(1, length(m)))  # Yのじゃんけんの手の確率　alpha10で真ん中に手厚い分布
yy <- t(y)
# じゃんけんごとに前進できる歩数行列
M <- matrix(c(0, -m[1], m[3], m[1], 0, -m[2], -m[3], m[2], 0), nr = 3, nc = 3)
cl <- makeCluster(20, "SOCK")  # 並列計算しているので，不用意にやると計算が終わらない。
registerDoSNOW(cl)
win <- foreach(i = seq(niter), .combine = "c") %dopar% {
    mean( (x[i, ] %*% M %*% yy) >= 0)  #ここを行列計算に変更する
}
stopCluster(cl)
colcut <- 100
# cols <- rainbow(colcut, start=.4, end=.35) #こっちの色の方が好き
cols <- colorpanel(colcut, low = "blue", mid = "white", high = "red")
idx <- cut(as.numeric(win), seq(0, 1, length = colcut + 1))
layout(matrix(c(1, 2, 1, 1), nr = 2), width = c(1, 2.5), heights = c(6, 1))
sc <- scatterplot3d(x, color = cols[idx], pch = 16, angle = 150, xlab = "Rock", 
    ylab = "Scissors", zlab = "Paper", cex.lab = 2, cex.axis = 1.2, mar = c(4, 
        3, 2, 3))
sc$points3d(x[which.max(win), 1], x[which.max(win), 2], x[which.max(win), 3], 
    pch = "*", col = "yellow", cex = 3)
sc$points3d(x[which.min(win), 1], x[which.min(win), 3], x[which.min(win), 3], 
    pch = "*", col = "green", cex = 3)
par(mar = c(4, 1, 0, 0))
image(as.matrix(sort(win)), col = cols, frame = FALSE, yaxt = "n")
mtext("Win probability", side = 3, line = 0.5, cex = 1.5)

行列計算化して、オリジナルより速度はかなりはやくなった。

コードをなぞりながら改変することにより、並列計算なるものを

はじめて経験出来た。これは勉強になった。

dirichlet分布も面白い。alpha=1で一様分布のような感じ。

alphaを大きくすると、中心に固まる。全部グーとか、全部パーとかは
あまりなさそうなので、アルファを大きくした場合が現実に近いかも
しれないと思った。