モンティ・ホール問題とその派生形による意思決定のシミュレーション

めっちゃ有名なあれ

　モンティ・ホール問題はヒューリスティックな意思決定と確率論的な正解が異なる有名な例題で、以下のようなクイズ形式のゲームを考える(リンクより引用)。

モンティ・モール問題:
0: 三つの扉がある。一つは正解。二つは不正解。
1:挑戦者は三つの中から一つ扉を選ぶ。
2:司会者(モンティ)は答えを知っており、残り二つの扉の中で不正解の扉を一つ選んで開ける。
3:挑戦者は残りの二つの扉の中から好きな方を選べる。このとき扉を変えるべきか？変えないべきか？
モンティ・ホール問題とその解説 | 高校数学の美しい物語

　この問題の直観的な回答は「残った二つの扉のいずれかが正解なのだから、その確率は変えても変えなくても正解の確率は1/2だ。だから変えるべきとも変えないべきとも言えない」となると思う*1。ところがこの問題に対する確率論的な正解は「扉を変えずに正解である確率＝1/3、扉を変えて正解である確率=2/3」であり、より良い意思決定は「扉を変えるべき」となる。この正解は直観では納得しづらく、実際この正解がサヴァント*2によって世に広く知れ渡った際、多くの数学者から反論が寄せられ大論争となったらしい。データ分析、あるいは統計・確率論的なアプローチに期待されることの一つに「直観に反する洞察」があるので、今回はこのモンティ・ホール問題から派生させた実際に展開できる問題を設定し、意思決定の流れを検討してみたい*3。

モンティ・ホール問題の正解の確認

　この問題の解説は既に多くのWebサイトでも示されているが、せっかくなので式展開とシミュレーションで簡単に確認しておく。まずは式展開から。

ベイズの定理による確認

　 $A, \ B, \ C$ の三つの扉について、それぞれが正解である事象を $A_{当} }, \ B_{当}, \ C_{当$ とする。いま挑戦者は $A$ を選び、モンティは残り二つの扉のうち $C$ を開けたとしよう(これを事象 $C_{開}$ とする)。問題は、「 $C$ が開けられたという条件の下で $A$ が正解である確率 $p(A_{当} | C_{開})$ 」と「 $C$ が開けられたという条件の下で $B$ が正解である確率 $p(B_{当} | C_{開})$ 」をそれぞれ見積もることである。それぞれ以下のように求まる*4。

${\displaystyle \begin{eqnarray} p(A_{当}|C_{開}) &=& \frac{p(C_{開}|A_{当})p(A_{当})} {p(C_{開})} \\ &=& \frac{p(C_{開}|A_{当})p(A_{当})} {p(C_{開}|A_{当})p(A_{当}) + p(C_{開}|B_{当})p(B_{当}) + p(C_{開}|C_{当})p(C_{当})} \\ &=& \frac {\frac{1}{2}・ \frac{1}{3}} {\frac{1}{2} ・ \frac{1}{3} + 1 ・ \frac{1}{3} + 0 ・ \frac{1}{3}} \\ &=& 1/3 \end{eqnarray} }$

${\displaystyle \begin{eqnarray} p(B_{当}|C_{開}) &=& \frac{p(C_{開}|B_{当})p(B_{当})} {p(C_{開})} \\ &=& \frac{p(C_{開}|B_{当})p(B_{当})} {p(C_{開}|A_{当})p(A_{当}) + p(C_{開}|B_{当})p(B_{当}) + p(C_{開}|C_{当})p(C_{当})} \\ &=& \frac {1・ \frac{1}{3}} {\frac{1}{2} ・ \frac{1}{3} + 1 ・ \frac{1}{3} + 0 ・ \frac{1}{3}} \\ &=& 2/3 \end{eqnarray} }$
　なお、上記では各扉に正解が割り当てられる確率や、挑戦者が最初に選んだ扉が正解だった場合残った二つの扉のどちらかが開けられる確率はランダムかつ等確率であるという前提の下で数値を代入していることに注意されたい。これより、挑戦者が最初に選んでいた扉 $A$ よりも扉 $B$ の方が当たる可能性が高いことが示された。

シミュレーションによる確認

　コードはこんな感じ。簡単のため挑戦者が最初に選ぶ扉は $A$ で固定している。クイズを1,000回やってみて扉を変更しなかった際の正解率と扉を変更した際の正解率を試算している。

set.seed(123)

N <- 1000

door_open <- function(N){
  door <- c("A", "B", "C")
  result <- data.frame(matrix(NA, nrow = N, ncol = 5))
  # 順に正解、モンティが開いた扉、挑戦者が最初に選んだ扉(Aで固定)、挑戦者が扉を変更すると判断した際に選ばれる扉、正解か否か
  colnames(result) <- c("answer", "opened_door", "select", "re_select", "correct_answer")
  
  # 正解の扉
  for(i in 1:nrow(result)){
    result$answer[i] <- sample(door, size = 1, prob = c(1/3, 1/3, 1/3))
  }
  
  # シミュレーションでは回答者は最初常にAを選択するものとする
  result$select <- "A"
  
  # 回答後開かれる扉をＢかＣから決定する
  for(i in 1:nrow(result)){
    if(result$answer[i] == "A"){
      result$opened_door[i] <- sample(c("B", "C"), size = 1, prob = c(1/2, 1/2))
    }else if(result$answer[i] == "B"){
      result$opened_door[i] <- "C"
    }else{
      result$opened_door[i] <- "B"
    }
  }
  
  # 扉が開かれた後、回答者が回答を変更した際の扉
  for(i in 1:nrow(result)){
    if(result$opened_door[i] == "B"){
      result$re_select[i] <- "C"
    }else{
      result$re_select[i] <- "B"
    }
  }
  
  # 最初の回答が正解の場合1、変更した回答が正解の場合2を返す
  for(i in 1:nrow(result)){
    if(result$answer[i] == result$select[i]){
      result$correct_answer[i] <- 1
    }else{
      result$correct_answer[i] <- 2
    }
  }
  
  return(result)
}

sim_result <- door_open(N)
prop.table(table(sim_result$correct_answer))

> prop.table(table(sim_result$correct_answer))

    1     2 
0.332 0.668

　扉を変更しなかった際の正解率(1:0.332)に対し変更した際の正解率(2:0.668)が約2倍と、こちらでも「扉を変えるべき」という結論が得られた。

モンティ・ホール問題の派生形の設定

　今回は実際の意思決定に応用するために、モンティ・ホール問題に対し以下の変更・追加設定を行い検討を進めてみる。

=======================================

クイズの運営は挑戦者やモンティとは別の組織が担っており、モンティは司会役として雇われている身である
そのため実際のクイズではモンティも正解の扉を知らず、司会進行に努めるのみである。公平性の観点から扉の正否や過去実績に関する一切の情報は組織から提供されない
挑戦者はモンティが正解の扉を知らないということを知らない
モンティは運営組織から「挑戦者の正解率を下げる」ことを求められており、一定期間ごとに評価が行われる。評価期間内での正解率が1/2を下回っていた場合にそれに応じた報酬を受け取り、逆に1/2を上回っていた場合、報酬は受け取れず解雇リスクが増大する(すなわちモンティには正解率を最小化するインセンティブが働く)。
つい最近、コラムニストのサヴァントが公の場で「モンティ・ホール・クイズは扉を変更したほうが当たる確率が高い！」と発表した
発表後、モンティはどうも挑戦者が扉を変更する可能性が高まっているように感じている
もし本当に挑戦者が扉を変更する可能性が高まっているのなら、モンティは司会の中で扉を変更しないよううまく誘導し、正解率を下げたい
モンティはサヴァントの発表前後で行われたクイズ各10回、計20回の結果を記憶している

=======================================

　この設定の下でモンティの気持ちを考えてみる。モンティも正解の扉は知らされていないため、モンティにとっては常に「挑戦者に扉を変えないよう仕向ける」ことが最善手となる。持っているデータを使って、モンティは答えるべき問いに対する回答を得て司会進行のあり方を検討する、というストーリー。
　なお、今回の検討では古典的な仮説検定等ではなく、ベイズ的なアプローチをとる。かなりおおらかに説明すると、頻度論的なアプローチでは推定したい母数は固定された真の値を持っているとするのに対し、ベイズ的なアプローチでは母数そのものが確率変数として分布していると考える。推定したい母数に対して設定した分布(事前分布)と与えられたデータをもって母数を繰り返しサンプリングし、その分布をもって母数を表現しよう(事後分布)というのが基本的なアプローチである。いくつか利点があるが、意思決定の観点から言えば、まどろっこしい手続きを一旦脇に置いて*5、人が確率に対して直観的に持つ感覚に近い形で情報を得られるという点が挙げられる。ベイズ統計に関する参考文献は以下などが分かりやすい。

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門

作者: 豊田秀樹
出版社/メーカー: 朝倉書店
発売日: 2015/06/25
メディア: 単行本
この商品を含むブログ (6件) を見る

StanとRでベイズ統計モデリング (Wonderful R)

作者: 松浦健太郎,石田基広
出版社/メーカー: 共立出版
発売日: 2016/10/25
メディア: 単行本
この商品を含むブログ (10件) を見る

また、古典的な統計学とベイズ統計の違いについては以下の記事が分かりやすい。
healthpolicyhealthecon.com

　以上の前提、設定の下で、モンティは以下のように検討を進めた。

モンティの検討①挑戦者が扉を変更する可能性は高まっている(≒対策を講じるべき)と考えて良いか？

　モンティはとにかくも今現在の扉を変更する確率( $p$ とする)を把握し、それが以前よりも高くなっているとどの程度信じて良いのかどうか知る必要がある。設定条件によるとモンティはサヴァント発表前後各10回のクイズ結果を覚えているため、これを使って答えを出したい。なお、当然最終的に関心のある結果は挑戦者のクイズの正否であるが、クイズに正解する確率を $q$ とすれば、

${\displaystyle \begin{eqnarray} q &=& (最初の選択が正解である確率) × (扉を変更しない確率) + (最初の選択が外れである確率) × (扉を変更する確率) \\ &=& \frac{1}{3}(1-p) × \frac{2}{3}p \\ &=& \frac{1}{3}(1 + p) \end{eqnarray} }$

　であるので、扉を変更する確率と正に対応していることが分かる。
　さて、手元にあるデータ( Monty_Bayes - Google ドライブ )を使って推定を試みる。サンプリングにはstanファイルを別途用意する。コードはこんな感じ。
●Stanのコード

data{
  int<lower=0> N[2];                                  // サヴァント発表前後のクイズ回数。ともに10回
  int n[2, 2];                                        // 発表前後×正否数の配列
}

parameters{
  simplex[2] p[2];                                    // 発表前後×正否確率の配列
}

model{
  for(i in 1:2){                                      // モンティ・ホール解説前後 
    for(j in 1:2){                                    // 回答を変更する/しない
      n[i, j] ~ binomial(N[i], p[i][j]);              // 一回の試行における正否は二項分布に従うと仮定
    }
  }
}

generated quantities{
  real delta;                                         // 発表前後における扉変更確率の差
  real<lower=0, upper=1> p_delta_over;                // 発表後の変更確率が発表前の変更確率を上回る確率
  real q_before;                                      // 発表前の正解確率
  real q_after;                                       // 発表後の正解確率 
  delta = p[2][1] - p[1][1];
  p_delta_over = step(delta);
  q_before = (1.0/3) * p[1][2] + (2.0/3) * p[1][1];
  q_after = (1.0/3) * p[2][2] + (2.0/3) * p[2][1];
}

●Rのコード

library(tidyverse)
library(readxl)
library(rstan)

dat_before <- read_excel("./data.xlsx", sheet = "phase1")  # 発表前の実績
dat_after <- read_excel("./data.xlsx", sheet = "phase2")   # 発表後の実績

N <- c(nrow(dat_before), nrow(dat_after))
n <- matrix(c(sum(dat_before$change), sum(dat_after$change),
              (N[1] - sum(dat_before$change)), (N[2] - sum(dat_after$change))), nrow = 2)

data <- list(N = N, n = n)

fit <- stan(file = "./code1.stan", data = data, iter = 12000, warmup = 2000, chains = 4, thin = 1)

traceplot(fit)
fit

●推定結果

> fit
Inference for Stan model: code1.
4 chains, each with iter=12000; warmup=2000; thin=1; 
post-warmup draws per chain=10000, total post-warmup draws=40000.

               mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
p[1,1]         0.41    0.00 0.10   0.22   0.34   0.41   0.48   0.62 36078    1
p[1,2]         0.59    0.00 0.10   0.38   0.52   0.59   0.66   0.78 36078    1
p[2,1]         0.68    0.00 0.10   0.48   0.62   0.69   0.75   0.86 37469    1
p[2,2]         0.32    0.00 0.10   0.14   0.25   0.31   0.38   0.52 37469    1
delta          0.27    0.00 0.14  -0.01   0.18   0.28   0.37   0.54 36611    1
p_delta_over   0.97    0.00 0.17   0.00   1.00   1.00   1.00   1.00 28345    1
q_before       0.47    0.00 0.03   0.41   0.45   0.47   0.49   0.54 36078    1
q_after        0.56    0.00 0.03   0.49   0.54   0.56   0.58   0.62 37469    1
lp__         -29.68    0.01 1.04 -32.45 -30.07 -29.36 -28.94 -28.67 17717    1

Samples were drawn using NUTS(diag_e) at Wed Sep 26 11:57:32 2018.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

　10,000回×4回で計40,000回パラメータをサンプリングした。これによるとサヴァント発表前の扉変更確率が41％であるのに対し、発表後は68％となっている。シミュレーションの97％では発表後が発表前を上回る変更確率であると推定されたので、かなりの確信をもって「サヴァント発表後、挑戦者は扉を変更する可能性が高まっている」と結論付けて良いように思える*6。モンティはこの結果を材料にして、司会中に施策を打つことを運営組織に願い出ることにした。

モンティの検討②「扉を変えない方が良いよ」と助言することはモンティにとって有益か？

　何等かの対策が必要であることが分かった。前述したようにモンティも正解の扉を知らない(ゆえに挑戦者の最初の回答の正否が分からない)ため、常に扉を変更しないよう仕向けることが最善である。ここで、挑戦者が「モンティは正解の扉を知らない、ということを知らない」のであれば、挑戦者がモンティの言うことを信じる確率の大小によっては、「扉を変えない方が良いよ」と助言することが有効に働く可能性が出てくる*7。そこでモンティは運営組織を説得し、司会中にこの助言を行うことが許されることとなった。助言開始後10回の結果をもって挑戦者がモンティの助言を信じる確率( $θ$ とする)を推定し、この施策が有効かどうか確かめたい。先ほどと同様に推定できる。
●Stanのコード

data{
  int<lower=0> N;
  int n[2];
}

parameters{
  simplex[2] theta;
}

model{
  for(i in 1:2){
    n[i] ~ binomial(N, theta[i]);
  }
}

generated quantities{
  real delta;
  real<lower=0, upper=1> delta_over;
  real<lower=0, upper=1> q;
  delta = theta[1] - 0.32;
  delta_over = step(delta);
  q = (1.0/3) * (2 - theta[1]);
}

●Rのコード

dat_lie <- read_excel("./data.xlsx", sheet = "phase3")  # 助言開始後10回の実績

N2 <- nrow(dat_lie)
n_trust <- c(sum(dat_lie$trust), N2 - sum(dat_lie$trust))

data2 <- list(N = N2, n = n_trust)

fit2 <- stan(file = "./code2.stan", data = data2, iter = 12000, warmup = 2000, chains = 4, thin = 1)

traceplot(fit2)
fit2

●推定結果

> fit2
Inference for Stan model: code2.
4 chains, each with iter=12000; warmup=2000; thin=1; 
post-warmup draws per chain=10000, total post-warmup draws=40000.

             mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
theta[1]     0.41    0.00 0.10   0.22   0.34   0.41   0.48   0.62 13806    1
theta[2]     0.59    0.00 0.10   0.38   0.52   0.59   0.66   0.78 13806    1
delta        0.09    0.00 0.10  -0.10   0.02   0.09   0.16   0.30 13806    1
delta_over   0.80    0.00 0.40   0.00   1.00   1.00   1.00   1.00 16819    1
q            0.53    0.00 0.03   0.46   0.51   0.53   0.55   0.59 13806    1
lp__       -15.39    0.01 0.72 -17.44 -15.56 -15.12 -14.94 -14.88 17975    1

Samples were drawn using NUTS(diag_e) at Wed Sep 26 13:45:43 2018.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

　挑戦者がモンティの助言を信じる確率は41％と推定された。ここで0.09と推定されているデルタ $δ$ は助言実施時に挑戦者がそれを信じる確率と、助言実施前の挑戦者が扉を変更しない確率の差であり、 $δ\_over$ はその差が正となる(助言の実施が実施しない場合に比べて有効である)確率である。導出過程は以下のようになる。
${\displaystyle 助言実施時の挑戦者の正解確率q'は、 \\ \begin{eqnarray} q' &=& (最初の選択が正解である確率)×(助言を信じる確率) + (最初の選択が外れである確率)×(助言を信じない確率) \\ &=& \frac{1}{3} θ + \frac{2}{3} (1 - θ) \\ &=& \frac{1}{3} (2 - θ) \end{eqnarray} }$
${\displaystyle 助言により挑戦者の正解率が下がれば良い(q' < q)ので、 \\ \begin{eqnarray} \frac{1}{3} (2 - θ) &<& \frac{1}{3} (1 + p) \\ \therefore 1 - p &<& θ \end{eqnarray} }$
${\displaystyle これより、δを以下で定義する。 \\ δ = θ - (1 - p) }$
　今回は検討①で得た扉を変更しない確率の平均値0.32をそのまま代入した。推定結果によると、挑戦者がモンティの助言を信じる確率は41％で、これにより挑戦者の正解確率は多少下がるものの、依然53％と確率1/2を上回っている。モンティは助言をより多くの挑戦者に信じさせる追加の対策か、まったく別の対策を講じる必要があることが分かった。

おわりに

　今回は有名な問題を通してベイズ統計の基本的なアプローチを追ってみた。このようにベイズ的アプローチでは、頻度論的なアプローチにとってつく「帰無仮説は採択できない」「信頼区間はその中に真の値が入る確率ではない」等のまどろっこしい概念やサンプルサイズを気にするといったことがなく、今あるデータから得られるとりあえずの解を人間の感覚に沿った確率の概念で求められる、というところが意思決定時には魅力的だと思う。一方でよく言われるように、ベイズ統計は主観を入れ込む余地があるため(先の例で言えば「挑戦者の80％はモンティの言うことを信じる気がする」などと言って正解率を推定することも出来る)、利用には注意が必要である。

*1:モンティ・ホール問題と直観的な思考過程について分かりやすく図解されたものは例えば以下: analytics-notty.tech

*2:史上最もIQが高い人として有名らしい

*3:ただしデータ分析による意思決定支援においては、「直観でもわかる知識はいらない」とされる一方で「直観に反する知識は受け入れられない」というジレンマが頻発する。この問題は実際上極めて重要な問題なので、また稿を改めて検討したい

*4:条件付き確率やベイズの定理等の解説は例えばこちらの記事が分かりやすいatarimae.biz

*5:仮説検定や信頼区間は初学時にはまずもって誤解する

*6:今回は議論の簡単化のため信用区間の検討は行わない

*7:モンティの助言と扉の正否が独立であるため、挑戦者はモンティの助言の有無から正解を推定できない、純粋にモンティを信じるかどうかが問われる、というタテツケ