Qiitaに書いてどーすんだ説はありますが。
ディリクレ分布の理解を進めるにはβ分布をガツンと理解しておく事に限る、と思います。
その足跡記事です。
結果的に、黄色い本の2.1をネチネチと説明した、という感じの文章になってしまいました。
1回のコイントス
ここに、{表1, 裏0}がある確率分布に従って出るコインがある。
次にコイントスをした際に、表裏どちらが出るかで私の人生が決まる。
やばいっすね。
計算しましょう。
コインのどちらが出るか、という確率は、何かμに支配されているはずです。
\begin{eqnarray}
p(x=1|μ)&=&μ,&\ where\ 0≤μ≤1 \\
p(x=0|μ)&=&1-μ&
\end{eqnarray}
このμのもとでxがどんな値をとるかという確率は、
μ^x(1-μ)^{1-x} \equiv Bern(x|μ)
これは、独立した1回の試行でコインがとる値の確率分布です。
（これをベルヌーイ分布 Bernoulli distributionと呼びます。）
例えばμ=0.5の時、xはBern(x|μ)に従うので、
p(x = 1|μ = 0.5) = 0.5^1\times(1-0.5)^{1-1} =0.5\\
p(x = 0|μ = 0.5) = 0.5^0\times(1-0.5)^{1-0} =0.5
と計算できます。
最尤推定
ベルヌーイ分布を支配しているμ。
我々は、その正体を明らかにせねばなるまい。
シスの暗黒卿は、パルパティーン議員だったのだ。
手元のデータ
D \{ x_1, x_2, ...\ x_i\ ... , x_N \}
があるとして、μを最尤推定してみる。
Dを独立なベルヌーイ試行の集合とすると、
尤度は、μのもとでそれぞれのデータが得られる確率を単純に掛ければよいから、
p(D|μ)=\prod_{i=1}^N p(x_i|μ)=\prod_{i=1}^N μ^{x_i}(1-μ)^{1-x_i}
尤度の最大化は対数尤度の最大化と等価なので、
\begin{eqnarray}
ln\ p(D|μ)
&=& ln\ \prod_{i=1}^N p(x_i|μ) \\
&=& \sum_{i=1}^N ln\ p(x_i|μ) \\
&=& \sum_{i=1}^N \big( x_i ln\ μ+(1-x_i)ln(1-μ) \big)
\end{eqnarray}
これが最大となるのは、微分=0の時なので、
\begin{eqnarray}
\frac{d}{dμ} \Big( ln\ p(D|μ) \Big) &=&
\frac{d}{dμ}\sum_{i=1}^N \big( x_i ln\ μ+(1-x_i)ln(1-μ) \big) \\
0&=&\frac{1}{μ_{ML}}\sum x_i\ - \frac{1}{1-μ_{ML}} \sum (1-x_i) \\
μ_{ML}&=& \frac{\sum x_i}{\sum x_i+\sum(1-x_i)} \\
&=& \frac{1}{N}\sum_{i=1}^N x_i
\end{eqnarray}
こ、これは…平均ですね(標本平均)。
つ、つまり、これまでのN回の試行Dにおいて、
p(x=1|μ)=μ=\bar{x}
これまで100回コインを投げて40回表が出たら、これまでの実績から考えた最も尤もらしいコインが表になる確率は0.4。
うーむ、そらそうだ。
事前分布
ただですね、これまで3回の実績しかないとします。
それが全て表だった場合、平均は1なので、以降永遠に表が出続けるというのが最も尤もらしいということになります。
これは3回で全てが表されているという前提が怪しい。
その限定された実績で過学習が起きているという事です。
今更ですが、コインは世界にこれ1枚では無いわけです。
もちろん、完全に同一なコインは1枚だけですが、大雑把に見れば、似たようなコインは山ほどあり、Homo族の猿の一種はすぐコイントスをする性質がありますよね。で、概ね50%50%ぐらいの割合なんじゃないっすかね？
これを前提として、今起きた「3回連続表」という現象を評価すうると、やや表に偏っているかな、ぐらいの評価になるはずです。これは「以降、永遠に表が出続ける」よりもずいぶんと穏やかで妥当な推測ですよね。
こうした前提の事を、「事前分布」と言います。
という事で、人生が掛かってますからね。
大先輩ケダリオンの如く、事前分布の肩の上に乗りましょう。
求めたいモノ
我々は既に決意を固めました。
もはや次の1回の(そしてこれまで行われてきたあらゆる)コイントスを、独立なベルヌーイ試行とは考えません。
次の1回のコイントスで表がでる確率は、これまでのデータに縛られる。
p(x=1|D,μ)
0 ≤ μ ≤ 1より、
\begin{eqnarray}
p(x=1|D,μ) &=& \int_{0}^{1} p(x=1|μ)\ p(μ|D)dμ \\
&=& \int_{0}^{1} μ\ p(μ|D)dμ
\end{eqnarray}
右辺の形は、「データに基づいてμが発生する確率」と「μの実現値」の積の総和ですわな。
これは、期待値と呼ばれるやつです。
p(x=1|D,μ)=\mathbb{E}[μ|D]
データが与えるμの期待値を求めればよろしい。
データとは、「N回の試行の結果、k回の{1}が得られた」という情報です。
（あとでここに戻ります）
N回試行の尤度
今、我々は、B系のモノの考え方を採用しています。
事後分布 \propto 尤度関数 \times 事前分布
推定したいのはμなので、事前分布・事後分布はp(μ|D)という形になるはずです。
従って、尤度関数は、p(D|μ)の形になる。
このDには、「N回投げたらk回の{1}が得られた」という情報が入るはずです。
kの分布を考えると、Nは変数になるので、
p(k|N,μ)= {}_NC_k\ μ^k(1-μ)^{N-k} \equiv Bin(k|N,μ)
これを二項分布(Binomial distribution)とモウします。
ここで確率の左側に来ているのが、コインを何度か投げたうちの{1}が出た回数、という事に注意してください。
さて、μに依存する要素だけで書き直すと、
\begin{eqnarray}
p(k,l|μ)&=& {}_NC_k\ μ^k(1-μ)^{l},\ where\ N=k+l \\
& \propto & μ^k(1-μ)^{l}
\end{eqnarray}
事前分布、どーする
再掲。
\begin{eqnarray}
事後分布 &\propto& 尤度関数 &\times 事前分布 \\
& \propto & μ^k(1-μ)^{l} &\times 事前分布
\end{eqnarray}
もう一度回り道をしますが、ベイズ的にやりたい事は、事前分布を元に事後分布を計算し、それを事前分布として事後分布を計算し、...という情報のアップデートを積み重ねてより正確なμの事後分布を知りたい、という事です(逐次計算)。
そのためには、正の実数a,bを用いて、
事前分布 \propto μ^{a}(1-μ)^{b}
という格好になっているとスゴく便利ですよね。
\begin{eqnarray}
事後分布 &\propto& 尤度関数 \times 事前分布 \\
& \propto & μ^k(1-μ)^{l} \times μ^{a}(1-μ)^{b}\\
& \propto & μ^{k+a}(1-μ)^{l+b}
\end{eqnarray}
となって、これはすぐに次の事前分布として利用できます。
正規化項を加えて、
p(μ|a,b)=\frac{\Gamma(a+b+2)}{\Gamma(a+1)\Gamma(b+1)}μ^{a}(1-μ)^{b}
(20160819. 7of9様よりご指摘をいただき修正しました。)
としますが、ちょっとお行儀を良く書き直してβ分布(beta distribution)と呼びます。
\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}μ^{a-1}(1-μ)^{b-1} \equiv Beta(μ|a,b) \\
ただし、
\Gamma(x)\equiv \int_0^\infty u^{x-1}e^{-μ}du
なんでワザワザaをa-1に書き直したのかというと、期待値をすっきり書けるからです。
\mathbb{E}[Beta(μ|a,b)]= \frac{a}{a+b}
Γ関数や期待値の詳細をAppendixにすっ飛ばして先に進みます。
事後分布は、
p(μ|k,l,a,b) \propto μ^{k+a-1}(1-μ)^{l+b-1}
正規化項をしっかり計算すると、
p(μ|k,l,a,b) = \frac{\Gamma(k+a+l+b)}{\Gamma(k+a)\Gamma(l+b)} μ^{k+a-1}(1-μ)^{l+b-1}
これ自体がBeta分布になっています。（共役事前と言います）
ので、
\mathbb{E}[Beta(μ|k,l,a,b)]= \frac{k+a}{k+a+l+b}
（あとでここに）を再掲
p(x=1|D,μ)=\mathbb{E}[μ|D]
試合終了。
p(x=1|D)=\frac{k+a}{k+a+l+b}
Beta分布とは何者だったのか。
祭りのあとに残った感慨は、結局、一体全体、Beta分布とは何か、という話です。
コインを1回振って出た結果の分布が、ベルヌーイ分布。
コインをN回振って表が出る回数の分布が、二項分布。
そして、コインそのものの性質を決めるのが、Beta分布です。
（これがだけが言いたかった）
変数a, bを様々に取る事で、完璧なコインからイカサマコインまで、あらゆるコインのμの事前分布＆事後分布を与える事ができる(ようにデザインされている)。
こうした、実際の試行に先立って与えられたコインの性質を規定している変数a, bの事を、ハイパーパラメータと呼んだりします。
Appendix
1. Γ関数
\Gamma(x)\equiv \int_0^\infty u^{x-1}e^{-μ}du
β分布の正規化項で暗闇を破る雷鳴の如く登場したこの関数。
復習：部分積分
\int f(x)g'(x) dx = f(x)g(x)-\int f'(x)g(x)dx
ですよな。
f(x)=u^{x-1}\\
g(x)=-e^{-u}
として、見くらべて頂きながら、
\begin{eqnarray}
\Gamma(x)&=& \int_0^\infty u^{x-1}e^{-μ}du \\
&=& \Big( -u^{x-1}e^{-u} \Big)_0^\infty + (x-1)\int_0^\infty u^{x-2}e^{-u}du \\
&=& (x-1)\int_0^\infty u^{x-2}e^{-u}du \\
&=& (x-1)\Gamma(x-1)
\end{eqnarray}
という事で、これの正体は一般化された階乗関数です。
Nが自然数であれば、
\Gamma(N+1) = N!
2. β分布の正規性
本文より、
事前分布 \propto μ^{a}(1-μ)^{b}
となって、いるとスゴく便利で、
正規化項を加えて、行儀良く書き直してβ分布(beta distribution)と呼びます。
\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}μ^{a-1}(1-μ)^{b-1} \equiv Beta(μ|a,b) \\
え〜と、ですね。正規化項がこの形なのは、立派な理由がありまして。
立派すぎて、書くのを躊躇われます。
\int_0^1μ^{a-1}(1-μ)^{b-1}dμ\ \equiv\ B(a,b)
と置きます。β関数と呼びます。
μ=sin^2θ
と置きます。
\begin{eqnarray}
B(a,b)&=&\int_0^1μ^{a-1}(1-μ)^{b-1}dμ\\
&=&\int_0^{\frac{\pi}{2}} \big(sin^2θ \big)^{a-1} \big( 1-sin^2θ \big)^{b-1}\ 2sinθcosθdθ \\
&=&2\int_0^{\frac{\pi}{2}} sin^{2a-1}θ\ cos^{2b-1}θ\ dθ \\
\end{eqnarray}
となりますね。
一方、Γ関数においてu=s^2と置いて、
\begin{eqnarray}
\Gamma(x) &=& \int_0^\infty u^{x-1}e^{-μ}du
&=&2\int_0^\infty s^{2x-2}e^{-2s^2}ds
\end{eqnarray}
さぁて、うまくいくかな…!?
\begin{eqnarray}
Γ(a)Γ(b)&=&
\Big( 2\int_0^\infty x^{2a-1}e^{-x^2}dx \Big)
\Big( 2\int_0^\infty y^{2b-1}e^{-y^2}dy \Big) \\
&=& 4\int_0^\infty\int_0^\infty x^{2a-1} y^{2b-1}e^{-(x^2+y^2)}dxdy \\
&=& 4\int_0^\infty\int_0^{\frac{\pi}{2}}
(r\ cosθ)^{2a-1} (r\ sinθ)^{2b-1}e^{-r^2}dθdr \\
&=& 4\int_0^\infty\int_0^{\frac{\pi}{2}}
r^{2(a+b)-2}e^{-r^2}cos^{2a-1}θ\ sin^{2b-1}θ\ dθdr \\
&=&
\Big( 2\int_0^\infty r^{2(a+b)-2}e^{-r^2} dr \Big)
\Big( 2\int_0^{\frac{\pi}{2}} cos^{2a-1}θ\ sin^{2b-1}θ\ dθ \Big) \\
&=&Γ(a+b)B(a,b) \\
\end{eqnarray}
ぶふー。極座標変換ですね。x=r cosθ, y=r sinθ
xもyも0~∞つまり正の値だから、θの積分範囲は0~pi/2になる。
で、要するに
B(a,b)=\frac{Γ(a)Γ(b)}{Γ(a+b)}
つまり、
\begin{eqnarray}
Beta(μ|a,b)&=&\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}μ^{a-1}(1-μ)^{b-1} \\
&=&\frac{1}{B(a,b)}μ^{a-1}(1-μ)^{b-1}
\end{eqnarray}
と書けますね。
そして、
\begin{eqnarray}
\int_0^1 Beta(μ|a,b)dμ
&=&\int_0^1 \frac{1}{B(a,b)}μ^{a-1}(1-μ)^{b-1}dμ \\
&=&\frac{1}{B(a,b)}\int_0^1 μ^{a-1}(1-μ)^{b-1}dμ \\
&=&\frac{B(a,b)}{B(a,b)}=1
\end{eqnarray}
おお、正規化されていました。
3. β分布の期待値
息切れしてきました。
Texのトレーニングはもう…。
まず、B(a-1,b)を考えます。
\begin{eqnarray}
B(a-1,b)&=&\frac{Γ(a-1)Γ(b)}{Γ(a+b-1)} \\
&=& \frac{Γ(a)}{a-1} \frac{Γ(b)}{Γ(a+b-1)} \\
&=& \frac{a+b-1}{a-1} \frac{Γ(a)Γ(b)}{Γ(a+b)} \\
&=& \frac{a+b-1}{a-1} B(a,b) \\
\end{eqnarray}
いざ、期待値は、実現値と確率の掛け算の総和っと叫んで突撃。
\begin{eqnarray}
\mathbb{E}[Beta(μ|a,b)] &=& \int_0^1 Beta(μ|a,b)\ μ\ dμ \\
&=&\int_0^1 \frac{μ^{a-1}(1-μ)^{b-1}}{B(a,b)}\ μ\ dμ \\
&=&\int_0^1 \frac{μ^{a}(1-μ)^{b-1}}{B(a,b)}\ dμ \\
\end{eqnarray}
ここで、a=α-1と置くと、
\begin{eqnarray}
\mathbb{E}[Beta(μ|a,b)]
&=&\int_0^1 \frac{μ^{α-1}(1-μ)^{b-1}}{B(α-1,b)}\ dμ \\
&=&\frac{α-1}{α+b-1}
\int_0^1 \frac{μ^{α-1}(1-μ)^{b-1}}{B(α,b)}\ dμ \\
\end{eqnarray}
正規化されているから、積分の中身＝１。
α=a+1に戻して、
\mathbb{E}[Beta(μ|a,b)]=\frac{a}{a+b}
ぶ、分散は…気が向いたら追加します。(永久に向かない可能性大)
参考文献
書籍
黄色い本
webサイト
http://bio-info.biz/statistics/distribution_beta_distribution.html
https://lecture.ecc.u-tokyo.ac.jp/~nkiyono/2006/miya-gamma.pdf