library(conflicted) # 関数名の衝突を警告
conflicts_prefer( # 優先的に使う関数を指定
dplyr::filter(),
dplyr::select(),
dplyr::lag(),
)
library(tidyverse)統計データ解析I
第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] "Bob"[1] "Alice"確率シミュレーションを行ってみる．
#' 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
6659 3341 table(my_data)/mc # 勝率の計算 (全実験回数で除算)my_data
Alice Bob
0.6659 0.3341 この賭けは先手が有利であることが確認できる．
確率シミュレーションの例題
以下の確率的な事象のシミュレーションを考えてみなさい．
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] FALSE詳細な情報を表示する場合 (絵を描くときに用いる) は以下のように実行する．
buffon_trial(d, l, verbose = TRUE) x theta cross
2.2400693 0.8048767 1.0000000 buffon_trial(d, l, verbose = TRUE) x theta cross
-0.01850161 5.00294411 1.00000000 試行の結果を図にする．
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
55220 44780 table(buffon_data)/mc # 確率(推定値)buffon_data
FALSE TRUE
0.5522 0.4478 print((2*l)/(pi*d)) # 針の交わる確率 (理論値)[1] 0.44563382*l/d/(table(buffon_data)/mc)["TRUE"] # π の推定値 TRUE
3.126396 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] "stay"詳細な情報を表示する場合 (絵を描くときに用いる) は以下のようになる． 表示は順に，賞品の位置(prize)，最初の選択(choice)，開かれたドア(door)を示している．
monty_trial(verbose = TRUE) prize choice door
1 2 3 monty_trial(verbose = TRUE) prize choice door
1 2 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
33239 16761 table(monty_data)/mc # 確率(推定値)monty_data
change stay
0.66478 0.33522 選択を変えた場合の方が賞品を獲得する確率が高いことがわかる．
試行を模擬する関数は，論理式を使って書くこともできる． 関数 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.66344 0.33656 table(replicate(mc, monty_trial(change = TRUE)))/mc # 選んだドアを変えた場合
loss win
0.33532 0.66468 秘書問題 (最適停止問題)
以下の条件のもと秘書を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] 6 3 4 7 8 2 5 1 9 10
$employed
[1] 3secretary_trial(n, 3, verbose = TRUE) # 3人目から採用を考える$applicants
[1] 5 2 1 7 4 8 10 6 3 9
$employed
[1] 1secretary_trial(n, 4, verbose = TRUE) # 4人目から採用を考える$applicants
[1] 1 6 9 4 5 2 3 10 8 7
$employed
[1] 7secretary_trial(n, 5, verbose = TRUE) # 5人目から採用を考える$applicants
[1] 3 8 10 9 7 4 5 6 2 1
$employed
[1] 2secretary_trial(n, 6, verbose = TRUE) # 6人目から採用を考える$applicants
[1] 5 8 7 9 10 6 2 3 1 4
$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
387 303 239 224 182 165 142 139 137 130 88 103 95 99 71 80 63 66 49 44
21 22 23 24 25 26 27 28 29 30
51 27 18 27 23 10 17 13 7 1
採用開始: 6
foo
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
940 525 371 240 141 113 91 71 61 33 40 30 38 20 23 16 19 24 20 23
21 22 23 24 25 26 27 28 29 30
14 17 11 15 15 19 21 13 18 18
採用開始: 10
foo
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
1128 501 239 138 112 65 63 51 46 42 38 25 32 36 27 33
17 18 19 20 21 22 23 24 25 26 27 28 29 30
31 43 26 36 22 25 31 28 32 24 25 22 38 41
採用開始: 14
foo
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
1134 388 158 88 56 61 54 46 39 43 55 55 37 52 42 49
17 18 19 20 21 22 23 24 25 26 27 28 29 30
28 42 42 57 37 55 45 50 53 40 44 46 55 49
採用開始: 18
foo
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
1006 267 114 77 71 64 51 48 62 55 53 60 60 59 75 52
17 18 19 20 21 22 23 24 25 26 27 28 29 30
56 56 72 66 54 56 57 58 56 58 66 59 57 55
採用開始: 22
foo
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
761 178 86 78 75 77 56 57 69 93 78 74 86 61 81 81 79 70 62 73
21 22 23 24 25 26 27 28 29 30
85 61 74 69 59 78 71 67 80 81 各 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をタイトルに表示