Monti Hall 問題をPythonで実験する

会社で採用の時に確率とか統計に関連した問題を出してみるみたいな話をしていたときにMonti Hall問題の話題がでて、以前これを理解するため（納得するため？）にPythonで実験してみて面白かったので記事にしてみました。

## Monti Hall問題
Monti Hall問題とは以下のような問題です。

```
あなたは景品がもらえるゲームに参加しています。
３つのドアがあって、１つには新車が、残りの２つにはやぎが入っています。
仮にあなたが１番のドアを選んだとしましょう。
すると司会者が３番のドアを開けてそこにはやぎが入っていました。
そして彼はこう問います。「（最初に選んだ１番から）２番のドアに変えますか？」
さて、２番のドアに変えたほうが良いのでしょうか？
（英語Wikiの訳）
```

[英語Wiki](https://en.wikipedia.org/wiki/Monty_Hall_problem)
[日本語Wiki](https://ja.wikipedia.org/wiki/%E3%83%A2%E3%83%B3%E3%83%86%E3%82%A3%E3%83%BB%E3%83%9B%E3%83%BC%E3%83%AB%E5%95%8F%E9%A1%8C)

この問題の答えは「２番に変えた場合の新車を選ぶ確率は2/3で、変えない場合の確率は1/3なので、必ず２番に変えたほうが良い」です。
これが答えなのですが正直答えを聞いても「なるほど確かに」とはなりにくい。少なくとも僕は直感的には理解できませんでした。実際多くの人が直感的に理解しにくいそうで、それゆえにパラドックスとして知られています。
数学的には条件付き確率、ベイズ統計を使った説明がWikiでされていますが、そんなことより実験してみようというのがこの記事の趣旨です。

## 実験用Pythonコード

以下のようなPythonのプログラムで実験してみます。

```python
"""

monti hall problem simulator

"""

import random
from functools import reduce

CAR = 'tesla model 3'
GOAT = 'goat'

def set_doors():
    doors = [CAR, GOAT, GOAT]
    random.shuffle(doors)
    return doors


def open_door(doors, choice):
    return doors[choice]


def game_without_monti_hall():
    "pick one door randomly"
    doors = set_doors()
    choice = random.randint(0, 2)
    return (doors, choice)


def game_with_monti_hall():
    "pick one door randomly then switch after monti choice"
    doors = set_doors()
    first_choice = random.randint(0, 2)
    monti_choice = [n for n in range(len(doors)) if n is not first_choice and
                    open_door(doors, n) is not CAR][0]
    second_choice = [n for n in range(len(doors))
                     if n is not first_choice and n is not monti_choice][0]

    return (doors, second_choice)


def eval_game(game, verbose=False):
    "evaluate game and return prize"
    game_res = game()
    prize = open_door(*game_res)
    if verbose:
        print(f'picked door #{game_res[1]}: prize was {prize}')
    return prize


def calculate_game_stats(game, n_try, verbose=False):
    "try game n_try times and calculate stats"
    prize_lst = (eval_game(game, verbose) for _ in range(n_try))
    def calc_gamestat(stat, prize):
        if prize is CAR:
            return (stat[0] + 1, stat[1])
        else:
            return (stat[0], stat[1] + 1)

    init_stat = (0, 0)  # car=0, goat=0
    res = reduce(calc_gamestat, prize_lst, init_stat)
    return res
```

プログラム自体はたいして複雑なものではないので、説明は省きます。`game_without_monti_hall`は単純に３つのドアから１つ選ぶ場合で、`game_with_monti_hall`がMonti Hallがドアを開けた後に、最初のドアから別のドアに変える場合です。

## 実験結果

先ほど作ったPythonのプログラムを実行してみました。

```python
In [5]: calculate_game_stats(game_without_monti_hall, 100000)
Out[5]: (33180, 66820)

In [6]: calculate_game_stats(game_with_monti_hall, 100000)
Out[6]: (66695, 33305)
```

１０万回ゲームをしてみて、車を当てた回数とヤギを当てた回数を計算しています。tupleの0個目が車を当てた回数、tupleの1個目がヤギを当てた回数です。
これで実験として妥当と言い切れるのかわかりませんが、最初の答えとよく合致した結果になっています。みなさんはどうかわかりませんが、僕はこれを計算してみた後は最初の答えにとても納得できるようになりました。（理解の方は置いておいて、気持ち的な問題です。）

## 余談
僕がこの問題を最初に知ったのは「ラスベガスをぶっつぶせ（原題：21）」というブラックジャックのカードカウンティングでカジノ荒らしをするMITの学生の映画でした。かなり脚色されていると思いますが（映画の主人公はもちろん白人なのですが、実際は台湾人のグループだったそうです）、実際にそういうグループがあったらしくそのメンバーの１人が"The House Advantage: Playing the Odds to Win Big in Business"という本を書いています。映画の方は数学とか確率は小道具程度の扱いですが、こちらの本は確率やそれをビジネスやギャンブルにどう活かすかといったことに興味がある人にはとてもおすすめです。