こんにちは。mutopsyです。この記事は，Stan Advent Calendar 2018の18日目のエントリー記事です。この記事では，森永製菓が発売しているおなじみのチョコレート菓子「チョコボール」で金のエンゼルが当たる確率の推定方法を例に，Stanによるベイズ推定のいろはをお伝えすることを目指します。チョコボール50周年を記念して行われた「金のエンゼル2倍キャンペーン」のデータを有効活用しつつ，事前分布としてベータ分布を用いた二項分布モデルで金のエンゼルが当たる確率を推定します。以前，ポン・デ・リングのネタ記事を書いたことがありますが，今度はチョコボールのガチ記事です。
2018/12/22追記：続編にあたる記事を書きました。→「銀のエンゼルが当たる確率のベイズ推定：金のエンゼルの何倍当たりやすいか？」
2019/12/15追記：この記事の内容をさらに発展させた内容が書籍になりました。→「たのしいベイズモデリング2：事例で拓く研究のフロンティア」の第1章
1. 目的と概要
改めまして，こんにちは。先日ぶじに博士論文を提出することができたmutopsyです。面白いことに（？），博士論文を提出したその日から風邪を引いてしまい，まだ熱が引いておりませんが，昔から「風邪にはMCMC」「Stanは百薬の長」とも言いますので，パラメータを収束させつつ体の熱を発散させようという意気込みでこの記事を執筆することにいたします。
タイトルの通り，この記事では「金のエンゼルが当たる確率」を主題として取り上げます。ご存知ない方はおられないかと思いますが，金のエンゼルとは，森永製菓から発売されているチョコレート菓子「チョコボール」の箱の上部にある「クチバシ」にごくまれに印刷されているもので，いわゆる「当たり」です。これ1枚と，夢がつまった「おもちゃのカンヅメ」を交換することができます。ちなみに銀のエンゼルも存在し，こちらは5枚でおもちゃのカンヅメと交換できますが，今回は銀のエンゼルについては触れないことにします（2018/12/22追記：続編で銀のエンゼルも扱いました。→「銀のエンゼルが当たる確率のベイズ推定：金のエンゼルの何倍当たりやすいか？」）。
筆者はチョコボールの大ファンで，子どものころから今に至るまで何度となく食べ続けてきました。その甲斐あって，銀のエンゼルには何度かお目にかかり，おもちゃのカンヅメを入手することもできましたが，未だかつて金のエンゼルをこの目で見たことはありません。では一体，金のエンゼルが当たる確率は一体どのくらいなのでしょうか。ベイズ推定という武器を手に入れたいま，この謎に取り掛かろうと思います。
とはいえもちろん，先行研究は多数存在します。例えばチョコボール統計さんは，日々データを蓄積し，さまざまな分析を通じてチョコボールの神秘を解き明かそうと試みておられます。また，後述するように，テレビやYoutubeの企画などでチョコボールを大量買いした結果を報告している先行研究も存在します。これらの先人の知恵をお借りし，巨人の肩に乗り，金のエンゼルが当たる確率をなるべく正確に推定するということが本記事の主たる目的です。
基本的なモデルは非常にシンプルです。まず，あるチョコボールの箱のクチバシに金のエンゼルが存在する確率をpと置きます。このpが今回推定したいパラメータです。ここでは，確率pはチョコボールの味や開封者，購入時期の違いなどに依らず一定であると仮定します。さて，いま，目の前にチョコボールの箱があるとします。この箱のクチバシを確認して，金のエンゼルがあれば1，なければ0と記録することにします。なお，ここでは銀のエンゼルがあった場合も0として記録することにします。このとき，金のエンゼルがあったか否かを表す2値変数Yは確率pのベルヌーイ分布に従います。
ベルヌーイ分布は非常に単純な分布で，確率pでY = 1，確率(1-p)でY = 0となることを表現しているだけに過ぎません。このような試行をCountTrial回繰り返したとき，金のエンゼルが当たる回数CountGoldは次の二項分布に従います。
これは，表が出る確率がpであるコインをCountTrial回投げたときに表が出る回数の分布と同じです。以上より，「チョコボールの箱を開けた回数」と「その中で金のエンゼルが見つかった回数」のデータがあれば，金のエンゼルが当たる確率pを推定することができるということです。ところが，次節で紹介するような事情を考慮すると，もう少しモデルを柔軟にしたほうが良さそうであることが分かります。
2. 金のエンゼル2倍キャンペーン
チョコボール誕生50周年を記念して，2017年11月〜2018年2月の期間，金のエンゼルが当たる確率が通常の2倍になるというキャンペーンが行われました。そこで，このキャンペーン期間中のデータを有効活用できるように，先ほどの二項分布モデルを次のように拡張します。
変更箇所はWという重みを表す観測変数です。通常の商品であればW = 1，金のエンゼル2倍キャンペーン対象商品であればW = 2 として確率に重みを付けることにより，このキャンペーンの内容をモデルに反映させることができます。それでは，このモデルを使って金のエンゼルが当たる確率を推定してみましょう。
3. 使用させて頂くデータ
分析で使用するデータは以下の出典から拝借いたしました。今回は使用しませんが，銀のエンゼルの集計結果も一緒に掲載しています。なお，金のエンゼル2倍キャンペーン期間中は銀のエンゼルが出ないようになっていたようです。
- チョコボール統計
※「第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 | CountGold | CountSilver | CountTrial | W |
| 1 | 0 | 13 | 289 | 1 |
| 2 | 1 | 0 | 114 | 2 |
| 3 | 1 | 42 | 915 | 1 |
| 4 | 2 | 40 | 1000 | 1 |
| 5 | 0 | 0 | 1000 | 2 |
（2018/12/23追記：「月曜から夜ふかし」のデータは，あらかじめ開封する箱の数を決めていたのではなく，金のエンゼルが当たるまで箱を開封し続けた結果得られたデータです。このような場合には二項分布を用いるのは適切ではなく，本当は幾何分布を用いるべきなのですが，今回は目をつむることにします。この辺りの事情を考慮したモデルはまた別の機会にお伝えします。）
今回の分析ではこれらのデータを使用します。複数の研究で得られたデータを統合して分析するという意味ではいわゆるメタ分析に近いかもしれませんが，今回はデータ収集者などに由来するバイアスは存在しないと仮定しているので，モデル自体はとてもシンプルです。もしかしたら出版バイアスぐらいは考慮したほうが良いのかもしれません（e.g., 金のエンゼルが1つも出なかった場合，どこにも報告されずにその結果がお蔵入りになる，など）。
4. 事前分布
ベイズ推定を行う際は，あらかじめパラメータの事前分布を設定しておく必要があります。今回の場合であれば，金のエンゼルが当たる確率pの事前分布を決めなければなりません。一番単純な決め方は，次のように一様分布を仮定することでしょうか。
⇔ p ~ Beta(1, 1)
つまり，pは確率なので0─1 (i.e., 0%─100%) の範囲に収まるが，実際に確率がどの程度なのかは見当も付かない，という信念を反映させた事前分布となります。もちろんこれでも良いのですが，筆者としてはあまり腑に落ちません。上述した通り，筆者はこれまで幾度となくチョコボールを食べてきましたが，一度も金のエンゼルを見たことがありません。きちんと計測した訳ではありませんが，少なく見積もっても100箱以上のチョコボールは確認したはずです。したがって，金のエンゼルが当たる確率は相当低いに違いありません。少なくとも50%─100%という範囲に収まるという事はあり得なさそうですが，一様分布を仮定した場合，pが0%─50%の範囲に存在するのと同じ程度の確率で50%─100%の範囲に存在するということになってしまいます。筆者のこの信念を事前分布として反映させるために，
というベータ分布を事前分布として用いることにします。ベータ分布は2つのパラメータαとβによって形状が決まり，実現値は0─1の範囲をとります (0と1は含まない)。(α, β) = (1, 1) のとき，一様分布と一致します。2つのベータ分布を図示すると以下のようになります。
図1. ベータ分布の説明。
ちなみに今回使用するBeta(1, 101)は，pの事前分布を一様分布とした場合に，100箱のチョコボールを確認したがどのクチバシにも金のエンゼルが印刷されていなかったときのpの事後分布と一致します。これは筆者の実際の経験をうまく反映した分布であるといえます (実際のβはもう少し大きいかもしれない)。この情報は正確にカウントした結果ではなく，あくまでも筆者の信念を反映したものであるという意味で，データではなく事前分布として扱います。
なお，データを見てから「確率は低そうだ」と判断して事前分布を設定し，その上で同じデータを使って分析してはいけません（情報が重複してしまうため）。もし「100箱開封したが金のエンゼルが1つも当たらなかった」という経験をデータとして扱うのであれば，事前分布には一様分布を用いるのが妥当でしょう（この場合，経験を事前分布として扱った場合と同じ結果が得られます）。もしくは，「当たり」なのだから確率は十分低いはず，という常識に即した信念を，緩い事前分布の形で表現するのも良いでしょう。どの事前分布を選ぶかは分析の目的と分析者の信念次第ですね (例えば，科学的な研究であればより客観的かつ恣意的でない事前分布を採用するべきでしょう)。
5. StanとRのコード
これまでに説明したモデルを実装するためのStanコードは以下の通りです。
Stanコード01: 金のエンゼルが当たる確率の推定
data {
int N; //データセットの数
int CountGold[N]; //金のエンゼルが当たった数
int CountTrial[N]; //総数
int W[N]; //重み (1か2)
}
parameters{
real<lower=0,upper=1> p; //金のエンゼルが当たる確率
}
model {
//モデル
for (n in 1:N){
CountGold[n] ~ binomial(CountTrial[n],p * W[n]);
}
//事前分布
p ~ beta(1,101);
}
generated quantities{
real p_zero[1000];
for (i in 1:1000){
p_zero[i] = (1-p)^i; //i箱見て1箱も金のエンゼルが入っていない確率
}
}
最後のgenerated quantitiesブロックでは，i箱のチョコボールを確認して1つも金のエンゼルが入っていない確率p_zeroを計算しています。今回の場合はStan内で生成量を計算するよりも，要約統計量を用いて事後的に生成量を計算したほうが効率的ですが，念のため書き方を示しておきました。
続いて，データを上記Stanコードに渡すためのRコードを示します。
Rコード01: チョコボールデータをStanに渡す
library(rstan)
d <- read.csv("金のエンゼルデータ.csv")
stanmodel01 <- stan_model(file="Stancode/chocoball.stan")
data <- list(N = nrow(d), CountGold = d$CountGold,
CountTrial = d$CountTrial, W = d$W)
fit1 <- sampling(stanmodel01,data=data,seed=123,
chains=4,iter=10500,warmup=500,thin=1)
分析にはrstan 2.16.2およびR 3.4.2を使用しました。長さ10,500のチェインを4本発生させ，ウォームアップ期間を500とし，得られた40,000個のMCMCサンプルを用いて事後分布を近似しました。Rhatは1.01未満で，目視からも十分に収束していることが確認できました。
6. 分析結果
それでは分析結果をお示しします。次の図は，金のエンゼルが当たる確率pの事後分布です (破線で事前分布も示しています)。
図2. 金のエンゼルが当たる確率pの事後分布 (破線は事前分布)。MAPおよび99%最高密度区間も記してある。
この結果から，金のエンゼルが当たる確率の点推定値はおよそ0.091%であると結論づけることができそうです。10,000箱に9箱の割合です。なかなか厳しいですね。また，99%最高密度区間 (HDI) の下限は0.018%，上限は0.262%という結果でした (95%HDIは[0.027%, 0.209%])。つまり，どんなに当たる確率が低かったとしても0.018%未満であるということはなさそうですし，どんなに確率を高く見積もったとしても0.262%を上回ることはなさそうです。
最後に，チョコボールをn箱空けた時に一度も金のエンゼルが出ない確率の事後予測分布を示します。これは先ほどのStanコードのgenerated quantitiesブロックで計算した生成量に基づく分析です。開封する箱の数を横軸に取ると，グラフは次のようになります (ちなみに，縦軸の値に金のエンゼルが当たる確率pを乗算すると，幾何分布と呼ばれる確率分布になります)。
図3. 金のエンゼルが一度も当たらない確率のMAPと95%最高密度区間を，開封する箱の数ごとに示したグラフ。
例えば，100回チョコボールを開封して金のエンゼルが1つも出ない確率はMAPで91.3% (95%HDI[81.1%, 97.3%]) となります。200回でも83.1% (95%HDI[65.4%, 94.3%]) です。たまにチョコボールを食べている人がたまたま金のエンゼルを当てられる確率はそこまで高くなさそうですね。筆者がこれまで一度も金のエンゼルを当てたことがないのも不思議ではなさそうです（その情報を使って分析しているので若干トートロジックではありますが）。一方，今回の事後予測分布から，762回チョコボールを開封すれば，金のエンゼルが1つも出ない確率のMAPが50%を下回るようです。つまり，762箱のチョコボールを確認すれば，少なくとも1箱から金のエンゼルが見つかる確率が50%を超えるであろうことが見込まれます。また，3,291箱開封すれば，少なくとも1つは金のエンゼルが見つかる確率が95%を超え，5,059箱で99%を超えることも示されました。本気で金のエンゼルを狙うのであれば，3,000〜5,000箱程度のチョコボールを用意すればほぼ確実でしょう (ただし，確信区間の広さを考慮するとそれでも足りないかもしれません)。ちなみに，金のエンゼル2倍キャンペーン時に1,000箱のチョコボールを開封して金のエンゼルが1つも手に入らない確率はMAPで16.2%，95%最高密度区間はおよそ1.5%─58.6%の範囲となります。先行研究のYoutube動画ではキャンペーン期間中に1,000箱買っても1枚も金のエンゼルが見つからなかったことが報告されていますが，これは十分ありうることと言えるでしょう。
今回の分析において留意すべき点はいくつかありますが，そもそも金のエンゼルが観測される頻度が極めて少ないことは気に留めるべきところかもしれません。実際，合計3,318箱ものチョコボールを開封しても金のエンゼルはたったの4枚しか出ていません。これだけ起こりにくい現象の生起確率を高い精度で推定するためにはもっともっと多くのデータが必要でしょう。銀のエンゼルが当たる確率についても同様の分析を行うことが可能ですので，結果を比較してみるのも面白いかもしれません (事前分布は変更する必要があります)。本当は「おまけ」として，多項分布を使って金のエンゼルが当たる確率と銀のエンゼルが当たる確率を同時に推定した分析結果も掲載する予定でいたのですが，残念ながら間に合いませんでした。申し訳ございません。いつかリベンジします（2018/12/22追記：リベンジしました→「銀のエンゼルが当たる確率のベイズ推定：金のエンゼルの何倍当たりやすいか？」）。おもちゃのカンヅメに対する人類の挑戦はこれからも続く……。今回は自分のデータではなく人様のデータを拝借させて頂きましたが，筆者自身もチョコボールデータの収集に貢献できればと考えています。けどその前にまず風邪を治します！おやすみなさい！
※2018/12/18 16:55追記：最後の分析の結果に誤りがあったのでグラフと数値を訂正し，少し加筆しました。
