ベイズ推定でエンゼルの出現確率を予測する

はじめに

チョコボールといえばエンゼル。
エンゼルがどのくらいの確率で入っているのか？というのは 全国のチョコボールファンが常に気になることです。 今回は、ついにこの永遠の謎に挑戦してみます。

とはいえ、結局まだデータがない状態なので、出現確率はX％ですということは言えませんのであしからず。

レアな現象の発生確率予測

事故や地震、工場での生産ライン異常など、これらは発生確率に差はあるものの一般的に稀な現象です。 稀な現象を予測するというニーズは常にあるのですが、 古典的な解析ではこれらの確率を予測するというのは難しい問題です。 稀なのでデータが少なく、予測確率には常に大きな誤差を含んでいる可能性があります。

チョコボールにおいても、エンゼルの出現は稀な現象です。 実際、このブログを初めてから４３箱開封してきましたが、まだエンゼルは現れていません。

頻度でエンゼルの出現確率を表現しようとすると、現在では0/43=0.0%$0/43=0.0\mathrm{\%}$となってしまい、 エンゼルは存在しないという予測になってしまいます。 エンゼルは本当はいるはずで、手元のデータではまだたまたま現れていないというだけなので、 0.0%という予測結果は望ましくありません。

そこで、ベイズ推定を使います。

なお、現在のキャンペーンでは銀のエンゼルは入っていないらしいので (この記事)、 本当に0%で正しいです。 なので、以下の記事では「金のエンゼル」の出現確率の予測を扱うことにします。

【トップに戻る】

最尤推定でのエンゼル出現確率予測

ベイズ推定での予測に入る前に、 まずは、古典的な方法である最尤推定を使って、エンゼルの出現確率の予測をしてみます。

以前の記事に書いたように、 最尤推定とは、尤度関数が最大になるパラメータを求めるというものです。

エンゼル出現の問題は、エンゼルが出現するかしないかの２通りの結果しかありません。 そこで、事象の発生モデルには二項分布を用います（下の式）。

f(x)=(nx)px(1−p)n−x${\displaystyle f(x)=\left(\begin{array}{c}n\\ x\end{array}\right)p^x(1-p{)}^{n-x}}$

ここで、nは試行回数（データ総数）、xは事象の発生回数を表しています。 pは事象の発生する確率を表しており、今回知りたいのは、この値です。

二項分布でモデル化すると尤度関数はそのまま二項分布になるので、 対数をとって対数尤度が0になるパラメータpが求めたい最尤推定量です。

logL(p)=log(nx)+xlogp+(n−x)log(1−p)=0${\displaystyle \log L(p)=\log \left(\begin{array}{c}n\\ x\end{array}\right)+x\log p+(n-x)\log (1-p)=0}$

この式を展開すると、最尤推定量はp=xn$p=\frac{x}{n}$ となり、頻度が最尤推定量ということになり、 0%つまり金のエンゼルは存在しないという推定結果になります。

そんなわけないですね。 金のエンゼルは存在するはずです！（見たことないけど）

【トップに戻る】

ベイズ推定でのエンゼル出現確率予測

前節では、最尤推定に依って金のエンゼルは存在しないという結論が導かれてしまいましたが、 ベイズ的にはどのような確率が予測されるのかを試してみましょう。

ベイズ推定の基本

ベイズ推定というのは、以前の記事に書いたように、 ベイズの公式に基づいて、観測データを入力しながらパラメータθ$\theta$の事後分布p(θ|X)$p(\theta |X)$を更新していく方法でした。

p(θ|X)=p(X|θ)p(θ)p(X)${\displaystyle p(\theta |X)=\frac{p(X|\theta )p(\theta )}{p(X)}}$

事後分布が、知りたいパラメータθ$\theta$（この場合はエンゼルの出現確率）の分布を表しているのでした。

モデル設計

まずは、ベイズ推定するためにデータ生成モデルの設計をします。

エンゼルの出現はするかしないかの2値だというのは前節で述べたとおりです。 ここで、エンゼルが出現する確率をp$p$とします。 チョコボールを買う毎に確率p$p$でエンゼルが存在するか決まっており、 この過程をシミュレートするために、データ生成過程にはベルヌーイ分布を仮定します。 ベルヌーイ分布というのは、N=1のときの二項分布のことです。

続いて、事前分布を考えます。
金のエンゼルの出現確率は、公式には公表されている数値は見つけられていません（公表されてない？）。 このような場合、事前知識が全く無いとして一様分布を仮定する（無情報事前分布）か、 過去の先人のデータを参考にして事前分布を設計する（経験ベイズ）かという方法がありますが、 まずは無情報事前分布を仮定してみます。

MCMCでの予測

作成したスクリプトは以下の通りです。

Sorry, something went wrong. Reload?

Sorry, we cannot display this file.

Sorry, this file is invalid so it cannot be displayed.

view raw estimate-engel-rate.ipynb hosted with ❤ by GitHub

4番目のセルで予測を行う関数(getMCMCResult)を作成しています。

    # 出現確率pの事前分布
    p = pm.Uniform('p', lower=0, upper=1)
    # 観測を結びつける
    obs = pm.Bernoulli('obs', p, value=data, observed=True)

この部分で、事前分布に一様分布を設定し、 データの生成過程としてベルヌーイ分布を設定しています。

    model = pm.Model([p, obs])
    mcmc = pm.MCMC(model)
    mcmc.sample(n_sample, n_burn)

続いて、モデルを作成してシミュレーションしています。

予測結果

これまでに得られている４３点のデータで予測してみた結果は以下のようになりました。

ヒストグラムがエンゼル出現確率p$p$の事後分布です(bin幅は約1%)。 赤い破線が95%信用区間（約7%）を示しています。
この結果から、これまでのサンプルでは、エンゼルの出現確率は0~1%にピークが出ています。 また、95%信用区間が約7%なので、7%以下であると強く確信していることがわかります。

次に、データが得られる毎に事後分布がどのように推移していくのか確認してみます。

左上から、観測データが2, 5, 10, 20, 30, 43個だったらどのような予測になるのかを示しています。 データが集まる毎に予測分布が０に収束していく様子が見えています。

【トップに戻る】

終わりに

エンゼルの出現確率予測という永遠の課題に挑戦してみました。
結果としては、現在のところ、出現確率は7%以下であると強く確信していることがわかりました。

このようなレアな現象を扱うにはまだまだデータが少なく、 そんなに面白い結果は得られていないですが、 頻度（最尤推定量）だけでは見えなかったものが見えているのがベイズ推定の面白いところです。

今後もデータを貯めていって、エンゼルの出現確率の予測精度を高めていきたいと思います。

補足ですが、 普通の問題では、ベイズ推定で事後分布が得られたら、 精度を高めるためにデータを集める運動をしたり、 現在予測される期待値に基づいて戦略を考えたりという動きをすると思います。 ですが、金のエンゼルを当てて「おもちゃのカンヅメ」を獲得するという課題については、 今すぐデータの取得をやめてメルカリで買うのが最良の選択です！
1,500円くらいで買えちゃいます。。。

【トップに戻る】

参考文献

定番の参考書。 ただ、ベイズ推定については載っていない。

上記の書籍の続き。 ベイズ推定など統計学の醍醐味はこっちに書いてある。

pymcの使い方というところに焦点を当てており、細かい理論は書いてないが、 pythonでMCMCしたいときにはとても参考になる。

【トップに戻る】

森永製菓 チョコボール<ピーナッツ> 28g×20箱

SIMERST 改良版　携帯タイプはかり ポケットデジタルスケール(秤)　0.001g-100g精密 業務用(プロ用)　デジタルスケール 電子天秤