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] "Alice"[1] "Alice"[1] "Bob"確率シミュレーションを行ってみる．
#' set.seed(8888) # 実験を再現する場合は適当なシードを指定する
mc <- 10000 # 回数を設定
my_data <- replicate(mc, my_trial())
glimpse(my_data) # 結果が mc 次元のベクトルになっていることが確認できる chr [1:10000] "Alice" "Alice" "Alice" "Alice" "Alice" "Alice" "Alice" ...簡単な集計を行う．
table(my_data) # 頻度を表示するmy_data
Alice Bob
6729 3271 table(my_data)/mc # 勝率の計算 (全実験回数で除算)my_data
Alice Bob
0.6729 0.3271 この賭けは先手が有利であることが確認できる．
確率シミュレーションの例題
以下の確率的な事象のシミュレーションを考えてみなさい．
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] TRUE[1] FALSE詳細な情報を表示する場合 (絵を描くときに用いる) は以下のように実行する．
buffon_trial(d, l, verbose = TRUE) x theta cross
1.123461 4.945397 0.000000 buffon_trial(d, l, verbose = TRUE) x theta cross
-1.542298 2.978795 1.000000 試行の結果を図にする．
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
55565 44435 table(buffon_data)/mc # 確率(推定値)buffon_data
FALSE TRUE
0.55565 0.44435 print((2*l)/(pi*d)) # 針の交わる確率 (理論値)[1] 0.44563382*l/d/(table(buffon_data)/mc)["TRUE"] # π の推定値 TRUE
3.15067 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] "stay"[1] "change"詳細な情報を表示する場合 (絵を描くときに用いる) は以下のようになる． 表示は順に，賞品の位置(prize)，最初の選択(choice)，開かれたドア(door)を示している．
monty_trial(verbose = TRUE) prize choice door
1 1 2 monty_trial(verbose = TRUE) prize choice door
2 1 3 実験の状況を絵にする．記号は前述のとおり．
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
33204 16796 table(monty_data)/mc # 確率(推定値)monty_data
change stay
0.66408 0.33592 選択を変えた場合の方が賞品を獲得する確率が高いことがわかる．
試行を模擬する関数は，論理式を使って書くこともできる． 関数 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.66416 0.33584 table(replicate(mc, monty_trial(change = TRUE)))/mc # 選んだドアを変えた場合
loss win
0.33148 0.66852 秘書問題 (最適停止問題)
以下の条件のもと秘書を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 1 4 3 2 6 8 5 7 10
$employed
[1] 1secretary_trial(n, 3, verbose = TRUE) # 3人目から採用を考える$applicants
[1] 7 8 1 9 10 3 4 2 6 5
$employed
[1] 1secretary_trial(n, 4, verbose = TRUE) # 4人目から採用を考える$applicants
[1] 4 9 5 10 1 7 2 3 6 8
$employed
[1] 1secretary_trial(n, 5, verbose = TRUE) # 5人目から採用を考える$applicants
[1] 7 10 2 9 5 8 3 6 4 1
$employed
[1] 1secretary_trial(n, 6, verbose = TRUE) # 6人目から採用を考える$applicants
[1] 7 9 4 5 8 6 10 2 1 3
$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
407 303 277 211 194 160 168 124 138 127 93 93 87 84 76 65 53 52 42 36
21 22 23 24 25 26 27 28 29 30
45 31 37 29 22 14 9 15 4 4
採用開始: 6
foo
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
922 549 312 233 176 133 92 74 60 58 34 33 26 13 17 21 21 21 18 22
21 22 23 24 25 26 27 28 29 30
17 21 22 13 16 11 11 24 15 15
採用開始: 10
foo
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
1135 493 257 164 97 68 48 44 37 30 39 38 31 23 30 21
17 18 19 20 21 22 23 24 25 26 27 28 29 30
31 29 28 30 26 36 40 36 28 31 35 31 32 32
採用開始: 14
foo
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
1121 406 153 110 66 53 42 47 54 43 55 46 46 40 44 53
17 18 19 20 21 22 23 24 25 26 27 28 29 30
39 44 48 42 43 44 38 50 41 46 33 45 62 46
採用開始: 18
foo
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
994 320 126 87 56 57 64 53 47 65 49 55 61 52 62 59 48 57 55 59
21 22 23 24 25 26 27 28 29 30
67 51 42 69 57 68 44 49 66 61
採用開始: 22
foo
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
815 181 97 74 70 86 61 63 70 65 69 67 70 63 54 72 79 67 63 65
21 22 23 24 25 26 27 28 29 30
76 71 90 70 76 76 87 70 62 71 各 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をタイトルに表示