ベイズマックスへようこそ。
このブログでは,心理学と統計が大好きなmutopsyが,
統計に関するあれやこれやを自由気ままに投稿します。

2018年12月22日

銀のエンゼルが当たる確率のベイズ推定:金のエンゼルの何倍当たりやすいか?

 こんにちは。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枚分)

【お詫び】当初の予定では,帰無仮説検定の代わりとしてのベイズファクターについて初学者向けの解説記事を書く予定でおりましたが,不測の事態が重なり準備が間に合わなかったため,急きょ内容を変更させて頂きました。楽しみにしていた方には申し訳ございません。いつかちゃんと書きます!

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つ以上存在する場合は,ベルヌーイ分布ではなく次のようにカテゴリカル分布を仮定するのが自然です。

Y ~ Categorical(p)
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 ~ Multinomial(N, p)
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 ~ Multinomial(N, p)
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は,

ratio = p_gold / p_silver

と表現することができます。

 最後に,金のエンゼルが当たる確率と銀のエンゼルが当たる確率の事前分布について考えましょう。ベルヌーイ分布や二項分布のパラメータpの事前分布としては前回の記事でも使用したベータ分布がよく用いられますが,カテゴリカル分布や多項分布のパラメータpの事前分布としてはディリクレ分布が用いられることが多いです。ディリクレ分布はベータ分布を多値に拡張したものと考えることができます。ただし,今回の記事ではディリクレ分布は使用せず,金のエンゼルが当たる確率と銀のエンゼルが当たる確率のそれぞれについて事前分布を考えることにします (これらの事前分布が決まれば,いずれのエンゼルも当たらない確率の事前分布も自動的に決まります)。

 まず,金のエンゼルが当たる確率の事前分布に関しては,前回の記事と同様,「今まで100箱以上のチョコボールを確認したが,金のエンゼルは1度も当たったことがない」という筆者の経験を反映させて,(α, β) = (1, 101) のベータ分布を考えようと思いますが,今回はそれに加えて,「金のエンゼルが当たる確率は50%未満である」という常識的な情報も付与した切断ベータ分布を金のエンゼルが当たる確率の事前分布として仮定します。

p_gold ~ Beta(1, 101)[0.0, 0.5]

この切断ベータ分布は,値が0.5以上となる確率が0になるようにベータ分布を修正したものです。もっとも,もともとBeta(1, 101) から0.5以上の値が得られる確率はもともとほぼ0%なので,Beta(1, 101)を事前分布として用いるのとほとんど変わりません。

 続いて銀のエンゼルが当たる確率の事前分布ですが,こちらに関してはこれといった事前知識がありません。銀のエンゼルは筆者自身何度か当てたことがあり,おもちゃのカンヅメと交換してもらったこともあるのですが,どのくらい当たりやすいのかについては実際のところよく分かりません。ただ,常識的に考えて50%よりも高いということは無さそうですので,

p_silver ~ Uniform(0.0, 0.5)

という一様分布を事前分布として仮定することにします。ところで,金のエンゼル・銀のエンゼルともに確率の上限が50%となるように事前分布を設定しましたが,これにより,金のエンゼルが当たる確率と銀のエンゼルが当たる確率の和が100%を超えてしまうのを未然に防ぐことができます。確率の和が100%を超えたり0%を下回ったりすることを許容するモデルでは推定がうまくいかないことがあるので,このように何らかの工夫を施すことが必要になる場合があります。他の方法としては,各事象の生起確率の事前分布をディリクレ分布などを用いて同時に設定したり,softmax関数やmin関数・max関数などを使って確率の和が1になるように補正したりする方法などがあります。Stanであれば,simplex型で変数を宣言することで,成分の和が1になるベクトルとしてパラメータを宣言することもできます (今回のモデルではsimplex型ではうまくいかないのでvecter型を指定します)。

3. 使用するデータおよびStanとRのコード

 前回の記事と同様,分析で使用するデータは以下の出典から拝借いたしました。

 これらの結果を次のようなデータとしてまとめました。

表1. 今回の分析対象となるチョコボールデータ
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) も示しています。

multinomial_result.png
図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の事後分布を見れば分かります。

goldangel_posterior.png
図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枚分の価値があるということになります。金のエンゼルが当たっても,すぐにおもちゃのカンヅメと交換してしまうのはもったいないのかもしれません。

 なお,前回の記事でも同じことが言えますが,確信区間の広さからも分かるように,エンゼルが当たる確率についてより正確な推定を行うにはまだまだデータが足りません。さらに多くのチョコボールデータを収集することができれば,同じモデルを使ってより精度の高い推定を行うことが可能となります。観測とデータ収集は科学の基本。データをどんどん蓄積していけば,エンゼルたちの微笑みの秘密にもっともっと近づけるはずです。筆者も貢献できるように頑張ります。ちなみにチョコボールデータを用いた分析のアイディアはもう少し残っているので (離散パラメータの推定に関するネタと,もう少し実用的なネタ),また機会があれば報告します。

posted by mutopsy at 00:00 | 統計

2018年12月18日

金のエンゼルが当たる確率のベイズ推定

 こんにちは。mutopsyです。この記事は,Stan Advent Calendar 2018の18日目のエントリー記事です。この記事では,森永製菓が発売しているおなじみのチョコレート菓子「チョコボール」で金のエンゼルが当たる確率の推定方法を例に,Stanによるベイズ推定のいろはをお伝えすることを目指します。チョコボール50周年を記念して行われた「金のエンゼル2倍キャンペーン」のデータを有効活用しつつ,事前分布としてベータ分布を用いた二項分布モデルで金のエンゼルが当たる確率を推定します。以前,ポン・デ・リングのネタ記事を書いたことがありますが,今度はチョコボールのガチ記事です。

2018/12/22追記続編にあたる記事を書きました。→「銀のエンゼルが当たる確率のベイズ推定:金のエンゼルの何倍当たりやすいか?

1. 目的と概要

 改めまして,こんにちは。先日ぶじに博士論文を提出することができたmutopsyです。面白いことに(?),博士論文を提出したその日から風邪を引いてしまい,まだ熱が引いておりませんが,昔から「風邪にはMCMC」「Stanは百薬の長」とも言いますので,パラメータを収束させつつ体の熱を発散させようという意気込みでこの記事を執筆することにいたします。

 タイトルの通り,この記事では「金のエンゼルが当たる確率」を主題として取り上げます。ご存知ない方はおられないかと思いますが,金のエンゼルとは,森永製菓から発売されているチョコレート菓子「チョコボール」の箱の上部にある「クチバシ」にごくまれに印刷されているもので,いわゆる「当たり」です。これ1枚と,夢がつまった「おもちゃのカンヅメ」を交換することができます。ちなみに銀のエンゼルも存在し,こちらは5枚でおもちゃのカンヅメと交換できますが,今回は銀のエンゼルについては触れないことにします(2018/12/22追記:続編で銀のエンゼルも扱いました。→「銀のエンゼルが当たる確率のベイズ推定:金のエンゼルの何倍当たりやすいか?」)。

 筆者はチョコボールの大ファンで,子どものころから今に至るまで何度となく食べ続けてきました。その甲斐あって,銀のエンゼルには何度かお目にかかり,おもちゃのカンヅメを入手することもできましたが,未だかつて金のエンゼルをこの目で見たことはありません。では一体,金のエンゼルが当たる確率は一体どのくらいなのでしょうか。ベイズ推定という武器を手に入れたいま,この謎に取り掛かろうと思います。

 とはいえもちろん,先行研究は多数存在します。例えばチョコボール統計さんは,日々データを蓄積し,さまざまな分析を通じてチョコボールの神秘を解き明かそうと試みておられます。また,後述するように,テレビやYoutubeの企画などでチョコボールを大量買いした結果を報告している先行研究も存在します。これらの先人の知恵をお借りし,巨人の肩に乗り,金のエンゼルが当たる確率をなるべく正確に推定するということが本記事の主たる目的です。

 基本的なモデルは非常にシンプルです。まず,あるチョコボールの箱のクチバシに金のエンゼルが存在する確率をpと置きます。このpが今回推定したいパラメータです。ここでは,確率pはチョコボールの味や開封者,購入時期の違いなどに依らず一定であると仮定します。さて,いま,目の前にチョコボールの箱があるとします。この箱のクチバシを確認して,金のエンゼルがあれば1,なければ0と記録することにします。なお,ここでは銀のエンゼルがあった場合も0として記録することにします。このとき,金のエンゼルがあったか否かを表す2値変数Yは確率pのベルヌーイ分布に従います。

Y ~ Bernoulli(p)

ベルヌーイ分布は非常に単純な分布で,確率pでY = 1,確率(1-p)でY = 0となることを表現しているだけに過ぎません。このような試行をCountTrial回繰り返したとき,金のエンゼルが当たる回数CountGoldは次の二項分布に従います。

CountGold ~ Binomial(CountTrial, p)

これは,表が出る確率がpであるコインをCountTrial回投げたときに表が出る回数の分布と同じです。以上より,「チョコボールの箱を開けた回数」と「その中で金のエンゼルが見つかった回数」のデータがあれば,金のエンゼルが当たる確率pを推定することができるということです。ところが,次節で紹介するような事情を考慮すると,もう少しモデルを柔軟にしたほうが良さそうであることが分かります。

2. 金のエンゼル2倍キャンペーン

 チョコボール誕生50周年を記念して,2017年11月〜2018年2月の期間,金のエンゼルが当たる確率が通常の2倍になるというキャンペーンが行われました。そこで,このキャンペーン期間中のデータを有効活用できるように,先ほどの二項分布モデルを次のように拡張します。

CountGold ~ Binomial(CountTrial, p × W)

変更箇所はWという重みを表す観測変数です。通常の商品であればW = 1,金のエンゼル2倍キャンペーン対象商品であればW = 2 として確率に重みを付けることにより,このキャンペーンの内容をモデルに反映させることができます。それでは,このモデルを使って金のエンゼルが当たる確率を推定してみましょう。

3. 使用させて頂くデータ

 分析で使用するデータは以下の出典から拝借いたしました。今回は使用しませんが,銀のエンゼルの集計結果も一緒に掲載しています。なお,金のエンゼル2倍キャンペーン期間中は銀のエンゼルが出ないようになっていたようです。

 これらの結果を次のようなデータとしてまとめました。

表1. 今回の分析対象となるチョコボールデータ
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 ~ Uniform(0, 1)
p ~ Beta(1, 1)

つまり,pは確率なので0─1 (i.e., 0%─100%) の範囲に収まるが,実際に確率がどの程度なのかは見当も付かない,という信念を反映させた事前分布となります。もちろんこれでも良いのですが,筆者としてはあまり腑に落ちません。上述した通り,筆者はこれまで幾度となくチョコボールを食べてきましたが,一度も金のエンゼルを見たことがありません。きちんと計測した訳ではありませんが,少なく見積もっても100箱以上のチョコボールは確認したはずです。したがって,金のエンゼルが当たる確率は相当低いに違いありません。少なくとも50%─100%という範囲に収まるという事はあり得なさそうですが,一様分布を仮定した場合,pが0%─50%の範囲に存在するのと同じ程度の確率で50%─100%の範囲に存在するということになってしまいます。筆者のこの信念を事前分布として反映させるために,

p ~ Beta(1, 101)

というベータ分布を事前分布として用いることにします。ベータ分布は2つのパラメータαとβによって形状が決まり,実現値は0─1の範囲をとります (0と1は含まない)。(α, β) = (1, 1) のとき,一様分布と一致します。2つのベータ分布を図示すると以下のようになります。

beta_distribution.png
図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の事後分布です (破線で事前分布も示しています)。

goldangel_posterior.png
図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を乗算すると,幾何分布と呼ばれる確率分布になります)。

goldangel_posterior.png
図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追記:最後の分析の結果に誤りがあったのでグラフと数値を訂正し,少し加筆しました。

posted by mutopsy at 05:00 | 統計