エンゼル出現確率の予測分布の推移を可視化

はじめに
やること
実装
- 結果
- コード
おわりに
参考
広告

はじめに

前回の投稿では、チョコボールの「金のエンゼル」の出現確率を、ベイズ推定（MCMCでの近似計算）により予測してみました。

今回はこの投稿の続きで、予測分布の推移を可視化してやろうと思います。

内容としては、pythonでgifアニメーションを作るという話です。（一本の記事にするにはネタが軽いですが。。。）

やること

分析カテゴリでの前回の投稿で、ベイズ推定を用いて金のエンゼルの出現確率を予測しました。

ベイズ推定というのは、超簡単に振り返ると、予測するパラメータを”分布（事後確率分布）”として予測するものでした。分布として予測するので、最尤推定と違って以下のようなメリットがあります。

予測値がどのくらいの確信を持っているのかを表現できる
信用区間を使って、予測値の下限を直感的に表現できる
稀な事象の確率を予測できる

特に３つ目が実用的には重要だと思っています。
稀な事象を予測したいという需要はかなり多いと思いますが、稀なのでデータを集めるのに非常に苦労することが多いと思います。そんな中、わからないなりに確率を予測できるので、何らかの対策を考えることもできるわけです。

ということで今回は、データが増えるに従って事後確率がどのように変化していくのかをGIFアニメーションで可視化してみます。

実装

pythonのグラフ描画ライブラリであるmatplotlibには、いくつかアニメーション作成用の関数があるようですが、今回はmatplotlib.animation.FuncAnimationを利用してGIFを作成してみます。詳しくは下記参考に記載のページを参照ください。

結果

まず最初に結果を載せます。

12/17現在、49箱の開封が済んでいます（12/15更新分）。その中で、エンゼルの出現は0です。

このデータでベイズ推定すると、95%信用区間は0.055となりました（下図）。前回、43箱開封時での95%信用区間は約7%だったので、データが増えたことによりエンゼル出現確率がもっと低いぞという確信を得たことになります。

f:id:hippy-hikky:20171218004136p:plain — 金のエンゼルの出現確率の予測分布。赤破線は95%信用区間を示し、値は0.055。

では、データが増える毎に、エンゼル出現確率の予測分布がどのように変化していくか見てみましょう。

f:id:hippy-hikky:20171218004400g:plain

x軸は固定しています。データが増える毎に95%信用区間を示す赤破線が0に寄っていく様子が見えます。

コード

python スクリプトで作成しました。

main以外は、前回のコードと全く同じです。

update関数を定義し、この中でデータ数を制限してMCMC予測をします。 update関数の中身も前回の予測のメイン部分をほぼそのままコピーしています。

84行目で、FuncAnimation関数にfigureオブジェクトとupdaterを渡し、アニメーションを作成します。 intervalで各グラフの表示間隔を指定しています（単位はミリ秒）。 framesでフレーム総数を指定します。今回はデータを一つづつ増やして予測分布を描くので、データ数分のフレームになります。

その他細かいところは、まだまだ勉強中です。もっと効率的な書き方あるかもしれません。知っている方、教えてください。