こんにちは。mutopsyです。この記事は，Stan Advent Calendar 2018の22日目のエントリー記事です。この記事では，森永製菓が発売しているおなじみのチョコレート菓子「チョコボール」で金のエンゼルが当たる確率と銀のエンゼルが当たる確率を，多項分布モデルを用いて同時推定し，銀のエンゼルが金のエンゼルの何倍出やすいかを明らかにすることを試みます。ちなみにこの記事は，以前書いた記事「金のエンゼルが当たる確率のベイズ推定」の続編という位置づけですが，前回の記事を読んでいなくても理解できるように書いたつもりです。前回と同様，チョコボール50周年を記念して行われた「金のエンゼル2倍キャンペーン」のデータを有効活用してベイズ推定を行います。
ちなみに結論から言えば以下のような推定結果が得られました：
- 金のエンゼルが当たる確率は0.02%─0.26% (MAPは0.09%；前回の結果のreproduce)
- 銀のエンゼルが当たる確率は3.28%─5.51% (MAPは4.29%)
- 珍しさという点で，金のエンゼル1枚は銀のエンゼル約11─159枚分の価値がある (MAPは33枚分)
2019/12/15追記：この記事の内容をさらに発展させた内容が書籍になりました。→「たのしいベイズモデリング2：事例で拓く研究のフロンティア」の第1章
目次
1. 目的と概要
2. エンゼルが当たる確率の生成メカニズムとしての多項分布モデル
3. 使用するデータおよびStanとRのコード
4. 金と銀のエンゼルが当たる確率の同時推定の結果
5. 銀のエンゼルは金のエンゼルの何倍出やすいか
1. 目的と概要
改めまして，こんにちは。チョコボールと早めのMCMCが効いたのか，先週から患っていた風邪が治りつつあるmutopsyです。前回の記事「金のエンゼルが当たる確率のベイズ推定」では，先人たちが収集したチョコボールデータを用いて二項分布モデルによるなんちゃってメタ分析を行い，金のエンゼルが当たる確率がMAPで約0.09% (99% HDI [0.02%, 0.26%]) であることを明らかにしました。続編となる今回の記事では，金のエンゼルが当たる確率と銀のエンゼルが当たる確率を同時に推定し，さらに銀のエンゼルが金のエンゼルの何倍出やすいのかを区間推定した結果を報告します。
1つの推定方法としては，前回と同じように二項分布モデルを用いて銀のエンゼルが出る確率 (金のエンゼルが出ず，はずれでもない確率) を推定するという手が考えられます。その推定結果を，前回得られた金のエンゼルが当たる確率の推定結果の点推定値 (0.09%) で除することで，銀のエンゼルが金のエンゼルの何倍出やすいかを計算することができます。ところが，この方法では「99%の確率で，銀のエンゼルは金のエンゼルよりも〇〜△倍出やすい」といった幅をもった推定 (区間推定) を行うことは困難です。そこで本記事では，金のエンゼルが当たる確率p_goldと銀のエンゼルが当たる確率p_silverを1つのモデルで同時に推定し，さらにこれらの比 ratio = p_gold / p_silver の生成量を計算します。このようなモデルを用いることで，金のエンゼルが当たる確率と銀のエンゼルが当たる確率の不確実性に関する情報 を保持したまま，その比の区間推定を行うことが可能です。このように，パラメータの事後分布が持つ情報を損失することなく差や比などの生成量の事後分布を簡単に推定できるのは，ベイズ推定の強みの1つです。
2. エンゼルが当たる確率の生成メカニズムとしての多項分布モデル
では，金のエンゼルが当たる確率と銀のエンゼルが当たる確率を同時に推定するにはどのようなモデルを用いれば良いのでしょうか。前回の記事では，あるチョコボールの箱のクチバシを確認したときに「金のエンゼルが当たる事象」と「金のエンゼルが当たらない事象」という2つの事象が起こり得るというベルヌーイ分布を仮定し，この試行を独立に複数回繰り返したときに金のエンゼルが当たる回数が二項分布に従うことを利用しました。しかし，実際には「金のエンゼルが当たる事象」「銀のエンゼルが当たる事象」「いずれのエンゼルも当たらない事象」という3つの事象が存在します。このように起こり得る事象が3つ以上存在する場合は，ベルヌーイ分布ではなく次のようにカテゴリカル分布を仮定するのが自然です。
p = {p_gold, p_silver, (1 - p_gold - p_silver)}
ここで，Yはいずれの事象が観測されたかを表す名義変数 (e.g., 0 = はずれ，1 = 銀のエンゼルが当たった，2 = 金のエンゼルが当たった)，p はそれぞれの事象が生起する確率を成分とするベクトルを表します (e.g., p = {1%, 5%, 94%}; 和は必ず100%になる)。いずれのエンゼルも当たらない確率は，金のエンゼルが当たる確率p_goldと銀のエンゼルが当たる確率p_silverを用いて1 - p_gold - p_silverと表現できます。要するに，ベルヌーイ分布を3値以上に拡張した分布がカテゴリカル分布です。もちろん事象の数が4つ以上の場合もカテゴリカル分布が利用できます (e.g., ダイスを振っていずれの目が出るか)。カテゴリカル分布に従う試行を独立にN回繰り返したときにそれぞれの事象が生起する回数は多項分布に従います。
Count = {Count_Gold, Count_Silver, Count_Miss}
p = {p_gold, p_silver, (1 - p_gold - p_silver)}
多項分布は二項分布を3値以上に拡張したものに相当します。なお，Stanでモデルを記述するときにはNは省略します (Countの成分の和から自動的に計算されるため)。
続いて，前回の記事と同様，「金のエンゼル2倍キャンペーン」について考えます。チョコボール誕生50周年を記念して，2017年11月〜2018年2月の期間，金のエンゼルが当たる確率が通常の2倍になるというキャンペーンが行われました。また，この期間は銀のエンゼルは出現しないことが明言されていました。つまり，金のエンゼル2倍キャンペーン中は，金のエンゼルが当たる確率が2倍になり，銀のエンゼルが当たる確率が0倍になる，と表現することが可能です。そこで先ほどのモデルを少し修正し，
Count = {Count_Gold, Count_Silver, Count_Miss}
p = {p_gold × W_Gold, p_silver × W_Silver, (残り)}
と表現します。通常時には (W_Gold, W_Silver) = (1,1)，金のエンゼル2倍キャンペーン時には (W_Gold, W_Silver) = (2,0) という重みを付与することで，キャンペーン時の確率を表現することができます。さらに，銀のエンゼルが金のエンゼルの何倍当たりやすいかを表す生成量ratioは，
と表現することができます。
最後に，金のエンゼルが当たる確率と銀のエンゼルが当たる確率の事前分布について考えましょう。ベルヌーイ分布や二項分布のパラメータpの事前分布としては前回の記事でも使用したベータ分布がよく用いられますが，カテゴリカル分布や多項分布のパラメータpの事前分布としてはディリクレ分布が用いられることが多いです。ディリクレ分布はベータ分布を多値に拡張したものと考えることができます。ただし，今回の記事ではディリクレ分布は使用せず，金のエンゼルが当たる確率と銀のエンゼルが当たる確率のそれぞれについて事前分布を考えることにします (これらの事前分布が決まれば，いずれのエンゼルも当たらない確率の事前分布も自動的に決まります)。
まず，金のエンゼルが当たる確率の事前分布に関しては，前回の記事と同様，「今まで100箱以上のチョコボールを確認したが，金のエンゼルは1度も当たったことがない」という筆者の経験を反映させて，(α, β) = (1, 101) のベータ分布を考えようと思いますが，今回はそれに加えて，「金のエンゼルが当たる確率は50%未満である」という常識的な情報も付与した切断ベータ分布を金のエンゼルが当たる確率の事前分布として仮定します。
この切断ベータ分布は，値が0.5以上となる確率が0になるようにベータ分布を修正したものです。もっとも，もともとBeta(1, 101) から0.5以上の値が得られる確率はもともとほぼ0%なので，Beta(1, 101)を事前分布として用いるのとほとんど変わりません。
続いて銀のエンゼルが当たる確率の事前分布ですが，こちらに関してはこれといった事前知識がありません。銀のエンゼルは筆者自身何度か当てたことがあり，おもちゃのカンヅメと交換してもらったこともあるのですが，どのくらい当たりやすいのかについては実際のところよく分かりません。ただ，常識的に考えて50%よりも高いということは無さそうですので，
という一様分布を事前分布として仮定することにします。ところで，金のエンゼル・銀のエンゼルともに確率の上限が50%となるように事前分布を設定しましたが，これにより，金のエンゼルが当たる確率と銀のエンゼルが当たる確率の和が100%を超えてしまうのを未然に防ぐことができます。確率の和が100%を超えたり0%を下回ったりすることを許容するモデルでは推定がうまくいかないことがあるので，このように何らかの工夫を施すことが必要になる場合があります。他の方法としては，各事象の生起確率の事前分布をディリクレ分布などを用いて同時に設定したり，softmax関数やmin関数・max関数などを使って確率の和が1になるように補正したりする方法などがあります。Stanであれば，simplex型で変数を宣言することで，成分の和が1になるベクトルとしてパラメータを宣言することもできます (今回のモデルではsimplex型ではうまくいかないのでvecter型を指定します)。
3. 使用するデータおよびStanとRのコード
前回の記事と同様，分析で使用するデータは以下の出典から拝借いたしました。
- チョコボール統計
※「第107回 チョコボール計測」までのデータを使用しました（番外編のデータは未使用）。
→【通常期間】金：0, 銀：13，総数：289
→【2倍キャンペーン期間】金：1，銀：(0)，総数：114 - テレビ番組「月曜から夜ふかし」（日本テレビ，2015年9月15日放送）
→【通常期間】金：1, 銀：42，総数：915 - Youtube動画「【検証】チョコボール1000個開封して金のエンゼル、銀のエンゼルの確率を調べてみた」
→【通常期間】金：2, 銀：40，総数：1000 - Youtube動画「金のエンゼル確率2倍はウソだった？チョコボール1000個開封して確かめた結果…」
→【2倍キャンペーン期間】金：0, 銀：(0)，総数：1000
これらの結果を次のようなデータとしてまとめました。
| DataID | Count_Gold | Count_Silver | Count_Miss | W_Gold | W_Silver |
| 1 | 0 | 13 | 276 | 1 | 1 |
| 2 | 1 | 0 | 113 | 2 | 0 |
| 3 | 1 | 42 | 872 | 1 | 1 |
| 4 | 2 | 40 | 958 | 1 | 1 |
| 5 | 0 | 0 | 1000 | 2 | 0 |
（2018/12/23追記：「月曜から夜ふかし」のデータは，あらかじめ開封する箱の数を決めていたのではなく，金のエンゼルが当たるまで箱を開封し続けた結果得られたデータです。このような場合には二項分布や多項分布を用いるのは適切ではなく，本当は幾何分布を用いるべきなのですが，今回は目をつむることにします。この辺りの事情を考慮したモデルはまた別の機会にお伝えします。）
Stanの仕様に合わせて，今回は試行の総数は使用せず，「金のエンゼルが当たった回数 (Count_Gold)」，「銀のエンゼルが当たった回数 (Count_Silver)」，「いずれのエンゼルも当たらなかった回数 (Count_Miss)」を基本データとして渡します。また，金のエンゼル2倍キャンペーン期間の確率の重みをW_GoldとW_Slverの2つの変数を使って与えておきます (キャンペーン期間中は金のエンゼルが当たる確率が2倍，銀のエンゼルが当たる確率が0倍になる)。
前節で説明したモデルを実装するためのStanコードは以下の通りです。
Stanコード01: 金のエンゼルが当たる確率と銀のエンゼルが当たる確率の同時推定
data {
int N; //データセットの数
int Count_Gold[N]; //金のエンゼルが当たった回数
int Count_Silver[N]; //銀のエンゼルが当たった回数
int Count_Miss[N]; //いずれのエンゼルも当たらなかった回数
int W_Gold[N]; //金のエンゼルが当たる確率の重み (1か2)
int W_Silver[N]; //銀のエンゼルが当たる確率の重み (1か0)
}
transformed data{
int Count[N,3];
Count[,1] = Count_Gold;
Count[,2] = Count_Silver;
Count[,3] = Count_Miss;
}
parameters{
real<lower=0,upper=0.5> p_gold; //金のエンゼルが当たる確率
real<lower=0,upper=0.5> p_silver; //銀のエンゼルが当たる確率
}
transformed parameters{
vector[3] p[N];
for (n in 1:N){
p[n,1] = p_gold * W_Gold[n];
p[n,2] = p_silver * W_Silver[n];
p[n,3] = 1 - (p[n,1] + p[n,2]);
}
}
model{
//モデル
for (n in 1:N){
Count[n] ~ multinomial(p[n]);
}
//事前分布
p_gold ~ beta(1,101);
}
generated quantities{
real ratio;
ratio = p_silver / p_gold; //銀のエンゼルは金のエンゼルの何倍当たりやすいか
}
金のエンゼルが当たる確率と銀のエンゼルが当たる確率はどちらも50%以下である，という事前分布の情報は，パラメータを宣言するときに「upper = 0.5」と上限を指定することで表現できます。
続いて，データを上記Stanコードに渡すためのRコードを示します。
Rコード01: チョコボールデータをStanに渡す
library(rstan)
d <- read.csv("チョコボールデータ.csv")
stanmodel <- stan_model(file="Stancode/chocoball_multinomial.stan")
data <- list(N = nrow(d), Count_Gold = d$Count_Gold, Count_Silver = d$Count_Silver,
Count_Miss = d$Count_Miss, W_Gold = d$W_Gold, W_Silver = d$W_Silver)
fit2 <- sampling(stanmodel, data=data, seed=123,
chains=4, iter=50500, warmup=500, thin=1)
分析にはrstan 2.16.2およびR 3.4.2を使用しました。長さ50,500のチェインを4本発生させ，ウォームアップ期間を500とし，得られた200,000個のMCMCサンプルを用いて事後分布を近似しました。Rhatは1.01未満で，目視からも十分に収束していることが確認できました。
4. 金と銀のエンゼルが当たる確率の同時推定の結果
それでは，金のエンゼルが当たる確率と銀のエンゼルが当たる確率の推定結果をお示しします。ついでに，いずれのエンゼルも当たらない確率の事後分布も示します。次の図は，各々の確率の事後分布です。MAPと99%最高密度区間 (HDI) も示しています。
図1. 金のエンゼルが当たる確率p_gold，銀のエンゼルが当たる確率p_silver，およびいずれのエンゼルも当たらない確率の事後分布。MAPおよび99%最高密度区間も記してある。
まず，金のエンゼルが当たる確率のMAPは0.086% (99% HDI[0.017%, 0.259%]，95% HDI [0.027%, 0.208%]) でした。前回の記事の推定結果と比べると僅かに値が小さいですが，これは推定上の誤差であると思われます (もしかしたら切断ベータ分布を事前分布に用いた影響もあるかもしれません)。とはいえ小数第2位までの値は一致しているので，前回の結果が再現できたと言っても良いでしょう。続いて，今回新たに得られた結果として，銀のエンゼルが当たる確率のMAPは4.290% (99% HDI[3.281%, 5.511%], 95% HDI[3.506%, 5.201%])，いずれのエンゼルも当たらない確率のMAPは95.60% (99% HDI[94.38%, 96.62%], 95% HDI[94.66%, 96.37%]) であると推定されました。つまり結論としては，99%の確率で，金のエンゼルが当たる確率は0.02%─0.26%，銀のエンゼルが当たる確率は3.28%─5.51%，いずれのエンゼルも当たらない確率は94.38%─96.62%の範囲にあると主張することができそうです。
5. 銀のエンゼルは金のエンゼルの何倍出やすいか
最後に，「銀のエンゼルは金のエンゼルの何倍当たりやすいか」という疑問に答えることを試みます。この値は生成量のratioの事後分布を見れば分かります。
図2. 銀のエンゼルが金のエンゼルの何倍当たりやすいかを表す倍率ratioの事後分布。MAPと99%最高密度区間も示してある。
この結果から，MAPで見ると，銀のエンゼルは金のエンゼルの32.6倍当たりやすいということが分かります。ただし，99%最高密度区間は11.0─158.6倍 (95% HDIは13.6─102.4倍) とかなり広い点には注意が必要です。ところで，おもちゃのカンヅメは金のエンゼル1枚もしくは銀のエンゼル5枚と交換できますが，金のエンゼル1枚の価値は本当に銀のエンゼル5枚分と同等なのでしょうか？これを検証するために，MCMCサンプルを使って「銀のエンゼルが当たる確率は金のエンゼルが当たる確率の5倍よりも大きい」という事象の生起確率を計算したところ100%となりました (2万個のMCMCサンプルは全て5よりも大きな値でした)。したがって，少なくとも珍しさという点において，金のエンゼルは銀のエンゼル5枚分以上の価値があると考えて間違いはなさそうです。より具体的に言えば，金のエンゼル1枚は銀のエンゼル約11─159枚分の価値があるということになります。金のエンゼルが当たっても，すぐにおもちゃのカンヅメと交換してしまうのはもったいないのかもしれません。
なお，前回の記事でも同じことが言えますが，確信区間の広さからも分かるように，エンゼルが当たる確率についてより正確な推定を行うにはまだまだデータが足りません。さらに多くのチョコボールデータを収集することができれば，同じモデルを使ってより精度の高い推定を行うことが可能となります。観測とデータ収集は科学の基本。データをどんどん蓄積していけば，エンゼルたちの微笑みの秘密にもっともっと近づけるはずです。筆者も貢献できるように頑張ります。ちなみにチョコボールデータを用いた分析のアイディアはもう少し残っているので (離散パラメータの推定に関するネタと，もう少し実用的なネタ)，また機会があれば報告します。
