スクフェスで学ぶ確率・統計（ベイズ推定編）
今回はスクフェスのイベント「メドレーフェスティバル（以下メドフェス）」を題材にして，統計の基礎を学ぼうという企画です．ベイズ統計の対象としては，離散値を扱うので非常に単純ですが勉強にはちょうどいいでしょう．ちなみに私は統計に関しては初歩的な知識しかなく，ベイズ統計は今回初めて触りました．以下の説明でも間違っているところがあるかもしれない．
今回のお題
メドフェスで金報酬，銀報酬，銅報酬が出る確率p_g, p_s, p_cをそれぞれ見積もりたい．n番目の報酬E_nを収集した時点でのそれぞれの確率を推定せよ．
理論
今回のケースでは確率分布の形，即ち母数に対する推定を行いたいため，次に示す母数分布に対するベイズの定理を利用する．
今回の計算に利用する尤度関数f(E_n|θ)について，報酬は金(g)，銀(s)，銅(c)の３種であるから，3次の多項分布に従う．
事前に確率に関する情報を得られていないため，理由不十分の原則に従い，事前分布p(θ)は一様分布
に従うものとしたいところだが，これだと最尤推定法と変わらない結果になってしまう．上の定理に代入するとわかるが，事後分布と尤度関数が比例してしまう．何が問題かというと，最尤推定法では「金が全く出なかった場合」などの特殊な場合に対応できない（確率0となってしまう．この現象をゼロ頻度問題という）．そこでもっと常識的な推定をするため，事前分布と事後分布が変わらないような分布を仮定してあげることで解決する．これを自然共役分布といい，K=3の多項分布の場合には3次のDirichlet分布
となる．事後分布は尤度関数×事前分布であるから，
である．
（2015/9/7補足：読み返しているうちに「なんでガンマ関数の中に何の説明もなく多項分布のn_k!がねじ込まれてんだよおかしいだろ」と上式に疑問が生じましたが、よくよく考えてみたら重要なのは確率変数p_kの部分だけなので、あとのガンマ関数部分はDirichlet分布に無理矢理合わせるための単なる規格化定数と考えられます（多項分布のn_k部分が１段目に含まれていないのも同じ理由）。つまり上式の２段目は=ではなく∝で結ばれるべきです。）
このように自然共役分布を使うと事前分布と事後分布がパラメータのみ異なる同じ関数で表される．則ち，報酬が1個増えるごとにDirichlet分布の対応するパラメータがα_k→α_k +1とベイズ更新されることになる．たとえばE_i=Cだった場合，α_c→α_c +1という風に銅報酬に対応する分布のパラメータのみが更新される．なお，Dirichlet分布の主要な統計量は次の通り．
ちなみに，以上の話の結果だけを抽出すると，ゼロ頻度問題を解決するためのLaplaceスムージング
が，実は全ての初期値α_kが一定であるという仮定を入れたDirichlet分布を事前分布としたベイズ推定で算出した期待値と一致する，という結果が得られる．
やってみた
さて，実際に自分のログを元にExcelでベイズ推定してみる．最終結果だけを言うと，全132試行のうち，金が4回，銀が38回，銅が90回出現した．ただし自分のログでは，１回のメドフェスで２個以上貰える報酬も別々の試行としてカウントしており，友人サポートで貰える報酬も1試行としている．また当然ながら，金銀報酬UPのサポートは一切つけていない．
単純に割り算で確認すると，p_g=0.030, p_s=0.288, p_c=0.682である．まあこれでも悪い推定ではないのだが，せっかくなのでnが小さい時の途中途中の推定確率も出してみよう．
1. 事前分布パラメータの決定
パラメータが1より小さいとそのパラメータに付随する確率の最頻値がマイナス値をとってしまい，1だと最尤推定法と同様ゼロ頻度問題が生じるので，パラメータは1より大きくすべきである．また，金報酬は確率が非常に小さいこと，金銀銅の順に確率が大きくなることを考慮して，
α_g=1.05, α_g=2.35, α_g=3.65
とした．小数点以下は（両方の意味で）適当に決めた．
2. 報酬によるパラメータのベイズ更新
小難しい理論を並べたが，結局は上述したように，金が出たらα_g → α_g + 1，銀が出たらα_s → α_s + 1，銅が出たらα_c → α_c + 1とカウントしていくだけである．つまり今回のログでは，最終段階で
α_g=5.05, α_g=40.35, α_g=93.65
となる．
3. 各時点での統計量の計算
今回は事前分布，事後分布にDirichlet分布を仮定しているので，上述したDirichlet分布の統計量の公式に当てはめるだけで計算が終わる．則ち，最終段階での金報酬の平均値，最頻値，分散はそれぞれ，
である．銀，銅の場合も計算は同様．
4. グラフ化
事後確率の最頻値および平均値のベイズ推定値は以下のように収束していく．どうやら，金報酬のようにごく小さい確率では平均値が最頻値を上回る傾向にあり，金をとった瞬間に一気に上がってその後指数関数的に減少する，という動きになるようだ．当然ながら，最初の方は情報が少ないため推定値が荒れている．なお分散は，回数を重ねるたびに0に近づいていくという至極当たり前の結果が得られた．本当はルート取って標準偏差にすべきなんだけどだいたいの傾向を見るにはこれで十分．
番外編
さて先程は確率分布の統計量で推定したが，統計量だけだと全体像が見えづらいのでDirichlet分布を直接見てしまいたいという欲求に駆られる．幸いにして今回の調査対象は金・銀・銅の３種類だったので，全体の分布を直接２次元プロット可能である．「えっ３次元＋分布の４次元じゃんどうすんの？」って一瞬思ったかもしれないが，よくよく見ると確率の総和が１であるという拘束条件を使って３次元に落とし込める．今回特に知りたいのは金・銀の確率なので，これらを先頭行・列に持ってきて0~1まで自分の知りたい分解能で並べ，中のマス目をDirichlet分布の式
に従って計算する．ただしExcelで計算する場合は，オーバーフロー・アンダーフローが怖いので対数をとってから計算するのが無難．なお，計算結果は当然p_g + p_s ≧ 1の箇所とどちらかが0になる箇所は未定義となる（確率0の地点は対数を経由しなければ0として出る）．計算結果は以下のグラフの通り．等高線が高い場所ほどそれらしき場所である．
結論
金報酬が出る確率...約3%
銀報酬が出る確率...約30%
銅報酬が出る確率...約67%
が尤もらしいことが今回の結果から分かりました．次回は報酬UPに本当に差があるのかを統計的に明らかにしたいと思います．
参考文献
主に以下のサイトでベイズ推定について学びました．この場を借りて感謝します．
説明が丁寧で例題も多く，非常にわかりやすかったです．
従来の推定法とベイズ推定法の違い | Sunny side up!
ベイズ主義の考え方がよくわかりました．
歪んだサイコロでベイズ - Negative/Positive Thinking
ほぼ同じ問題だったので具体的な解き方について非常に参考になりました．
Bernoulli distribution and multinomial distribution (ベルヌーイ分布と多項分布)
Bernoulli分布と多項分布の関係がよくわかりました．
Dirichlet分布で感じた疑問に的確な答えを見せてくれました．
http://hosho.ees.hokudai.ac.jp/~kubo/ce/2013/kubo2013esj.pdf
統計的な態度がどうあるべきかが詳しく書いてあり，非常に勉強になりました．
（2015/9/7追記）
LDAでは何故ディリクレ分布を仮定するのか - SELECT * FROM life;
多項分布×Dirichlet分布の疑問が晴れました。
