【R】司会者の選ぶ箱が偏った時のモンティ・ホール問題

モンティ・ホール問題が統計検定1級で出題された

2023年の統計検定1級を受けてきましたー。残念ながら統計応用のみ合格ということで来年こそは1級合格を目指していきたいと思います...
驚いたことに今年の統計検定では、共通問題でモンティ・ホール問題が出題されました!
統計に詳しい方なら皆さんご存知とは思いますが、モンティ・ホール問題とは、

1. 3つの箱(A, B, C)があり、その中に当たりが1つ存在する。
2. 参加者が3つの中から1つを選ぶ。
3. 次に、当たりの選択肢を知っている司会者が、残りの2つの中からハズレのものを1つ選んで参加者に知らせる。
   この状況において、選択肢を変更した方が有利である(2/3の確率であたりを引くことができる)というのがモンティ・ホール問題です。

統計検定1級では、

1. 3つの箱(A, B, C)があり、その中に当たりが1つ存在する。
2. 参加者はAを選んだ。
   という状況の上で、司会者が開ける可能性がある箱が2つあり、Cを開ける確率が一般にwであるとき、箱Bが当たりである確率pをwを用いて表わせ。
   という問題が出題されました。基本的なシチュエーションではw=0.5であり、司会者側は開ける可能性がある箱が2つあるときにはランダムに選択するのですが、ここで箱に対する嗜好性がある場合に拡張されているわけです。
   解いている間、「よくまぁこんな設定を思いつくものだな」と思いつつ、なんやかんや論述していた記憶があります。(答え自体はあっていた)

ということで今回はこの統計検定の問題にちなんで、司会者がCを開けやすかったり開けにくかったりしたら、箱を変えるという選択自体常に最良と言えるか？ という問題を考えてみます。

シミュレータづくり

そこで、一般化された状態でのシミュレータを作っていきましょう。

まずはsetup

library(tidyverse)

そして、確率wと箱を用意

# 司会者が箱3を開ける確率w
w <- 0
boxes <- c("A", "B", "C")

次にざっと100000回ぐらい試行してもらいましょう。
普通の言語であればfor文を使うのですが、せっかくRを使っているのでベクトル化しましょう。
まずはあたりの箱と、参加者が最初に選んだ箱を用意します。

set.seed(42)
simulations <- tibble(
  id = 1:100000,
  あたりの箱 = map_chr(id, ~ sample(boxes, 1)),
  参加者が最初に選んだ箱 = map_chr(id, ~ sample(boxes, 1)),)

simulations

結果はこんな感じ。

> simulations
# A tibble: 100,000 × 3
    id あたりの箱 参加者が最初に選んだ箱
   <int> <chr>    <chr>
 1   1 A      B
 2   2 A      B
 3   3 A      C
 4   4 A      A
 5   5 B      A
 6   6 B      C
 7   7 B      B
 8   8 A      B
 9   9 C      A
10  10 C      B
# ℹ 99,990 more rows
# ℹ Use `print(n = ...)` to see more rows

さて、ここから司会者が箱を選択します。ここが少しややこしくいくつか場合分けする必要があります。

1. Cがあたりかつ参加者がCを選んだ。
   わかりやすいです。A or Bを0.5の確率で選択します。
2. Cがハズレかつ参加者がCを選んだ。
   AまたはBがあたりなので、あたりで無い方を確定で選択します。
3. Cがあたりかつ参加者がA or Bを選んだ。
   司会者はCを選ぶことができません。なので参加者が選ばなかったハズレの箱を確定で選択します。
4. Cがハズレかつ参加者がA or Bを選んだ。
   ここはさらに2つに分岐します
   1. あたりと参加者が選択した箱が同じ
      確率wで司会者が残った2つの箱を選択
   2. あたりと参加者が選択した箱が異なる
      Cの箱が自動的に選択されます

...といってもよくわからないので、こういうときは全部書くのが良い方法です。
全パターンを書いてみました。

あたりの箱	参加者が最初に選んだ箱	司会者の選択する箱
A	A	B or C (確率w)
A	B	C
A	C	B
B	A	C
B	B	A or C (確率w)
B	C	A
C	A	B
C	B	A
C	C	A or B

これを見るとかなりスッキリ記述できます。
つまり、

1. あたりの箱 ≠ 参加者が最初に選んだ箱 ならば 司会者は残りの箱を選択する
2. あたりの箱 = 参加者が最初に選んだ箱 = C ならば A or Bを0.5の確率で選択する
3. それ以外のときは確率wでCの箱を選択する。

ということで、これらをまとめた関数を準備します。

# 司会者が開ける箱を決定する関数
choose_host_box <- function(あたり, 参加者, boxes, w) 
{
  if (あたり != 参加者) 
  {
    setdiff(boxes, c(あたり, 参加者))
  } 
  else if (あたり == "C" && 参加者 == "C")
  {
    sample(setdiff(boxes, "C"), 1)
  } 
  else 
  {
    remaining_boxes = setdiff(boxes, c(あたり, 参加者))
    if (runif(1) < w) 
    {
      "C"
    } 
    else 
    {
      setdiff(remaining_boxes,"C")
    }
  }
}

simulations <- tibble(
  id = 1:100000,
  あたりの箱 = map_chr(id, ~ sample(boxes, 1)),
  参加者が最初に選んだ箱 = map_chr(id, ~ sample(boxes, 1)))

# シミュレーションデータに司会者が開ける箱を追加
simulations <- simulations |>
  mutate(司会者が開ける箱 = map2_chr(あたりの箱, 参加者が最初に選んだ箱, ~ choose_host_box(..1, ..2, boxes, w)))

simulations

# A tibble: 100,000 × 4
      id あたりの箱 参加者が最初に選んだ箱 司会者が開ける箱
   <int> <chr>      <chr>                  <chr>
 1     1 A          B                      C
 2     2 A          A                      B
 3     3 C          C                      A
 4     4 C          B                      A
 5     5 A          B                      C
 6     6 A          B                      C
 7     7 C          B                      A
 8     8 C          C                      A
 9     9 C          C                      A
10    10 A          A                      B
# ℹ 99,990 more rows
# ℹ Use `print(n = ...)` to see more rows

このあと、参加者が最初に選んだ箱からスイッチするかどうかを決め、最終的に当たり外れのフラグを立てます。

simulations <- simulations |> 
rowwise() |>
mutate(
    箱変更あり = setdiff(boxes, c(参加者が最初に選んだ箱, 司会者が開ける箱)),
    箱変更なし = 参加者が最初に選んだ箱,
    )|> 
mutate(
    箱変更あり_あたり = ifelse(あたりの箱 == 箱変更あり, TRUE, FALSE),
    箱変更なし_あたり = ifelse(あたりの箱 == 箱変更なし, TRUE, FALSE),
)

> simulations
# A tibble: 100,000 × 8
# Rowwise:
      id あたりの箱 参加者が最初に選んだ箱 司会者が開ける箱 箱変更あり
   <int> <chr>      <chr>                  <chr>            <chr>
 1     1 A          B                      C                A
 2     2 A          A                      B                C
 3     3 C          C                      A                B
 4     4 C          B                      A                C
 5     5 A          B                      C                A
 6     6 A          B                      C                A         
 7     7 C          B                      A                C
 8     8 C          C                      A                B
 9     9 C          C                      A                B
10    10 A          A                      B                C
# ℹ 99,990 more rows
# ℹ 3 more variables: 箱変更なし <chr>, 箱変更あり_あたり <lgl>,
#   箱変更なし_あたり <lgl>
# ℹ Use `print(n = ...)` to see more rows

これでシミュレーションができたので、関数としてまとめましょう

library(tidyverse)

# 司会者が開ける箱を決定する関数
choose_host_box <- function(あたり, 参加者, boxes, w) {
  if (あたり != 参加者) {
    setdiff(boxes, c(あたり, 参加者))
  } else if (あたり == "C" && 参加者 == "C") {
    sample(setdiff(boxes, "C"), 1)
  } else {
    remaining_boxes = setdiff(boxes, c(あたり, 参加者))
    if (runif(1) < w) {
      "C"
    } else {
      setdiff(remaining_boxes, "C")
    }
  }
}

# シミュレーションを実行する関数
run_simulation <- function(num_simulations, w) {
  set.seed(42)
  boxes <- c("A", "B", "C")

  simulations <- tibble(
    id = 1:num_simulations,
    あたりの箱 = map_chr(id, ~ sample(boxes, 1)),
    参加者が最初に選んだ箱 = map_chr(id, ~ sample(boxes, 1))
  ) |> 
    mutate(
      司会者が開ける箱 = map2_chr(あたりの箱, 参加者が最初に選んだ箱, ~ choose_host_box(..1, ..2, boxes, w))
    ) |> 
    rowwise() |> 
    mutate(
      箱変更あり = list(setdiff(boxes, c(参加者が最初に選んだ箱, 司会者が開ける箱)))[[1]],
      箱変更なし = 参加者が最初に選んだ箱
    ) |> 
    ungroup() |> 
    mutate(
      箱変更あり_あたり = ifelse(あたりの箱 == 箱変更あり, TRUE, FALSE),
      箱変更なし_あたり = ifelse(あたりの箱 == 箱変更なし, TRUE, FALSE),
    )

  return(simulations)
}

# 例：100000回のシミュレーションを実行
simulations_result <- run_simulation(100000, 0)

シミュレータの結果

ここで、パラメータwを変えていったとき、箱変更あり/なしの正解率を出してみましょう。

library(tidyverse)

# シミュレーション結果から正解率を計算する関数
calculate_accuracy <- function(simulations) {
  simulations %>%
    summarise(
      箱変更あり_正解率 = mean(箱変更あり_あたり, na.rm = TRUE),
      箱変更なし_正解率 = mean(箱変更なし_あたり, na.rm = TRUE)
    )
}

# 異なるwの値に対するシミュレーションの実行と正解率の計算
run_simulation_with_different_w <- function(num_simulations, w_values) {
  results <- map_dfr(w_values, ~{
    simulation_result <- run_simulation(num_simulations, .x)
    accuracy <- calculate_accuracy(simulation_result)
    tibble(w = .x, accuracy)
  })

  return(results)
}

# 例：異なるwの値でシミュレーションを実行し、正解率を計算
w_values <- seq(0, 1, by = 0.05)
simulation_results <- run_simulation_with_different_w(10000, w_values)
simulation_results

> simulation_results
# A tibble: 21 × 3
       w 箱変更あり_正解率 箱変更なし_正解率
   <dbl>             <dbl>             <dbl>
 1  0                0.655             0.345
 2  0.05             0.655             0.345
 3  0.1              0.655             0.345
 4  0.15             0.655             0.345
 5  0.2              0.655             0.345
 6  0.25             0.655             0.345
 7  0.3              0.655             0.345
 8  0.35             0.655             0.345
 9  0.4              0.655             0.345
10  0.45             0.655             0.345
# ℹ 11 more rows
# ℹ Use `print(n = ...)` to see more rows

意外ですが、司会者が箱Cを開けやすかろうが、開けにくなかろうが常に箱を変更したほうが良いという結果が得られます。
直感的には違和感がありますが、実は確率wが効いてくるパターンはあたりの箱 = 参加者が最初に選んだ箱のときのみであり、
このような結果が得られるわけです。

最後に

意外と知られている問題でも条件設定によってはまた違う趣があって面白いですね！