JAGS でモンティ・ホール問題を解く
JAGS と R 言語の私の入門・練習として、モンティ・ホール問題を JAGS (rjags) で解いてみる。
| ■ |
はじめに
|
モンティ・ホール問題については、↓が詳しい。この記事では、この本をただ「本」と言って参照する。
|
ちなみに以前に私はモンティ・ホール問題について詳しい記事(pdf)を書いている。そちらもご参照いただけたらと思う。
「モンティ・ホール問題」について詳しくは、これらをご参照いただくとして、早速、本題に入ろう。
その前に、R と JAGS、rjags のインストールはすでに済んでいるものとする。ちなみに Stan を将来インストールするつもりなら、R は C:/Program Files 下では、なくルート C:/ 以下にインストールしておいたほうがいい。
| ■ |
もっとも単純なモデル
|
モンティ・ホール問題は次のような問題とする。
「司会者モンティ・ホールのテレビのショー。ドアが三つありそのうちの一つに当たりが入っていて、残りははずれ。プレイヤーは一つのドアをまず選ばされる。プレイヤーはドア 1 を選んだ。しかし、ドアはまだ開けない。モンティは、残りのドアのうち、はずれのドアを選んでそれを開ける。その上で、プレイヤーは、最初に選んだドア 1 にするか、モンティが残したもう一つのドアを選びなおすか選べることにする。ドア 1 を選んだときの当たりの確率、選択を変えたときの当たりの確率はどれぐらいになるか。」
JAGS でもっとも単純にプログラムすると次のようになるだろう。 (monty_hall_1.jag とする。)
model {
d ~ dcat(pi)
pi[1] = 1.0/3
pi[2] = 1.0/3
pi[3] = 1.0/3
md ~ dcat(pi2)
pi2[1] = 1.0/2
pi2[2] = 1.0/2
m <- ifelse(d == 1, md + 1, ifelse(d == 2, 2, 3))
choose1 <- equals(d, 1)
choosem <- equals(d, m)
}
dcat はカテゴリカル分布で、最初のものは、pi 配列の割合に応じてここでは 1 から 3 までのどれかの数値を出してくる。
d は「当たり」が入っているドアを示す。プレイヤーは必ずドア 1 をまず選ぶとする。md は、ドア 1 が「当たり」のとき、モンティが 2 と 3 のどちらのドアを開けるかという「サイコロ」であり、1 のときドア 2 を 2 のときドア 3 を開ける。ちなみにドア 1 が外れでないときは、モンティは「当たり」でないほうのドアを開ける。モンティの残したドアが m になる。
choose1 はモンティが開けたあと、なおドア 1 をプレイヤーが選び続けた場合の「当たり」の確率を集計するためのもの、choosem はモンティが開けたあと、モンティが残したドア m をプレイヤーが選んだときの「当たり」の確率を集計するためのものである。
実行するには R 上で次のようにする。
library(rjags)
m <- jags.model('monty_hall_1.jag', list(), list(list(d = 1, md = 1)), n.chains=1)
update(m, 10000)
post.list <- coda.samples(m, c('choose1', 'choosem', 'pi'), n.iter=10000)
print(summary(post.list))
ここで update はバーンイン期間だそうだが、なぜこれがいるか、10000 に指定しているが、本当はどれぐらいやればいいのか、…といったあたりはまだ私はわかっていない。n.chains も普通は 4 だが、デバッグのしやすさもあってとりあえず 1 にしている。データを与える必要はないため list() を与えている。n.chains = 1 なので、1 個分、d と md の初期値を与えている。
結果は次のようになる。
Iterations = 10001:20000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 10000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
choose1 0.3347 0.4719 0.004719 0.004719
choosem 0.6653 0.4719 0.004719 0.004719
pi[1] 0.3333 0.0000 0.000000 0.000000
pi[2] 0.3333 0.0000 0.000000 0.000000
pi[3] 0.3333 0.0000 0.000000 0.000000
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
choose1 0.0000 0.0000 0.0000 1.0000 1.0000
choosem 0.0000 0.0000 1.0000 1.0000 1.0000
pi[1] 0.3333 0.3333 0.3333 0.3333 0.3333
pi[2] 0.3333 0.3333 0.3333 0.3333 0.3333
pi[3] 0.3333 0.3333 0.3333 0.3333 0.3333
モンティ・ホール問題の結果として正しく、choose1 の確率が 1/3、choose2 の確率が 2/3 におおよそなっている。とりあえず pi の値も見ているが、もちろん、それぞれ 1/3 のままである。
| ■ |
パターン２のモンティ・ホール問題
|
本の３章「ベイズ的モンティ」の最初に載っている「パターン２」のモンティ・ホール問題を考える。これは、モンティがドアを開けるとき、ランダムに選んでいて、たまたま開けたところが外れていたという問題である。次のような問題になる。
「プレイヤーは一つのドアをまず選ばされる。プレイヤーはドア 1 を選んだ。しかし、ドアはまだ開けない。モンティは、残りのドアのうち、ランダムに選んだドア 2 を開けると、それははずれだった。その上で、プレイヤーは、最初に選んだドア 1 にするか、モンティが残したもう一つのドア 3 を選びなおすか選べることにする。」
ここで大事なのは、逆にモンティがドアを開けて当たりだった場合については無視すべきこと。これをどう JAGS で表現するか？ 最初はその方法がわからず、data としてモンティがドアを開けたときはずれである例を 10000 個ぐらい与えて、そのとき推定された pi や pi2 から結果を導く…みたいなこともやってみたのだが、それは何か違う。
結局、思い付いたのが、sum の結果を observe するのを用い、pi や pi2 を変形させる方法である。
ソースは次のようになる。(monty_hall_2.jag とする。)
data {
zero <- 0
two <- 2
}
model {
d ~ dcat(pi)
md ~ dcat(pi2)
pi[1] ~ dbeta(1.0, 1.0)
pi[2] ~ dbeta(1.0, 1.0)
pi[3] ~ dbeta(1.0, 1.0)
pi2[1] ~ dbeta(1.0, 1.0)
pi2[2] ~ dbeta(1.0, 1.0)
zero ~ sum(equals(md + 1, d), 0) ## i.e. d !~ md + 1
m <- ifelse(md == 1, 3, 2)
two ~ sum(m, 0) ## i.e. 2 ~ m
choose1 <- equals(d, 1)
choosem <- equals(d, m)
}
dbeta というのはベータ分布だが、そのパラメータが α = 1.0, β = 1.0 というのは uniform 分布と同じである。pi や pi2 の事前分布はようするに「わからない」ということ。
次に、「当たり d はモンティの選んだドア md +1 ではない」という条件を zero ~ sum(…) の行が表している。
ここがいくつかのハックがある。まず、ある変数がある値になることを v1 ~ v2 などと書くことは JAGS ではできない。しかし、v1 ~ sum(v2, 0) ということはできる。sum は observe できると決まっているからである。
次に 0 ~ sum(v2, 0) と書くことはできない、~ の左側は変数でなければならないからである。しかし、一旦、v1 ~ sum(v2, 0) と書いた上で、上の data ブロックで v1 <- 0 と固定するのは受け容れられる。
さらに、d はすでに d ~ dcat(pi) で「代入」されているため、再び d ~ sum(v2, 0) という形式で書くことはできない。しかし、one ~ sum(equals(v2, d), 0)という形なら書くことができる。今回はそれを否定したいので、zero ~ sum(equals(v2 , d), 0) の形にした。
…というわけである。
あと、m は、md で選ばなかったほうである。
two ~ sum(m, 0) は上の説明でわかったと思うが、2 ~ m を表している。
実行するには R 上で次のようにする。
m <- jags.model('monty_hall_2.jag', list(), list(list(pi = c(0.1, 0.1, 0.1), pi2 = c(0.1, 0.1), d = 1, md = 2)), n.chains=1)
update(m, 10000)
post.list <- coda.samples(m, c('choose1', 'choosem', 'pi', 'pi2'), n.iter=10000)
print(summary(post.list))
結果は次のようになる。
Iterations = 11001:21000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 10000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
choose1 0.4969 0.5000 0.005000 0.006644
choosem 0.5031 0.5000 0.005000 0.006644
pi[1] 0.5288 0.2851 0.002851 0.004289
pi[2] 0.5374 0.2856 0.002856 0.004447
pi[3] 0.4440 0.2881 0.002881 0.004148
pi2[1] 0.4095 0.2868 0.002868 0.004362
pi2[2] 0.5982 0.2587 0.002587 0.003442
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
choose1 0.00000 0.0000 0.0000 1.0000 1.0000
choosem 0.00000 0.0000 1.0000 1.0000 1.0000
pi[1] 0.03240 0.2886 0.5450 0.7779 0.9788
pi[2] 0.03458 0.2987 0.5518 0.7871 0.9806
pi[3] 0.01748 0.1896 0.4169 0.6873 0.9658
pi2[1] 0.01353 0.1524 0.3687 0.6402 0.9593
pi2[2] 0.10269 0.3971 0.6248 0.8228 0.9813
「パターン２」のモンティ・ホール問題として正しく、choose1 の確率が 1/2、 choosem の確率が 1/2 におおよそなっている。
注目すべきは、pi と pi2 の事前確率が歪んでいることである。当たりを選ぶときは、ドア 3 が選ばれにくいように、モンティがドアを選ぶときはドア 3 が選ばれやすいようになっている。
それにもかかわらず、いや、そうだからこそ、choose1 と choosem の確率はほぼ等しく 1/2 となるのである。
| ■ |
標準的モンティ・ホール問題 再考
|
この「パターン２」と同様の戦略で「パターン１」を解いてみよう。モンティは規則正しくドアを選ぶが、選ばれたドアはドア 2 であったと固定してみる。
ソースは次のようになる。(monty_hall_3.jag とする。)
data {
two <- 2
}
model {
d ~ dcat(pi)
md ~ dcat(pi2)
pi[1] ~ dbeta(1.0, 1.0)
pi[2] ~ dbeta(1.0, 1.0)
pi[3] ~ dbeta(1.0, 1.0)
pi2[1] ~ dbeta(1.0, 1.0)
pi2[2] ~ dbeta(1.0, 1.0)
m <- ifelse(d == 1, md + 1, ifelse(d == 2, 2, 3))
two ~ sum(m, 0) ## i.e. 2 ~ m
choose1 <- equals(d, 1)
choosem <- equals(d, m)
}
実行するには R 上で次のようにする。
m <- jags.model('monty_hall_3.jag', list(), list(list(pi = c(0.1, 0.1, 0.1), pi2 = c(0.1, 0.1), d = 1, md = 2)), n.chains=1)
update(m, 10000)
post.list <- coda.samples(m, c('choose1', 'choosem', 'pi', 'pi2'), n.iter=10000)
print(summary(post.list))
結果は次のようになる。
Iterations = 11001:21000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 10000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
choose1 0.3367 0.4726 0.004726 0.007833
choosem 0.6633 0.4726 0.004726 0.007833
pi[1] 0.5009 0.2894 0.002894 0.004378
pi[2] 0.5586 0.2776 0.002776 0.004249
pi[3] 0.4391 0.2884 0.002884 0.004247
pi2[1] 0.5356 0.2831 0.002831 0.004445
pi2[2] 0.4646 0.2912 0.002912 0.004445
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
choose1 0.00000 0.0000 0.0000 1.0000 1.0000
choosem 0.00000 0.0000 1.0000 1.0000 1.0000
pi[1] 0.02807 0.2477 0.4960 0.7508 0.9742
pi[2] 0.04592 0.3294 0.5814 0.7997 0.9812
pi[3] 0.01607 0.1817 0.4106 0.6775 0.9645
pi2[1] 0.03605 0.2992 0.5488 0.7852 0.9787
pi2[2] 0.01780 0.2094 0.4473 0.7146 0.9718
正しく、choose1 は 1/3、choosem は 2/3 になっている。
注目すべきは pi と pi2 で、当たりはドア 2 が選ばれやすく、ドア 3 が選ばれにくくなり、モンティはドア 2 を選びやすくなっている。
| ■ |
結論
|
モンティ・ホール問題を「パターン１」「パターン２」について JAGS で解き、答えが正しいことを確かめた。
ただし、当たりを選ぶこと、モンティがドアを選ぶことの基準が「わからない」状況で、ある事実が発生したという段階を経れば、分布がすでに歪むことも確認された。
| ■ |
ライセンス
|
パブリックドメイン。 (数式のような小さなプログラムなので。)
自由に改変・公開してください。
| ■ |
配布物
|
今回の配布物は以下の zip ファイル。更新があれば下のリンクの中身は最新のものに置き変わっているはず。
GitHub にも登録してあるが、更新は、ここの更新のあと1日から3日遅れの予定。
過去のバージョンが欲しい場合は、GitHub から辿っていけばよい。
| 更新： | 2020-04-18 |
| 初公開： | 2020年04月18日 19:28:11 |
| 最新版： | 2020年04月18日 19:32:11 |
2020-04-18 19:28:13 (JST) in シミュレーション R | 固定リンク | コメント (1)
![[RSS]](/archive/rss.gif "コメント一覧のRSS"){kind=small}
コメント
初公開: jags_monty_hall-20200418.zip。バージョン 0.0.1。
jags_monty_hall-20200418.zip に相当するものは↓からどうぞ。
https://github.com/JRF-2018/jags_monty_hall/releases/tag/v0.0.1
[cocolog:91836101](↓)に感想も書いたのでよければそちらもお読みください。
http://jrf.cocolog-nifty.com/statuses/2020/04/post-6fbd1e.html
＞「JAGS でモンティ・ホール問題を解く」という記事を書いた。R 言語とベイズ統計モデリングの JAGS の練習として。まだ、JAGS の勉強をしようとしてぜんぜんしていない段階なので、これでいいのか自信がないが…。＜
投稿: JRF | 2020-04-18 20:10:16 (JST)