EMアルゴリズムについて勉強したので実装してみた（実装は後日...）
どうも、IS19erのTeraponです。この記事はISer Advent Calender 2018の12日目として書かれたものです。昨日はiwannattoさんによる「OCaml print tips: %a, Format module - 嵐が過ぎたあとに」でした。来年の関数論理型実験の時に改めて読み直したいと思います。
さて、本日の記事では、機械学習の著名な本「Pattern Recognition and Machine Learning」(通称PRML）の9章に登場するEMアルゴリズムについて説明した後、多変量ガウス分布に対してEMアルゴリズムを適用するクラスをPythonで簡単に実装したいと思いましたが、思いの外長くなりそうだったのと、枠が埋まりそうになかったのもあって、とりあえず説明を今日やって、実装パートは19日に再掲します。（上の画像はPRMLの上巻です。下巻もあります。9章は下巻です！）
ここで少し宣伝です。今私は6人のチーム（19erのよっちゃんさんもいます）で上述のPRML（リンクを踏むとAmazonに飛びます。が、最近PRMLの原著が無料公開されましたね...原著でいい人はわざわざ高い本を買う必要もありません。）の輪読会を執り行っており、過ぎたる夏休みに上巻を読み終えました。来たる春休みには下巻を読む予定で今のところ夏の６名に+２名の8人で行う予定ですが、心強い仲間を募っています。数学（特に線形代数、解析回り、最適化あたりなど）強い人がきてくれるととても助かります。それ以外にも興味のある人がもしいましたら話を聞くだけでも構いませんので私に声をかけてください！！
では本題へ...
最初に機械学習について軽く触れてからEMアルゴリズムの話に移りたいと思います。
（今回は一般的なEMアルゴリズムではなく、あくまで混合ガウス分布に応用することを主なテーマとしたEMアルゴリズムについて書いていきます。一般のEMアルゴリズムもPRMLの中で紹介せれています。）
よくわからなかったら機械学習ってこんな感じなんやくらいの気持ちで読んでください笑
※ここより下、私が機械学習について無知であるが故に嘘をたくさん書くかもしれません。お強い方々は恐縮ではございますが誤りにお気づきになられたら優しくご指摘いただけると光栄でございます。
簡単に概要の説明
機械学習の一目的として、たくさんの変数（たくさんとは限らないが）からとある値（もちろん複数でも構わない）を予測したいという野望があります。そこでまずはこれらの変数と予測したい値との間にはとある関数関係があると考えます。一般的には以下のような流れで機械学習を行います。（日本語おかしい？）
現在得られているたくさんの変数と予測したい値のデータセットの一部を使って予測モデルを構築（こちらを教師データという）
残りのデータを用いてそのモデルの信頼性を評価（こちらは検証データという）
どうしてこういうことを行うのかというと、できたモデルが教師データに適合し過ぎたあまりに新たに得られたデータに対する汎用的な予測能力を失う（これを過学習という）を防ぐためです。
機械学習で最も重要なのは構築される予測モデルの形状の骨の部分を制御する（ほんまか？）重みパラメータというものと、そこからの散らばりなどの肉付きの良さなどの部分を制御する（ほんまか？）ハイパーパラメータを適切な値に合わせる（チューニングという）ことです。要するに予測モデルの形状をうまくデータの分布形状にフィットするように変えていくということです。
（ただし、形状ではなく複雑さからアプローチするノンパラメトリックな方法もあります。また、チューニング以外のタスクとして、変数の追加や削除あるいは変数の一部に変数変換を行ってそれを新たな変数として組み込むこと（深層学習ではこれらの変数を自動抽出してくれる）、予測モデルに使われている分布を変更することなどがあります。）
そして最適な値を求める方法として、機械学習には大きく分けて２つの考え方があります。
頻度主義・・・現在得られているデータから尤もらしさを表す尤度関数を作り、それを求めたいパラメータについて最大化するという姿勢をとる。これだと、得られているデータが変なデータ（偏りがあるデータ）だった時にその影響を強く受けてしまい、新たに得られるデータへの汎用的な予測能力を失ってしまうという問題をはらんでいる。
ベイズ主義・・・そもそもデータ分布ってこんな感じだよねという想定をあらかじめしておく。（これを事前分布という）そのあと、得られたデータの意見も汲みつつ（事前分布と尤度関数を掛け合わせることにより）改めてデータ分布を予測する（これを事後分布という）という姿勢をとる。事後分布を求める際の計算はベイズの定理に基づいており、ここからベイズ主義と呼ばれている。事前分布を仮定している分、頻度主義に比べて汎化性能が高くなる。
そしてPRMLは後者のベイズ主義の方を明らかに支持している書物であり、それを読んで勉強したが故に私もベイズ主義の方が基本的に強いという認識を持っています。（ほんまか？）
概要の説明はこの程度にしておきます。もっと知りたいのであればPRMLの上巻に今述べたようなことは大体書いてありますので是非是非。
混合ガウス分布の話
データが単峰性を持つときはガウス分布に色々な工夫を凝らしてそのデータを表現しようとすることができましたが、データが多峰性を持つときは一つのガウス分布ではどれだけ頑張っても適切に表現することはできません。そのため、ガウス分布の線型結合でそれらを表現しようという気持ちになります。これが混合ガウス分布です。
$$ p(\boldsymbol{x}) = \sum_{k=1}^{K}\pi_{k}\boldsymbol{N}(\boldsymbol{x}|\boldsymbol{\mu_{k}},\boldsymbol{\Sigma_{k}}) $$
ただし、
$$ 0\leq\pi_{k}\leq1 $$
$$ \sum_{k=1}^{K}\pi_{k}=1 $$
の２式を満たすように$\pi_{k}$をとって正規化を施します。$\boldsymbol{N}(\boldsymbol{x}|\boldsymbol{\mu_{k}},\boldsymbol{\Sigma_{k}})$はk番目のガウス分布を表します。$\boldsymbol{\mu_{k}}$は平均ベクトル、$\boldsymbol{\Sigma_{k}}$は共分散行列を表しています。N個に分割された全データ集合$\boldsymbol{X}=(\boldsymbol{x_{1}},...,\boldsymbol{x_{N}})$について考えると、
$$ p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma})=\prod_{n=1}^{N}p(\boldsymbol{x_{n}})=\prod_{n=1}^{N}\sum_{k=1}^{K}\pi_{k}\boldsymbol{N}(\boldsymbol{x_{n}}|\boldsymbol{\mu_{k}},\boldsymbol{\Sigma_{k}}) $$
となります。これが先ほど説明した尤度関数になります。尤度関数の最大化は対数をとった対数尤度関数の最大化と等価なため、一般的には尤度関数を以下の対数尤度関数にしてから最大化を試みます。(頻度主義的な解法について考えています）
$$ ln\,p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma})=\sum_{n=1}^{N}ln\,\Biggl(\sum_{k=1}^{K}\pi_{k}\boldsymbol{N}(\boldsymbol{x_{n}}|\boldsymbol{\mu_{k}},\boldsymbol{\Sigma_{k}})\Biggr) $$
なお、ここでは詳しい説明は省略しますが、この尤度関数には特異性というものが存在し、１つのガウス分布要素の平均があるデータ点と一致した時に分散が極めて小さくなると、その項の値が発散してしまい、対数尤度関数が発散してとんでもない過学習を引き起こす恐れがあります。なおここでは触れませんが、この問題は、PRMLの10章で紹介されるEMアルゴリズムの一般化として得られる変分推論法のベイズ主義的なアプローチにより解決されます。
さてこの対数尤度関数からは残念ながら閉な関数としての解析解が得られません（混合ガウス分布が$\boldsymbol{x}$以外に$\boldsymbol{\pi}$に関する潜在変数を持っていることに関係がある、潜在変数の存在確認はPRMLのなかで行われているので、詳しく知りたい人は読んでみてください）。そのため、潜在変数を持つモデルの最尤解（尤度関数を最大化するパラメータの値）を求めるための、繰り返し的な手法を用いたEMアルゴリズムというエレガントかつ強力な手法を用います。
EMアルゴリズムの話
さて、だらだら書いてきましたが、ようやく本題であるEMアルゴリズムの説明に入ります。（ほとんどPRML9.2.2節の要約ですが...）
- 平均ベクトル$\boldsymbol{\mu_{k}}$の値の導出
尤度関数を最大化するために、まず平均ベクトル$\boldsymbol{\mu_{k}}$について尤度関数を微分して０とおきます。
$$ \nabla_{\boldsymbol{\mu_{k}}}ln\,p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma})=0 $$
$$ \Leftrightarrow\sum_{n=1}^{N}\frac{\pi_{k}\boldsymbol{N}(\boldsymbol{x_{n}}|\boldsymbol{\mu_{k}},\boldsymbol{\Sigma_{k}})}{\sum_{j=1}^{N}\pi_{j}\boldsymbol{N}(\boldsymbol{x_{n}}|\boldsymbol{\mu_{j}},\boldsymbol{\Sigma_{j}})}\boldsymbol{\Sigma_{k}^{-1}}(\boldsymbol{x_{n}}-\boldsymbol{\mu_{k}})=0 $$
ここで左辺のシグマの中にある分数部分は事後確率となっています。これを負担率$\gamma(z_{nk})$とおきます（この$z_{nk}$が先ほど述べた潜在変数に当たります）。つまり、
$$ \gamma(z_{nk})=\frac{\pi_{k}\boldsymbol{N}(\boldsymbol{x_{n}}|\boldsymbol{\mu_{k}},\boldsymbol{\Sigma_{k}})}{\sum_{j=1}^{N}\pi_{j}\boldsymbol{N}(\boldsymbol{x_{n}}|\boldsymbol{\mu_{j}},\boldsymbol{\Sigma_{j}})} $$
ということです。先ほど述べた特異性がないという仮定のもとで上の式を$\boldsymbol{\mu_{k}}$について解くと、
$$ \boldsymbol{\mu_{k}}=\frac{\sum_{n=1}^{N}\gamma(z_{nk})\boldsymbol{x_{n}}}{\sum_{n=1}^{N}\gamma(z_{nk})} $$
となります。簡単のために分母を以下のように$N_{k}$とおきます。
$$ N_{k}=\sum_{n=1}^{N}\gamma(z_{nk}) $$
この$N_{k}$は、後で対数尤度関数を$\pi_{k}$について微分するとわかりますが、k番目のガウス要素が表現するデータ集合の部分集合（クラスターという）に属するデータ点の数を表現しています。
- 共分散行列$\boldsymbol{\Sigma_{k}}$の値の導出
さて次に共分散行列$\boldsymbol{\Sigma_{k}}$について対数尤度関数を微分して０と置くことで、共分散行列についても最大化することを考えましょう。導出と書きましたが、途中計算が長いので省略して結果だけ書きます（ごめんなさい）（Texつらい）。
$$ \nabla_{\boldsymbol{\Sigma_{k}}}ln\,p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma})=0 $$
$$ \Leftrightarrow\boldsymbol{\Sigma_{k}}=\frac{1}{N_{k}}\sum_{n=1}^{N}\gamma(z_{nk})(\boldsymbol{x_{n}}-\boldsymbol{\mu_{k}})(\boldsymbol{x_{n}}-\boldsymbol{\mu_{k}})^{T} $$
これを見ると、負担率$\gamma(z_{nk})$を重みパラメータとする、各々のガウス分布の共分散行列に対する重み付き平均として混合ガウス分布の共分散行列が与えられていることがわかります。
- 混合係数（といいます）$\pi_{k}$の値の導出
さて、最後に混合係数$\pi_{k}$についても対数尤度関数を最大化させましょう。混合係数には、
$$ \sum_{k=1}^{K}\pi_{k}=1 $$
という制約があったので、制限付き最大化問題に帰結します。そのため、ラグランジュの未定乗数法を用いなければいけません。よって、
$$ \nabla_{\pi_{k}}\Biggl(ln\,p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma})+\lambda\biggl(\sum_{k=1}^{K}\pi_{k}-1\biggr)\Biggr)=0 $$
とします。これを解いていくと、
$$ \sum_{n=1}^{N}\frac{\boldsymbol{N}(\boldsymbol{x_{n}}|\boldsymbol{\mu_{k}},\boldsymbol{\Sigma_{k}})}{\sum_{n=1}^{N}\gamma(z_{nk})\boldsymbol{N}(\boldsymbol{x_{n}}|\boldsymbol{\mu_{j}},\boldsymbol{\Sigma_{j}})}+\lambda=0 $$
これを解くと、再び途中計算の記述は省略しますが、
$$ \lambda=-N $$
$$ \pi_{k}=\frac{N_{k}}{N} $$
を得ました。
- 上での結果を使ってEMアルゴリズムを作る
よくみてみると、負担率の計算には平均ベクトル、共分散行列、混合係数ベクトルの値が必要であり、平均ベクトル、共分散行列、混合係数ベクトルの計算には負担率の値が必要です。この相補的（日本語おかしい？）な関係を利用すれば良いのです。
つまり、負担率の計算と、３つのパラメータの計算を交互に行えば、値はどんどん最適解に近づいていくことがわかります。前者のステップにはEステップ、後者のステップにはMステップという名前がついています。そしてそれらを繰り返すからEMアルゴリズムという名前なのですね。
最終的には対数尤度関数が最大値に収束していくことになるので、閾値を作って、変化量がそれ以下になったらアルゴリズムを止めれば良いことになります。では最後にアルゴリズムを簡単にまとめましょう。
EMアルゴリズムまとめ
- ①パラメータを初期化する。
$$ \boldsymbol{\mu_{k}}=\boldsymbol{\mu_{k0}} $$
$$ \boldsymbol{\Sigma_{k}}=\boldsymbol{\Sigma_{k0}} $$
$$ \pi_{k}=\pi_{k0} $$
- ②Eステップ：負担率を計算する。
$$ \gamma(z_{nk})=\frac{\pi_{k}\boldsymbol{N}(\boldsymbol{x_{n}}|\boldsymbol{\mu_{k}},\boldsymbol{\Sigma_{k}})}{\sum_{j=1}^{N}\pi_{j}\boldsymbol{N}(\boldsymbol{x_{n}}|\boldsymbol{\mu_{j}},\boldsymbol{\Sigma_{j}})} $$
- ③Mステップ：負担率を用いて３つのパラメータの値を上から順番に計算する。
$$ \boldsymbol{\mu_{k}}^{new}=\frac{1}{N_{k}}\sum_{n=1}^{N}\gamma(z_{nk})\boldsymbol{x_{n}} $$
$$ \boldsymbol{\Sigma_{k}}^{new}=\frac{1}{N_{k}}\sum_{n=1}^{N}\gamma(z_{nk})(\boldsymbol{x_{n}}-\boldsymbol{\mu_{k}^{new}})(\boldsymbol{x_{n}}-\boldsymbol{\mu_{k}^{new}})^{T} $$
$$ \pi_{k}^{new}=\frac{N_{k}}{N} $$
ただし、
$$ N_{k}=\sum_{n=1}^{N}\gamma(z_{nk}) $$
である。
- ④収束の判定（初回はここには来ずに②に戻る）
$$ ln\,p(\boldsymbol{X}|\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma})=\sum_{n=1}^{N}ln\,\Biggl(\sum_{k=1}^{K}\pi_{k}\boldsymbol{N}(\boldsymbol{x_{n}}|\boldsymbol{\mu_{k}},\boldsymbol{\Sigma_{k}})\Biggr) $$
を計算して、その変化量と閾値$\delta$との大小を判定。$\delta$より大きければ②に戻る。
終わりに
さて説明（前提知識がないと読みにくいような不親切な説明で申し訳ないです...）もようやく終わったので、19日の記事ではEMアルゴリズムのクラスをPythonで実装していきます。そのあと、有名なクラス分類用のオープンデータなどを適当に２変数でとって図式化して、ステップを重ねるにつれて収束していく様を観察したいと思います。
結局なんでこんな記事を書いたかというと誰かが興味を持って輪読会に参加してくれないかなあと思ったからなんですね。こんな感じのことをひたすらやります！！下のリンクを踏むとPRMLの原著をPDFとしてダウンロード出来るところへ飛ぶので、よかったら一度ちょっとペラペラしてみたくださいね。
明日はるんおじさんの記事です！
ここまで長々と読んでいただいた方、どうもありがとうございました。