相変わらず、黄色い本の復習コーナーになってしまった。2.2。
二値変数の場合
前記事のサマリーです。
コイントスの結果は2値変数であった。
1回の試行の結果Xは、表{1}裏{0}のいずれかを取る。
{1}を取る確率μを用いて、
p(x|μ)=μ^x(1-μ)^{1-x}\equiv Bern(x|μ)
とした。(ベルヌーイ分布と呼ぶ。)
N回の試行の結果「N回のうちK回が表がでた」は、
p(k|N,μ)= {}_NC_k\ μ^k(1-μ)^{N-k} \equiv Bin(k|N,μ)
と表現されました。（二項分布と呼ぶ）
これを尤度関数に用いたベイズ推定で便利なように事前分布を、
\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}μ^{a-1}(1-μ)^{b-1} \equiv Beta(μ|a,b) \\
と置くと、
p(μ|D,a,b)\propto μ^k(1-μ)^{N-k} \times μ^{a-1}(1-μ)^{b-1}
となって、事後分布もBeta分布になって、めでたし。
コイントスの結果を与えるベルヌーイ分布・二項分布に対し、
コイン自体の性質を支配する事前分布を与えるのがβ分布。
多値変数
コインを拡張するとサイコロになります(断言)。
サイコロなので（暗に？）同時に2つの目が出ない、という制約を置いている。
また、必ず何かしらの目が出る、という制約も置いている。
1回の試行で1つの面が取りうる値は出る{1}出ない{0}の2値変数(これだけ見ればベルヌーイ分布)。
1回の試行でサイコロ(の全ての面が)取りうる状態は、例えば「5の目が出た」は、下記で表せる。
x=(0,0,0,0,1,0)^t
tは転置を表す(ベクトルは縦だよ、って事)。
つまり、全ての面の同時確率を考える必要がある。
K面体のサイコロを考え、$x_5$=1が生じる確率を$μ_5$として、$x_5=1$かつ他が全てゼロの確率は、
p(x|μ)=\prod_{k=1}^K μ_k^{x_k}
where
μ=(μ_1, μ_2, ..., μ_K)^t,\ μ_k≥0,\ \sum_k μ_k=1
（あとでここに戻ります）
ポイントは$x_k$乗になっているところ。
$x_k=0$となった面は、$μ_k$（$x_k=1$が生じる確率）がいくつであっても$μ_k$^x_kは1になる。
この様に互いに排他的な取りうる状態Kのうち1つを取るような変数(サイコロ)を多値変数と呼ぶ。
期待値
\mathbb{E}[p(x|μ)]=\sum_x x\prod_k μ_k^{x_k}=(μ_1, μ_2, ..., μ_K)^t=μ
正規化の確認
(いや、まぁ、サイは投げられるわけだから、何らかの面が出る、というのは前提なので、とりうる状態xについて確率を足し合わせたら絶対に1になるんだけれど)
\sum_x p(x|μ)=\sum_{k=1}^Kμ_k^{x_k}=1
確認ですが、Kは面ですよ。
xはK次元ベクトルです。
多値変数の最尤推定
コイントスと同様、N回の独立なサイコロ投げの集合Dを考えます。
D=\{x_1, x_2, ... x_n ..., x_N \}
Nは試行です。Kは面です。
ここで、このDに所属するxを説明できる最も尤もらしいμを推定したい。
p(D|μ)=\prod_{n} p(x_n|μ)
=\prod_{n} \prod_{k} μ^{x_{nk}}
=\prod_{k} μ^{\sum_{n} x_{nk}}
=\prod_{k}μ_k^{m_k}
where
m_k=\sum_n x_{nk}
x_nkとは、ある試行nにおいて面kが出たら{1}出なかったら{0}というヤツ。
それをnに対してsumするのでm_kは、N回の試行を通じて面kが出た数。
対数尤度をとって、
\begin{eqnarray}
\rm{ln}\big( p(D|μ) \big)
=\rm{ln}\big( \prod_{k}μ_k^{m_k} \big)
=\sum_k m_k\ \rm{ln}\ μ_k \\
\end{eqnarray}
定石通りμ_kについて微分して=0と置く。
\frac{m_k}{μ_{kML}}=0
あれ、μについての条件がもとまらない。
そこで、μについての制約条件=定数となる事を利用して、
関数=関数+λ\times (制約条件-定数)
と変換してからμについて微分すると、きっと係数λが残るからハッピー。
(ラグランジュ未定乗数法とかなんとか言う)
制約条件としては、(あとでここ)に書いていありますね。
\sum_k μ_k=1
よって、
\begin{eqnarray}
\rm{ln}\big( p(D|μ) \big)
=\sum_k m_k\ \rm{ln}\ μ_k +λ\big( \sum_k μ_k-1 \big)\\
\end{eqnarray}
μについて微分して=0とおくと、
\frac{m_k}{μ_{kML}}+λ=0 \Leftrightarrow μ_{kML}=-\frac{m_k}{λ}
制約条件より、
\sum_k μ_{kML}=1 \Leftrightarrow λ=-\sum_k m_k=-N
最後の一押しは、m_kは、N回の試行を通じて面kが出た数、だからです。
それをKについてsumしたら、試行回数と一致しますよね。
（1回の試行で必ず1つの面が出るため。）
従って、
μ_{kML}=-\frac{m_k}{λ}=\frac{m_k}{N}
これは、N回のうち面kが出た割合ですね。
前回、ベルヌーイ試行の集合に対して最尤推定をしたら平均が出てきましたね。
それと同じイメージです。
長かった…。
ベイズ化しよう。
以上が、「互いに独立な試行を説明するμ」の推定でした。
で、前回と同じ様に、事前分布という巨人の肩に乗る事にする。
つまり、「これまでの試行の結果を反映したμ」を推定したい。
このためには、最初に尤度関数としてp(D|μ)の形を考えよう。
それに合わせて共役事前な分布p(μ|D)をデザインすれば、きっと上手くいく。
多項分布
では、最初に尤度関数としてp(D|μ)の形を考えよう。
このDって何だろうか。
「N回の多値変数の試行」の結果を表現する確率分布が必要。
結果は「面1がm1回、面2がm2回、...、面KがmK回」として出てくるハズ。
例えば、
D=\{x_1, x_2, x_3, x_4, x_5\} \\
x_1=(1,0,0,0,0,0)\\
x_2=(1,0,0,0,0,0)\\
x_3=(0,1,0,0,0,0)\\
x_4=(0,0,1,0,0,0)\\
x_5=(0,0,1,0,0,0)\\
面1が出る{1}確率がμ1。これがN回の中でm1回起こる。
N回の中でm1回の取りうる組み合わせを考えて、
{}_NC_{m_1}μ_1^{m_1}={}_5C_2μ_1^{m_1}
面2が出る{1}確率がμ2。これが残りのN-m1回の中でm2回起こる。
N-m1回の中でm2回の取りうる組み合わせは
{}_{N-m_1}C_{m_2}μ_2^{m_2}={}_{5-2}C_1μ_2^{m_2}
以下同を全て掛け合わせるのだけれど、最初のコンビネーション部分は、
{}_{5}C_{2}\times{}_{(5-2)}C_{1}\times{}_{(5-2-1)}C_{2}\times{}_{(5-2-1)}C_{0}\times{}_{(5-2-1)}C_{0} \\
=\frac{5\times4}{2!}\times\frac{3}{1!}\times\frac{2\times1}{2!}\times1\times1
要するに、N席をサイコロK面で椅子取りゲームする組み合わせなので、分子は絶対にN!になる。
という事で、
\frac{N!}{m_1!m_2!...m_K!}\prod_k μ_k^{m_k}\equiv Mult(m_1, m_2, ..., m_K|μ,N)
多項分布(multinomial distribution)とモウします。
制約条件は、全ての試行で何かしらの面が必ず出る、すなわち、
\sum_k m_k=N
Nを既知として、コイツを尤度関数にすれば良さそうですな。
事前分布ドウする？
B系の呪文。
\begin{eqnarray}
事後分布
&\propto& 尤度関数 \times 事前分布\\
&\propto& \prod_k μ_k^{m_k} \times 事前分布\\
\end{eqnarray}
beta分布の時と同じ理屈でお行儀よく進めると、
事前分布\ p(μ|a)\propto \prod_k μ_k^{a_k-1}
だと共役事前分布になりますよね。
正規化項を付け加えて、
p(μ|a)=\frac{Γ(a_0)}{Γ(a_1)...Γ(a_K)}\prod_k μ_k^{a_k-1}\equiv Dir(μ|a)
where
a=(a_1, ...,a_K)^t,\ a_0 = \sum_k a_k\\
μ=(μ_1, μ_2, ..., μ_K)^t,\ 1≥μ_k≥0,\ \sum_k μ_k=1
ディリクレ分布(Dirichlet distribution)とモウします。
ふ〜。出てきおった。
事後分布ドウなる？
B系の呪文。
\begin{eqnarray}
事後分布\ p(μ|D,a)
&\propto& 尤度関数 \times 事前分布\\
&\propto& \prod_k μ_k^{m_k} \times \prod_k μ_k^{a_k-1}\\
&\propto& \prod_k μ_k^{m_k+a_k-1}\\
\end{eqnarray}
正規化項をキチンと計算すると、
\begin{eqnarray}
p(μ|D,a)
&=&\frac{Γ(N+a_0)}{Γ(m_1+a_1)...Γ(m_K+a_K)}\prod_k μ_k^{m_k+a_k-1} \\
&=&Dir(μ|m+a)
\end{eqnarray}
となり、やはりディリクレ分布。
期待値は、ごちょごちょ計算すると、
\mathbb{E}[p(μ|D,a)]=\frac{m+a}{\sum_k m_k+a_k}
に、なります。
結局ディリクレ分布とは何者か。
多項分布は、N回サイコロを投げてどの目がどれだけ出たか、を説明する確率分布でした。でもこの結果から、このサイコロを投げるとどの目がどれだけ出るのかを推定するためには、その事前分布が必要です。
そのために都合のイイようにデザインされているのがディレクレ分布です。
何が都合がいいかというと、共役事前分布になっている、という事。
二項分布とβ分布の関係にそっくりですね。
Appendixで、ごちょごちょ部分を詳述しても良いのですが、あらびき様よりも整頓して書く自信も無く、スカッと爽やかに割愛しようと思います。
(決してただただ面倒臭いとかTex組むのモウ疲れたという訳ではナイ)
ま、β分布がナントカなるんだから多値変数にしてもきっとナントカなるんでは、ぐらいの感触で許していただけませんでしょうかね。
いや、それより、この後の方が大変なんですよ。そっちに進みましょう。
この時は、まだ誰も、椅子取りゲームの概念が大事件を引き起こすとは思いもしなかったのだ… To be contonued..
参考文献
あらびき様: [ディリクレ分布まとめ]
(http://d.hatena.ne.jp/a_bicky/20100402/1270139105)
Wikipedia-En: [Dirichlet distribution]
(https://en.wikipedia.org/wiki/Dirichlet_distribution)
yamano357様： [トピックモデルことはじめ]
(https://speakerdeck.com/yamano357/tokyowebmining46th)
北大 鈴木様：PRML2.2解説