library(conflicted) # 関数名の衝突を警告
conflicts_prefer( # 優先的に使う関数を指定
dplyr::filter(),
dplyr::select(),
dplyr::lag(),
)
library(tidyverse)第6講 確率シミュレーション
乱数を用いた数値実験
準備
以下で利用する共通パッケージを読み込む．
コイン投げの賭け
問題
以下のようなコイン投げの賭けを考える．
- Alice と Bob の二人で交互にコインを投げ，最初に表が出た方を勝ちとする．
この賭けの勝率を求めるための確率シミュレーションを行いなさい．
解答欄
解答例
試行を行う関数を定義する．
my_trial <- function(){
while(TRUE){ # 永久に回るループ
if(sample(c("h","t"),1) == "h"){ # h:head，t:tail
return("Alice") # Aliceが表を出して終了
}
if(sample(c("h","t"),1) == "h"){
return("Bob") # Bobが表を出して終了
}
#' どちらも裏ならもう一度ループ
}
}
Note
この関数は rbinom を用いても実装できるので試みてみよう
試行を行ってみる. このとき勝った方の名前が表示される．
my_trial(); my_trial(); my_trial() # 3回行ってみる[1] "Bob"[1] "Alice"[1] "Bob"確率シミュレーションを行ってみる．
#' set.seed(8888) # 実験を再現する場合は適当なシードを指定する
mc <- 10000 # 回数を設定
my_data <- replicate(mc, my_trial())
glimpse(my_data) # 結果が mc 次元のベクトルになっていることが確認できる chr [1:10000] "Alice" "Bob" "Bob" "Bob" "Alice" "Bob" "Alice" "Alice" ...簡単な集計を行う．
table(my_data) # 頻度を表示するmy_data
Alice Bob
6703 3297 table(my_data)/mc # 勝率の計算 (全実験回数で除算)my_data
Alice Bob
0.6703 0.3297 この賭けは先手が有利であることが確認できる．
確率シミュレーションの例題
以下の確率的な事象のシミュレーションを考えてみなさい．
Buffon の針
2次元平面上に等間隔 d で平行線が引いてある． 長さ l の針を この平面上にランダムに落としたとき， 平行線と交わる確率はいくつか？ ただし l\leq d とする．
解答欄
解答例
針を投げる試行を定義する． 試行の周期性から1本の線の周りで問題を考えれば良い． 針の中心位置 x\;(-d/2<x<d/2) と向き \theta\;(0<\theta<2\pi) をランダムに生成する． 針の先と尾の座標の符号が異なれば位置0の線と交わっていると判断できる．
buffon_trial <- function(d, l, verbose = FALSE){ # dとlを指定
x <- runif(1, min = -d/2, max = d/2) # 位置
theta <- runif(1, min = 0, max = 2*pi) # 角度
cross <- FALSE # 交わったかどうかを示す変数
if((x+l*cos(theta)/2)*(x-l*cos(theta)/2) < 0){
cross <- TRUE # 交わった場合に書き換え
}
if(verbose == TRUE){ # 位置と角度も返す
return(c(x = x, theta = theta, # ベクトルで返すので
cross = as.numeric(cross))) # データ型を揃える
} else { # 交わったかどうかだけ返す
return(cross) # 論理値のまま返す
}
}試行を行ってみる．
d <- 10
l <- 7
buffon_trial(d, l); buffon_trial(d, l)[1] FALSE[1] TRUE詳細な情報を表示する場合 (絵を描くときに用いる) は以下のように実行する．
buffon_trial(d, l, verbose = TRUE) x theta cross
-2.8125240 0.6066636 1.0000000 buffon_trial(d, l, verbose = TRUE) x theta cross
1.5230809 0.5476163 1.0000000 試行の結果を図にする．
mc <- 10
replicate(mc, buffon_trial(d, l, verbose = TRUE)) |>
t() |> # 配列を転置
as_tibble() |> # tibble に変換
mutate( # 図に必要な情報を加える
y = runif(mc, min = -(d+l)/2, max = (d+l)/2), # 中心のy座標をランダムに生成
x1 = x-l/2*cos(theta),
x2 = x+l/2*cos(theta),
y1 = y-l/2*sin(theta),
y2 = y+l/2*sin(theta),
cross = as.logical(cross)) |> # 交わりを表すラベル(論理値)
ggplot() +
geom_vline(xintercept = c(-d,0,d)) + # 縦棒を描く
geom_segment(aes(x = x1, y = y1, # 始点
xend = x2, yend = y2, # 終点
colour = cross)) + # 交わりを表示
labs(x = "x", y = "y") +
coord_fixed() # 座標軸の比を1に固定
x と \theta の関係を図にする．
mc <- 3000
replicate(mc, buffon_trial(d, l, verbose = TRUE)) |>
t() |>
as_tibble() |>
mutate( # 図に必要な情報を加える
cross = as.logical(cross)) |> # 交わり表すラベル(論理値)
ggplot(aes(x = x, y = theta, colour = cross)) +
geom_vline(xintercept = c(-d/2,d/2)) + # 境界を描く
geom_hline(yintercept = c(0,2*pi)) + # 横線を描く
geom_point() + labs(y = expression(theta)) 
大規模な確率シミュレーションを行う．
mc <- 100000 # 回数を設定
buffon_data <- replicate(mc, buffon_trial(d, l)) 簡単な集計を行い，理論値と比較する．
table(buffon_data) # 頻度 (TRUEが針の交わった回数)buffon_data
FALSE TRUE
55473 44527 table(buffon_data)/mc # 確率(推定値)buffon_data
FALSE TRUE
0.55473 0.44527 print((2*l)/(pi*d)) # 針の交わる確率 (理論値)[1] 0.44563382*l/d/(table(buffon_data)/mc)["TRUE"] # π の推定値 TRUE
3.14416 Monty Hall 問題
ゲームの参加者の前に閉まった3つのドアがあって， 1つのドアの後ろには景品の新車が， 2つのドアの後ろには外れを意味するヤギがいる． 参加者は新車が置かれたドアを当てると新車がもらえる．
参加者が1つのドアを選択した後， 司会のモンティが残りのドアのうちヤギがいるドアを開けてヤギを見せる． ここで参加者は，最初に選んだドアを残っているドアに変更してもよいと言われる．
参加者はドアを変更すべきだろうか？
解答欄
解答例
クイズに答える試行を定義する． 以下は賞品の位置は固定して解答者をランダムに配置する場合を想定している．
monty_trial <- function(verbose = FALSE){
#' 賞品は1に置いてあるものとする
choice <- sample(1:3, size = 1) # 解答者の最初の選択
if(choice == 1) { # 変えないのが正解の場合
win <- "stay"
door <- sample(2:3, size = 1) # Monty が開ける扉 (2か3をランダムに選択)
} else { # 変えるのが正解の場合
win <- "change"
if(choice == 2) { # 2を選択したら3を開ける
door <- 3
} else { # そうでなければ(3を選択)2を開ける
door <-2
}
}
if(verbose == TRUE){ # 賞品，解答者の選択，Monty が開ける扉を返す
return(c(prize = 1, choice = choice, door = door))
} else { # 勝ち負けの条件を返す
return(win)
}
}賞品と解答者をランダムに配置する場合は以下のように定義すればよい． if 文で素朴に書くこともできるが， 集合を操作する関数 setdiff(), union() を利用している．
monty_trial <- function(verbose = FALSE){
prize <- sample(1:3, size = 1) # 賞品の置かれた扉
choice <- sample(1:3, size = 1) # 解答者の最初の選択
if(prize == choice) { # 変えないのが正解の場合
win <- "stay"
door <- sample(setdiff(1:3, prize), size = 1) # Monty が開ける扉 (ランダムに選択)
} else { # 変えるのが正解の場合
win <- "change"
door <- setdiff(1:3, union(prize, choice)) # Monty が開ける扉
}
if(verbose == TRUE){ # 賞品，解答者の選択，Monty が開ける扉を返す
return(c(prize = prize, choice = choice, door = door))
} else { # 勝ち負けの条件を返す
return(win)
}
}試行を行ってみる．賞品を獲得できた選択(stay/change)が返ってくる．
monty_trial(); monty_trial()[1] "change"[1] "change"詳細な情報を表示する場合 (絵を描くときに用いる) は以下のようになる． 表示は順に，賞品の位置(prize)，最初の選択(choice)，開かれたドア(door)を示している．
monty_trial(verbose = TRUE) prize choice door
1 1 2 monty_trial(verbose = TRUE) prize choice door
1 3 2 実験の状況を絵にする．記号は前述のとおり．
mc <- 15
replicate(mc, monty_trial(verbose = TRUE)) |>
t() |> # 配列を転置
as_tibble() |> # データフレームに変更
rowid_to_column(var = "index") |> # 試行の番号を追加
pivot_longer(!index) |> # index 以外を long format 化
ggplot(aes(x = index, y = value,
colour = name, shape = name)) +
geom_point(size = 4) + # サイズを大きめに表示
scale_y_continuous(breaks = 1:3)
実験とともに勝率がどのように変化するか図示する．
mc <- 400
replicate(mc, monty_trial()) |>
as_tibble_col(column_name = "win") |>
rowid_to_column(var = "index") |>
mutate(stay_win = win == "stay", # 留まって勝つ場合
ratio = cumsum(stay_win)/index) |> # 勝率の推移を計算
ggplot(aes(x = index, y = ratio)) +
geom_line(colour = "blue") +
geom_hline(yintercept = 1/3, colour = "orange") + # 理論値
ylim(c(0,1)) # y軸の描画範囲を指定
大規模な確率シミュレーションを実行する．
mc <- 50000 # 回数を設定
monty_data <- replicate(mc, monty_trial()) 簡単な集計を行う．
table(monty_data) # 頻度monty_data
change stay
33346 16654 table(monty_data)/mc # 確率(推定値)monty_data
change stay
0.66692 0.33308 選択を変えた場合の方が賞品を獲得する確率が高いことがわかる．
試行を模擬する関数は，論理式を使って書くこともできる． 関数 sample(x, size) では x に1つの数字 n を渡すと 1:n と解釈されるので，以下では x に文字を渡すようにしている．
monty_trial <- function(change = FALSE, # ドアを変えるか否か
verbose = FALSE){
doors <- LETTERS[1:3] # ドアの名前 A,B,C
prize <- sample(doors, 1) # 賞品のドア
choice <- sample(doors, 1) # 選んだドア
monty <- sample(doors[(doors != prize) & (doors != choice)], 1) # モンティが開いたドア
if(change){ # 選択を変えた場合の最後に選んだドア
final <- doors[(doors != choice) & (doors != monty)]
} else { # 選択を変えない場合
final <- choice
}
if(verbose){ # TRUEの場合
return(c(prize=prize,choice=choice,monty=monty,final=final))
} else {
return(ifelse(prize == final,"win","loss"))
}
}この関数を用いて，簡単な集計を行う．
table(replicate(mc, monty_trial()))/mc # 最初に選んだドアのままの場合
loss win
0.66494 0.33506 table(replicate(mc, monty_trial(change = TRUE)))/mc # 選んだドアを変えた場合
loss win
0.33244 0.66756 秘書問題 (最適停止問題)
以下の条件のもと秘書を1人雇うとする．
n 人が応募しており n は既知とする．
応募者には 1 位から n 位まで順位付けできる．
無作為な順序で1人ずつ面接を行う．
毎回の面接後その応募者を採用するか否かを決定する．
不採用にした応募者を後から採用することはできない．
“r-1 番までの応募者は採用せず， r 番以降の応募者でそれまで面接した中で最も良い者を採用する” という戦略を取るとき，最適な r はいくつだろうか?
解答欄
解答例
秘書の採用を模擬する試行を定義する．
secretary_trial <- function(n, r, verbose = FALSE){ # nとrを指定
applicants <- sample(1:n, size = n) # n人の順位をランダムに並べ替える
ref <- applicants[1:(r-1)] # r-1人目までは参照(最高順位が基準点)
test <- applicants[r:n] # r人目以降が採用候補
idx <- which(test < min(ref)) # 採用候補の中で基準点より良い人を選出
if(length(idx) == 0) { # 1人も基準より良い人がいない場合
employed <- applicants[n] # 最後の人を採用
} else {
employed <- test[idx[1]] # 最初に基準点を越えた人を採用
}
if(verbose == TRUE){ # 全順位も返す
return(list(applicants = applicants,
employed = employed)) # 2つの異なるベクトルをリストで束ねる
} else { # 採用した者の順位のみ返す
return(employed)
}
}条件を変えながら試行を行ってみる． applicants は候補者の順位， employed は戦略に従って採用された人の順位を示す．
n <- 10 # 候補者は10名
secretary_trial(n, 2, verbose = TRUE) # 2人目から採用を考える$applicants
[1] 9 8 3 5 10 4 7 1 2 6
$employed
[1] 8secretary_trial(n, 3, verbose = TRUE) # 3人目から採用を考える$applicants
[1] 6 2 4 8 1 3 9 7 10 5
$employed
[1] 1secretary_trial(n, 4, verbose = TRUE) # 4人目から採用を考える$applicants
[1] 10 6 2 9 5 1 4 3 8 7
$employed
[1] 1secretary_trial(n, 5, verbose = TRUE) # 5人目から採用を考える$applicants
[1] 3 10 6 4 5 8 2 1 9 7
$employed
[1] 2secretary_trial(n, 6, verbose = TRUE) # 6人目から採用を考える$applicants
[1] 4 7 8 6 9 5 2 10 3 1
$employed
[1] 2大規模なシミュレーションを実行する． 適当な r に対して，採用された人の順位の頻度表を作成している．
#' set.seed(8888) # 実験を再現したい場合はシードを指定
mc <- 3000
n <- 30 # 候補者数を変えて実験
secretary_data <- tibble(r = NULL,
employed = NULL)
for (r in 2:(n-1)) {
foo <- replicate(mc, secretary_trial(n, r))
if(r %in% c(2,6,10,14,18,22)) { # いくつか表示
cat("採用開始: ", r, "\n")
print(table(foo))
}
secretary_data <- bind_rows(secretary_data, # 前の実験結果に追加
tibble(r = rep(r,mc),
employed = foo))
}採用開始: 2
foo
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
400 324 237 220 202 183 155 119 109 105 112 99 74 84 76 83 59 59 51 43
21 22 23 24 25 26 27 28 29 30
33 39 34 38 14 16 12 7 8 5
採用開始: 6
foo
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
916 540 356 230 170 122 87 73 56 38 44 26 28 20 23 31 19 18 19 17
21 22 23 24 25 26 27 28 29 30
18 21 20 16 18 14 17 18 15 10
採用開始: 10
foo
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
1098 494 271 151 105 69 57 39 36 37 32 35 26 25 29 34
17 18 19 20 21 22 23 24 25 26 27 28 29 30
45 36 35 33 30 28 33 26 28 38 35 37 33 25
採用開始: 14
foo
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
1120 405 171 103 60 37 41 47 48 46 48 41 37 42 52 44
17 18 19 20 21 22 23 24 25 26 27 28 29 30
39 51 59 61 43 53 53 37 40 61 38 45 40 38
採用開始: 18
foo
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
993 283 132 66 60 58 67 56 64 62 56 56 49 55 63 41 59 52 61 50
21 22 23 24 25 26 27 28 29 30
47 57 61 72 60 55 69 69 60 67
採用開始: 22
foo
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
758 191 89 73 68 77 61 73 78 59 76 69 85 84 63 72 75 57 79 71
21 22 23 24 25 26 27 28 29 30
72 75 68 78 72 65 74 82 72 84 各 r でどのような順位が採用されたか分布を図示する． 上記の頻度表の箱ひげ図を r ごとに描いたものである．
secretary_data |>
ggplot(aes(x = r, y = employed)) +
geom_boxplot(aes(group = r), # rごとに集計
fill = "white", colour ="royalblue") +
labs(title = paste("n =", n)) # nをタイトルに表示
理論的に良いとされるrの値 (nが十分大きい場合に求めることができるので調べてみよ) を計算する．
n/exp(1) [1] 11.03638各 r で1位を採用できる確率を図示する．
secretary_data |>
group_by(r) |>
summarize(ratio = mean(employed == 1)) |> # 1位ならTRUE(1)
ggplot(aes(x = r, y = ratio)) +
geom_step(colour = "royalblue") + # 階段関数で描画
geom_vline(xintercept = n/exp(1), colour = "red") + # 理論値
labs(title = paste("n =", n)) # nをタイトルに表示