---
title: "統計データ解析I"
subtitle: "第6講 練習問題 解答例"
date: "`r Sys.time()`"
format:
    html: 
      toc: true
      html-math-method: katex
      self-contained: true
      grid: 
        margin-width: 350px
execute: 
  echo: true
  warning: false
reference-location: margin
citation-location: margin
tbl-cap-location: margin
fig-cap-location: margin
editor: visual
editor_options: 
  chunk_output_type: console
---

## 準備

以下で利用する共通パッケージを読み込む．

```{r}
library(conflicted)  # 関数名の衝突を警告
conflicts_prefer(    # 優先的に使う関数を指定
    dplyr::filter(),
    dplyr::select(),
    dplyr::lag(),
    )
library(tidyverse)
```

## コイン投げの賭け

### 問題

以下のようなコイン投げの賭けを考える．

-   Alice と Bob の二人で交互にコインを投げ，最初に表が出た方を勝ちとする．

この賭けの勝率を求めるための確率シミュレーションを行いなさい．

### 解答例

試行を行う関数を定義する．

```{r}
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が表を出して終了
        }
        #' どちらも裏ならもう一度ループ
    }
}
```

::: callout-note
この関数は rbinom を用いても実装できるので試みてみよう
:::

試行を行ってみる. このとき勝った方の名前が表示される．

```{r}
my_trial(); my_trial(); my_trial() # 3回行ってみる
```

確率シミュレーションを行ってみる．

```{r}
#' set.seed(8888) # 実験を再現する場合は適当なシードを指定する
mc <- 10000       # 回数を設定 
my_data <- replicate(mc, my_trial())
glimpse(my_data)      # 結果が mc 次元のベクトルになっていることが確認できる
```

簡単な集計を行う．

```{r}
table(my_data)    # 頻度を表示する
table(my_data)/mc # 勝率の計算 (全実験回数で除算)
```

この賭けは先手が有利であることが確認できる．

## 確率シミュレーションの例題

以下の確率的な事象のシミュレーションを考えてみなさい．

### Buffon の針

2次元平面上に等間隔 $d$ で平行線が引いてある． 長さ $l$ の針を この平面上にランダムに落としたとき， 平行線と交わる確率はいくつか？ ただし $l\leq d$ とする．

### 解答例

針を投げる試行を定義する． 試行の周期性から1本の線の回りで問題を考えれば良い． 針の中心位置 $x\;(-d/2<x<d/2)$ と向き $\theta\;(0<\theta<2\pi)$ をランダムに生成する． 針の先と尾の座標の符号が異なれば位置0の線と交わっていると判断できる．

```{r}
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) # 論理値のまま返す
    }
}
```

試行を行ってみる．

```{r}
d <- 10
l <- 7
buffon_trial(d, l); buffon_trial(d, l)
```

詳細な情報を表示する場合 (絵を描くときに用いる) は以下のように実行する．

```{r}
buffon_trial(d, l, verbose = TRUE)
buffon_trial(d, l, verbose = TRUE)
```

試行の結果を図にする．

```{r}
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$ の関係を図にする．

```{r}
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)) 
```

大規模な確率シミュレーションを行う．

```{r}
mc <- 100000 # 回数を設定 
buffon_data <- replicate(mc, buffon_trial(d, l)) 
```

簡単な集計を行い，理論値と比較する．

```{r}
table(buffon_data)      # 頻度 (TRUEが針の交わった回数)
table(buffon_data)/mc   # 確率(推定値)
print((2*l)/(pi*d))     # 針の交わる確率 (理論値)
2*l/d/(table(buffon_data)/mc)["TRUE"] # π の推定値
```

### Monty Hall 問題

ゲームの参加者の前に閉まった3つのドアがあって， 1つのドアの後ろには景品の新車が， 2つのドアの後ろには外れを意味するヤギがいる． 参加者は新車が置かれたドアを当てると新車がもらえる．

参加者が1つのドアを選択した後， 司会のモンティが残りのドアのうちヤギがいるドアを開けてヤギを見せる． ここで参加者は，最初に選んだドアを残っているドアに変更してもよいと言われる．

参加者はドアを変更すべきだろうか？

### 解答例

クイズに答える試行を定義する． 以下は賞品の位置は固定して解答者をランダムに配置する場合を想定している．

```{r}
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()` を利用している．

```{r}
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)が返ってくる．

```{r}
monty_trial(); monty_trial()
```

詳細な情報を表示する場合 (絵を描くときに用いる) は以下のようになる． 表示は順に，賞品の位置(prize)，最初の選択(choice)，開かれたドア(door)を示している．

```{r}
monty_trial(verbose = TRUE)
monty_trial(verbose = TRUE)
```

実験の状況を絵にする．記号は前述のとおり．

```{r}
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)
```

実験とともに勝率がどのように変化するか図示する．

```{r}
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軸の描画範囲を指定
```

大規模な確率シミュレーションを実行する．

```{r}
mc <- 50000 # 回数を設定 
monty_data <- replicate(mc, monty_trial()) 
```

簡単な集計を行う．

```{r}
table(monty_data)    # 頻度
table(monty_data)/mc # 確率(推定値)
```

選択を変えた場合の方が賞品を獲得する確率が高いことがわかる．

試行を模擬する関数は，論理式を使って書くこともできる． 関数 `sample(x, size)` では `x` に1つの数字 `n` を渡すと `1:n` と解釈されるので，以下では `x` に文字を渡すようにしている．

```{r}
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"))
    }
}
```

この関数を用いて，簡単な集計を行う．

```{r}
table(replicate(mc, monty_trial()))/mc              # 最初に選んだドアのままの場合
table(replicate(mc, monty_trial(change = TRUE)))/mc # 選んだドアを変えた場合
```

### 秘書問題 (最適停止問題)

以下の条件のもと秘書を1人雇うとする．

-   $n$ 人が応募しており $n$ は既知とする．

-   応募者には $1$ 位から $n$ 位まで順位付けできる．

-   無作為な順序で1人ずつ面接を行う．

-   毎回の面接後その応募者を採用するか否かを決定する．

-   不採用にした応募者を後から採用することはできない．

"$r-1$ 番までの応募者は採用せず， $r$ 番以降の応募者でそれまで面接した中で最も良い者を採用する" という戦略を取るとき，最適な $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` は戦略に従って採用された人の順位を示す．

```{r}
n <- 10 # 候補者は10名
secretary_trial(n, 2, verbose = TRUE) # 2人目から採用を考える
secretary_trial(n, 3, verbose = TRUE) # 3人目から採用を考える
secretary_trial(n, 4, verbose = TRUE) # 4人目から採用を考える
secretary_trial(n, 5, verbose = TRUE) # 5人目から採用を考える
secretary_trial(n, 6, verbose = TRUE) # 6人目から採用を考える
```

大規模なシミュレーションを実行する． 適当な `r` に対して，採用された人の順位の頻度表を作成している．

```{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))
}
```

各 `r` でどのような順位が採用されたか分布を図示する． 上記の頻度表の箱ひげ図を `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が十分大きい場合に求めることができるので調べてみよ) を計算する．

```{r}
n/exp(1) 
```

各 `r` で1位を採用できる確率を図示する．

```{r}
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をタイトルに表示
```