PRML 演習4.9(+Chapter4.2.2)

会社の同期でやっているPRML勉強会で演習問題4.9の解説とそれに伴い、4.2.2章の解説らしきものの担当になってQiita書いたので公開します。なおこのQiitaはFate/Stay Nightを見ながら作成されました。

##演習問題4.9 Introduction
演習問題4.9は次のような問題です。

-------
クラスの事前確率$p(C_k)=\pi_k$と一般的なクラスの条件付き確率密度$p(\phi|C_k)$によって定義されるKクラス分類問題の生成モデルを考える. ここで, $\phi$は入力特徴ベクトルである. 学習データ集合${\phi_n,t_n}$が与えられたと仮定する. ここで, $n=1,...,N$であり, $t_n$は, 1-of-$K$符号化法を使う長さ$K$の2値目的変数ベクトルである. つまり, パターン$n$のクラスが$C_k$である場合, 2値目的変数ベクトルは構成要素$t_{nj}=I_{jk}$を持つ. データがこのモデルから独立に抽出されると仮定すると, その事前確率に対する最尤解が以下の式で与えられることを示せ.
$$\pi_k=\frac{N_k}{N}$$
ここで, $N_k$はクラス$C_k$に割り当てられるデータの個数である.

-------

この問題がどういう問題か一言で言うと、多クラスに関する事前確率の最尤推定をしてくれ、と言う感じかなと。で、多クラスはいいけどそもそも二値分類の場合はどうすればいいのでしょうかと言う疑問が出てくるので、この問題を解く前にまずは該当する章を見てみますか、ということでChapter4.2.2をみてみます。

##Chapter4.2.2 最尤解
クラスの条件付き確率密度$p(X|C_k)$に対するパラメトリックな関数形を決めればクラスの事前確率$p(C_k)$とともにパラメータの値を最尤法を用いて決定できる、と前置きとともに2クラスの場合の話が展開されます。

データ集合${X_n,t_n}$が与えられている場合において、$n=1,...,N$、$t_n=1$はクラス$C_1$、$t_n=0$はクラス$C_2$を表すものとします。つまり、二値分類ですね。また、この時の事前確率を$p(C_1)=\pi$、$p(C_2)=1-\pi$と定義します。このような前提のもとで、クラス$C_1$からのデータ$X_n$に対して$t_n=1$、つまりクラス1になる確率は次のようになります。(正規分布を仮定)

$$
 p(x_n,C_1)=p(C_1)p(x_n|C_1)=\pi N(x_n|\mu_1,\Sigma)
$$

またクラス2に関しても同様に定式化できます。
$$
 p(x_n,C_2)=p(C_2)p(x_n|C_2)=(1-\pi) N(x_n|\mu_2,\Sigma)
$$

ではこれらを用いて尤度関数を定義しましょう。尤度関数はある確率$p$に関するパラメータによって構成される関数で、この関数を最大化することで「尤もらしい」パラメータの推定値を求められるわけですね。(PRML上巻p.22)

ここでは尤度関数はいくつかのパラメータを導入して次のように定式化することができます。($t,X$は与えられるデータ集合)
$$
p(t,X|\pi,\mu_1,\mu_2,\Sigma)=\prod_{n=1}^N{[\pi N(x_n|\mu_1,\Sigma)]^{t_N}}{[(1-\pi) N(x_n|\mu_2,\Sigma)]^{1-t_N}}
$$

ということで、この尤度関数を最適化すれば確率$p$を尤もらしくするパラメータがわかるので頑張りましょう。各パラメータについて最適化していきたいと思います。また、最適化は基本的に対数尤度関数について行います。

* $\pi$について

$\pi$に依存する対数尤度関数の項を抜き出してきます。
$$
\sum_{n=1}^N{t_n\ln\pi+(1-t_n)\ln(1-\pi)}
$$
これについて、$\pi$に関する微分=0として解いていきます。

```math
\begin{align}
\frac{\partial}{\partial \pi}{\sum_{n=1}^N{t_n\ln\pi+(1-t_n)\ln(1-\pi)}}=0 \\
\sum_{n=1}^N({t_n\frac{1}{\pi}-(1-t_n)\frac{1}{1-\pi}})=0 \\
\sum_{n=1}^N({(1-\pi)t_n}-\pi (1-t_n))=0 \\
\sum_{n=1}^N({t_n-\pi})=0 \\
\pi=\frac{1}{N}\sum_{n=1}^N{t_n}=\frac{N_1}{N}=\frac{N_1}{N_1+N_2}
\end{align}
```

ここで$N_1$はクラスC_1内のデータ総数、$N_2$はクラス$C_2$内のデータの総数を表しているので、$\pi$に関する最尤推定は、大方の予想通りクラス$C_1$内のデータの個数の割合となります。

さて、というわけでニクラスの場合における事前確率に関する最尤推定を無事に終えることができました。これを一般化することで演習4.9も解けるようになりますねー。([演習問題4.9 解答](##　演習問題4.9 解答))

この章では他の部分の項に関する最尤推定についても、ある程度ですが追っていきます。

* $\mu_1,\mu_2$に関する最尤推定

同じように関連する項を拾ってくると次のようになります。$N(x_n|\mu_1,\Sigma)$が正規分布であることに留意して式を展開します。

```math
\begin{align}
&\sum_{n=1}^N t_n\ln N(x_n|\mu_1,\Sigma) \\
&=\sum_{n=1}^N t_n \ln(\frac{1}{{2\pi}^{D/2}})(\frac{1}{{\Sigma}^{1/2}})exp\{-\frac{1}{2}(x_n-\mu_1)^T \Sigma^{-1} (x_n-\mu_1) \} \\
&=-\frac{1}{2}(x_n-\mu_1)^T \Sigma^{-1} (x_n-\mu_1)+(\mu_1に関係しない定数)
\end{align}
```

これについて、$\mu_1$に関する微分を先ほどの操作と同様に行い、式整理を行って次のようになります。(正規分布の最尤推定なので、p.91などをご参照ください)
$\mu_2$に関しても同様になるので一緒に書きます。

```math
\mu_1=\frac{1}{N_1} \sum_{n=1}^N t_n x_n \\
\mu_2=\frac{1}{N_2} \sum_{n=1}^N t_n x_n
```

* 共有共分散行列$\Sigma$に関する最尤推定

時間が足りなかったので割愛。参照すべきページは同様にp.91などです。

```math
S = \frac{N_1}{N}S_1 + \frac{N_2}{N}S_2 \\
S_1 = \frac{1}{N_1} \sum_{n\in C_1}(x_n-\mu_1)(x_n-\mu_1)^T \\
S_2 = \frac{1}{N_2} \sum_{n\in C_2}(x_n-\mu_2)(x_n-\mu_2)^T
```

ということで2クラスの場合における事前確率に関する最尤推定が完了しました。これを一般化に拡張していくこともそれほど難しくない、とPRMLではいっていますね...

### そもそもなぜ最尤推定をいまやっているのか？
ではなぜこの章で最尤推定をやっているのかという疑問ですが、4.2章が確率的生成モデルを扱う章だからですね。

>分類を確率的な視点から考え、線形決定境界を持つモデルがデータ分布に関する簡単な家庭を設けることでどのように生成されるかを示す. 1.5.4節では分類に対する識別的アプローチと生成的アプローチについて議論した. ここでは, クラスの条件付き確率密度$p(x|C_k)$とクラスの事前確率$p(C_k)$をモデル化するクラスの条件付き確率密度分布$p(x|C_k)$とクラスの事前確率$p(C_k)$からベイズ定理を使って, 事後確率$p(C_k|x)$を計算する.

これが4.2節の冒頭の文章ですが、ここにあるようにこの章では分類問題を確率的に考えようという章になっています。

そのあとロジスティックシグモイド関数・ソフトマックス関数などの話が続き、4.2.1節では連続値の入力がある場合に事前確率・事後確率をどのように定義することができるのかを確認しその挙動などをグラフなどで見ている、という感じですね。

これを受けて、4.2節では最尤推定を行いパラメータ最適化を行うという流れになっているのではないかと思われます。

さて、長々前置きをして来ましたが演習問題4.9の解答です。

##演習問題4.9 解答
まず多クラスの場合における尤度関数を定義します。2クラスの時の尤度関数を参考にしつつ、それを一般化すればよいですね。注意点としては、2クラスの時は二項分布みたいな形でかけますが一般化した際にはクラスが$k$存在するということで、$t_n,1-t_n$みたいに書けていた部分が、$t_{n1},t_{n2},...,t_{nk}$みたいになることにご注意ください。

いろいろ考えた末に次のように表現できると思います。

```math
\prod_{n=1}^N \prod_{k=1}^K [\pi_k N(x_n|\mu_k,\Sigma)]^{t_{nk}}
```

演習4.9で求められているのは事前確率に対する最尤解、つまり$\pi$についてなので、Chapter4.2.2の最尤推定で示した一番最初のことをするだけですね。

では手始めに対数にします。

```math
\begin{align}
&\ln \prod_{n=1}^N \prod_{k=1}^K [\pi_k N(x_n|\mu_k,\Sigma)]^{t_{nk}} \\
&= \sum_{n=1}^N \sum_{k=1}^K \ln[\pi_k N(x_n|\mu_k,\Sigma)]^{t_{nk}} \\
&= \sum_{n=1}^N \sum_{k=1}^K \ln[{\pi_k}^{t_{nk}} N(x_n|\mu_k,\Sigma)^{t_{nk}}] \\
&= \sum_{n=1}^N \sum_{k=1}^K (t_{nk}\ln\pi_k+t_{nk}\ln N(x_n|\mu_k,\Sigma))
\end{align}
```

$\pi$に依存する項に注目します。

```math
f(\pi_k) = \sum_{n=1}^N (\sum_{k=1}^K t_{nk} \ln\pi_k)
```

これについて、$\pi$に関して最適化していきます。ここではラグランジアンを使って解くことを考えます。

制約条件として$\sum_{k=1}^K \pi_k = 1$を使い、次のように解くべき問題を定義します。

```math
L(\pi_k,\lambda)=\frac{\partial}{\partial \pi}f(\pi_k)+\lambda(\sum_{k=1}^K \pi_k -1)
```

では$\pi$について微分します。

```math
\begin{align}
\sum_{n=1}^N (\sum_{k=1}^K t_{nk} \frac{1}{\pi_k})+\frac{\partial}{\partial \pi}(\lambda\sum_{k=1}^K \pi_k - \lambda)=0& \\
\sum_{n=1}^N \sum_{k=1}^K \frac{t_{nk}}{\pi_k}+\lambda=0& \\
N_k = \sum_{n=1}^N t_{nk} = -\lambda \pi_k& \\
\end{align}
```

ここで$\lambda$を明らかにしないといけないのですが、両辺をぐっと睨むとそもそもこいつらは合計したらなんとかなる奴らでした。$N_k$はkについて合計すると$N$になり、$\pi_k$はkについて合計すると確率なので1になります。
よって
$$
\lambda = -N
$$
です。
なので、これをさっきの式に戻して、$\pi_k$について整理すると
$$
\pi_k=\frac{N_k}{N}
$$

を得ることができます。ふー。