【かゆい所に手が届く】EMアルゴリズム解説とPythonによるGMM(混合ガウス分布)への実装。インターン趣味カメラレビューアカデミック映画/ドラマお問い合わせ管理人激推しの海外旅行便利グッズ!アカデミック【かゆい所に手が届く】EMアルゴリズム解説とPythonによるGMM(混合ガウス分布)への実装。2019年12月12日EMアルゴリズムを実装しながら学びたい!数式ではなんとなく理解したけど実際どうやって使うの...?今回は,確率モデルの最尤解を求めるための繰り返し近似法であるEM(Expectation–Maximization)アルゴリズムの解説とPythonで実装する方法をお伝えしていこうと思います。本記事はpython実践講座シリーズの内容になります。その他の記事は,こちらの「Python入門講座/実践講座まとめ」をご覧ください。本記事の実装はPRML「パターン認識と機械学習<第9章>」に基づいています。演習問題は当サイトにて簡単に解答を載せていますので,参考にしていただければと思います。解説したページになります。...コーディングに関して未熟な部分がたくさんあると思いますので,もし何かお気づきの方は教えていただけると幸いです。誤りについてもご指摘していただけると非常に助かります。全体のコードEMの気持ち機械学習を学ばれている方であれば,EMアルゴリズムが一番最初に大きく立ちはだかる壁だとも言えます。何をしたいのか,そもそも何のための手法なのかが見えなくなってしまう場合が多いと思います。そこで,今回は実装の前に,簡単にEMアルゴリズムの気持ちをお伝えしてから,ザッと数学的な背景をおさらいして,最後に実装を載せていきたいと思います。一問一答形式でEMアルゴリズムに関してみていきたいと思います。疑問1:何のための手法なの?EMアルゴリズムは,確率モデルのパラメータに関する最尤解を求めるための手法です。噛み砕いていきましょう。確率モデルというのは,簡単に言えば「現象の裏側に何か適当な分布を仮定すること」です。適当な分布を仮定したところで,その分布の形状を決定するパラメータを定めなくては,現象を確率モデルで説明することはできません。そこで,得られた情報(現象)に尤も良くフィットするようなパラメータを求める操作が「最尤推定」なのです。ちなみに,こちらの記事(【初学者向き】ベイズ推論とは?事前分布や事後分布をド素人向けに分かりやすく解説してみます!)では,尤度,事前分布,事後分布等のベイズ推論に登場する用語を簡単に説明していますので,ぜひ参考にして見てください。詳しく説明すると,EMアルゴリズムは潜在変数をもつ確率モデルの最尤解を求める手法です。多くのモデルの場合,対数尤度関数を計算すると「log-sum」が出現してしまいます。そのため,潜在変数を含むモデルの最尤解を解析的に求めることは難しくなっているのです。EMアルゴリズムでは,この問題を「下界を最大化していく」ことで解決していきます。疑問2:EMってどういう意味?EMアルゴリズムのEMは「Expectation–Maximization」の頭文字を取ったものです。アルゴリズム」といったところでしょうか。実は,あとで説明するように,EMアルゴリズムは二つのステップを繰り返します。それぞれのステップは「Eステップ」「Mステップ」と呼ばれていて,それぞれ「期待値を取る」「偏微分をして最尤解を求める」という操作に相当しています。EMアルゴリズムのEMというのは,近似手法の各ステップの操作を表しているのです。疑問3:K-meansとの関わりは?EMアルゴリズムを学んでいると,必ずと言っていいほど出てくる話題がK-meansとの関連性です。実は,K-meansは混合ガウス分布に関するEMアルゴリズムの分散を0に近づけたときの極限として与えられます。K-meansは,あるデータ点がどこのクラスタに属するかを「ハード」に割り当てます。4番目のデータ点は2個目のクラスタに属する,といった具合にです。混合ガウス分布に基づくクラスタリングでは,$n$番目のデータが$k$番目のクラスタに属する割合を$\gamma_{nk}$で表します。,$\gamma_{42}=0.3$であれば,4番目のデータ点が2番目のクラスタに属する割合は0.3というようになります。このようにクラスタリングすることを,「ソフト」な割り当てを行う,と言います。ソフトからハードな割り当てにするためには,ガウス分布において分散(共分散行列の対角成分)を0に近づければOKということです。疑問5:下界ってなに?対数尤度関数を2つの項に分解したときの片っぽうです。ちなみに,もう一方はKLダイバージェンスです。KLダイバージェンスが非負値なために,対数尤度関数が下回らない下限として「下界」と呼ばれています。Mステップでは,パラメータに関して下界を最大化します。疑問4:Q関数ってなに?一般化EMアルゴリズムの導出で出てくるのが,Q関数です。結論から言うと,Q関数は「尤度関数の潜在変数に関する期待値」のことです。何のことだかサッパリだと思います。数学的な背景は後ほど説明するとして,とりあえずの解釈をお伝えしておくと,下界をパラメータに関して整理しただけの式のことです。logの性質を利用して下界を和になおしたときに,パラメータに依存しない部分を除けばQ関数の出来上がりです。疑問5:一般化アルゴリズムのEステップでは事後確率を求めるのになぜ「E」ステップと呼ばれているの?クリティカルな疑問だと思います。後ほど説明するように,一般化EMアルゴリズムの説明では,Eステップでは潜在変数の事後確率を求めます。Eステップでは「期待値」を求めるはずでしたよね...?ここでは,私の見解をお伝えしようと思います。Eステップでは,潜在変数の事後確率を計算しますが,この理由は「Mステップの目的関数であるQ関数を計算するのに必要だから」です。数学的な背景としては,KLダイバージェンスを最小化する潜在変数の確率分布を求める操作に相当しています。Q関数自体が期待値を取る操作に相当していましたよね。Q関数値という「期待値」を計算するために必要な事後確率をEステップでは計算しているのです。注意が必要なのですが,混合ガウス分布のEステップでは事後確率を求めるのではなくて,負担率(事後確率に関する潜在変数の期待値)を計算します。実はこの計算は,事後確率を求める操作と本質的には同じです。混合ガウス分布のEステップで負担率を求める理由は,Mステップで求める最尤解を計算するのに必要だからです。疑問6:K(クラスタの個数)はどうやって決めるの?EMアルゴリズムでは,クラスタリングするクラスの個数は自分で決めなくてはなりません。ここが,問題点とも言えるでしょう。ちなみに,EMアルゴリズムをベイズ的に拡張した「変分ベイズ法」では,この問題点をエレガントに解決してくれます。疑問7:EMの解は必ず正しいの?実は,EMアルゴリズムで導出される最尤解は,局所解の可能性があります。そのため,初期値をうまく与えてあげる必要があります。初期値としてクラスタの平均値は標準正規分布から生成し,共分散行列は単位行列として設定するとよいでしょう。疑問8:潜在変数の事後確率って計算できるの?Eステップでは,潜在変数の事後確率を計算することになりますが,この計算って実際に行えるのでしょうか。第一に仮定している確率モデルを利用することができるからです。ガウス分布やポアソン分布を仮定していれば,式を知っていますから,計算することができますね。事後分布ですから,パラメータは既知であるという条件です。潜在変数以外の値は分かっている状態ですので,分布の式に当てはめることで事後確率を計算することができます。EMアルゴリズムまずは,一般のEMアルゴリズムについて確認していったうえで,混合ガウス分布に適用していくという流れにしましょう。以下では数式を用いていきます。p(\boldsymbol{X}|\boldsymbol{\theta})$でした。潜在変数$\boldsymbol{Z}$を噛ませることで,以下のように分解することができます。KL(q\|p)\end{eqnarray}第一項目は下界と呼ばれており,第二項目はKLダイバージェンスを表しています。それぞれの定義は以下の通りです。\}\end{eqnarray}ちなみに,式(1)の導出には2つの方法があります。1つは,対数尤度関数にイェンゼンの不等式を利用して下界を導出する方法です。これによれば,対数尤度関数は下界によって下からおさえられます。その差は$q$と$p$のKLダイバージェンスになります。2つ目の方法は,式(2)と式(3)を足すことで対数尤度関数になることを確認する方法です。これは,かなり天下り的な発想ですが,確認できることに変わりはありません。この2つの分解を利用して,対数尤度関数を最大化していきましょう。最適化するパラメータの対象は,$\theta$と$q$です。高校数学でも習ったように,$\theta$と$q$をそれぞれ固定しながらその和を最大化していくという方法をとりましょう。それでは,下界とKLダイバージェンスはどちらから注目していけばいいのでしょうか。それは,$q$のみに依存しているKLダイバージェンスです。一点注意点があります。感覚的には「KLダイベージェンスを最大化する$q$」を求めたくなりますが,実はこれが反対で,KLダイバージェンスを最小化する$q$を求めることになります。理由は,KLダイバージェンスを最小化することで,下界が最大化されるからです。このポイントを理解するためには,固定されている変数に注意する必要があります。p(\boldsymbol{X}|\boldsymbol{\theta})$も固定されます。和の結果が固定されているため,一方(KLダイバージェンス)が最小化されれば,もう一方(下界)は最大化されるということなのです。qに関する固定化こちらの考え方は簡単で,KLダイバージェンスが最小となるのは,KLダイバージェンスが0となるとき,つまり$q=p$のときです。これを数式で表してみましょう。\boldsymbol{\theta})\end{eqnarray}まずはKLダイバージェンスから注目しますこの操作が,EMアルゴリズムのEステップに相当します。,$\theta$を固定してKLダイバージェンスを最小化する$q$を求めるということです。このとき,固定化されている下界も$\theta$が変化するために変化します。続くMステップでは,下界をパラメータに関して最大化します。これは,対数尤度関数を底上げしているイメージで,実際に対数尤度関数は$q$には依存せず$\theta$に依存しているため,このMステップで増加していくことになります。今度は下界を最大化します実際に,数式で確認していきましょう。下界の定義から和に分解します。$\theta$に関係しない項は定数項にしてしまいます。const\end{eqnarray}ここで,Eステップの結果を思い出してみましょう。上の式では,まだKLダイバージェンスを0にした結果が反映されていませんね。\boldsymbol{\theta})$を代入してみましょう。ここでは,更新されていない$\theta$は$\theta^{old}$と書くことにします。\boldsymbol{\theta^{old}})$を代入することになります。const\end{eqnarray}はい,出ました!これが噂のQ関数です。const\end{eqnarray}上でもお伝えした通り,Q関数はしっかりと「尤度関数の潜在変数に関する期待値」となっていますね。一般化EMアルゴリズムをまとめておきましょう。混合ガウス分布とEMアルゴリズムさて,本題に入っていきたいと思います。今回は,EMアルゴリズムを適用する具体例として,混合ガウス分布(GMM)のパラメータ推定を挙げたいと思います。混合ガウス分布にEMアルゴリズムを適用するって...?となると思います。具体的なイメージをしていきましょう。あるデータが与えられます。そのデータを,いくつかのクラスタに分類したいと思っているとします。そのときに,与えられたデータ点は「いくつかのガウス分布から生成されたものだ」と仮定するのが,混合ガウス分布を利用したクラスタリングです。EMアルゴリズムでは,いくつかの混合ガウス分布のパラメータ(平均,共分散行列,混合係数)を推定するという流れになります。与えられたデータの生成をよく裏付けられるようないくつかのガウス分布の形状を求めることが,今回のミッションです。Eステップからみていきましょう。本来,事後確率を求めるのがEステップでした。GMMの場合は,「事後確率に関する潜在変数の期待値」を求めるのがEステップに相当します。Eステップで期待値まで計算してしまえば,そのままMステップに利用できるからです。実際に,数式をみてみましょう。ベイズの定理より,潜在変数の事後確率は以下のように変形できます。\}^{z_{nj}}}\end{eqnarray}Eステップでは,この期待値を計算してしまうのでした。実際に,式にしてみましょう。\boldsymbol{\Sigma}_j)}\end{eqnarray}出てきました。この結果が,負担率と呼ばれるものです。数式では$\gamma(z_{nk})$と定義されます。\gamma(z_{nk})\end{eqnarray}この式の解釈を考えてみましょう。分母は,混合ガウス分布の各クラスタに関する和です。分子は,混合ガウス分布のある特定のクラスタ1つの値です。負担率とは,あるクラスタ$k$が,データ$\boldsymbol{x}_n$をどれだけよく説明しているかを,全体に占める割合で表したものなのです。この負担率を使えば,「どのデータ」が「どのクラスタ」に「どの程度」属しているのかを知ることができます。ちなみに,負担率は「$z_{nk}=1$となる事後確率」とも等価です。,$\boldsymbol{z}$は二値変数のベクトルであり,事後確率に関する潜在変数の期待値を取れば,$z_{nk}=1$となる事後確率しか残らないからです。1$ですので消えてしまい,期待値と事後確率は完全に等価になります。以下に数学的な証明を載せておきます。\gamma(z_{nk})\end{eqnarray}他にも,以下のようにして「$z_{nk}=1$の事後確率と期待値は等価」を示すことができます。\gamma(z_{nk})\end{eqnarray}さて,まとめていきます。上述の通り,Eステップでは「事後確率に関する潜在変数の期待値(=負担率)」を計算するのでした。負担率は,混合ガウス分布の式を使って,割合のような形で定義されていましたね。\boldsymbol{\Sigma}_j)}\end{eqnarray}さて,次はMステップです。Mステップでは,Q関数を最大化するパラメータを偏微分を利用して求めていきます。1$という制限がありますので,ラグランジュの未定乗数法を利用するということです。1$の制限をつけたラグランジュ関数を$G$とおくことにします。1)\end{eqnarray}それぞれ,実際に計算をしていきましょう。\frac{N_k}{N}\end{eqnarray}ただし,式変形の途中(※)で以下の公式を利用しました。X^{-T}\end{eqnarray}さて,まとめていきましょう。★GMM-EMのEステップ★Q関数を最大化するパラメータを偏微分,もしくはラグランジュの未定乗数法を利用して求める。\frac{N_k}{N}\end{eqnarray}以上で,理論的な背景は終了です。以下では,早速実装を確認していきたいと思います。最初に細かく関数を見ていって,最後に全体のコードをお伝えしたいと思います。本当であれば,しっかりとクラスを作って見通しの良いコーディングをするべきなのですが,ここでは関数をとりあえず定義していって,それらを組み合わせてEMという関数を作る方針でいきたいと思います。今回は以下のような用法を想定しています。全体のコードを「em.py」として保存し,コマンドライン上で以下のように引数を与えてあげれば,python上のコードでsysが拾ってくれます。負担率とパラメータはファイルに書き出して保存するようにしています。Axes3D標準的なライブラリ等をインポートしています。データを読み込むためにcsvもインポートしています。laは逆行列や行列式を求めるためにインポートしています。sys.argv[3]用法にも書いた通り,最初の引数はデータファイル名,次の引数は負担率をcsvで保存するためのファイル名,最後の引数はパラメータの値を保存するためのファイル名です。np.array(l)csvを読み込んで,numpyのarrayとして定義しておきます。今回は,10000個の3次元データを利用します。以下では「(サンプル数)×(次元数)」で保存されているとします。np.exp(exp)/denominここでは,ある1点のデータを多次元ガウス分布に突っ込んだときに出される出力を計算しています。ポイントなのですが,ここでは共分散行列の逆行列と行列式を引数として設定しています。この理由は,この関数内で逆行列を行列式を計算してしまうと,10000回も(同じ)その計算を繰り返すことになるからです。outputここでは,データ全体をガウス分布に突っ込んだ時の出力を行列で返しています。ここで逆行列を行列式を計算することにより,無駄な計算を省くことができます。0)[None,:]上で求めた多次元ガウス分布の出力を,$\pi_k$で重み付けして出力します。出力は,後のために次元を拡張しておきます。Pi平均の初期値は標準正規分布から,共分散の初期値は単位行列,混合係数の初期値はクラス数で等分した確率としています。np.sum(logs)収束判定のために,対数尤度関数を計算する関数を定義しています。対数尤度関数は,混合正規分布のクラス数に関する和の対数を取って,それらのサンプル数に関する和を計算します。1閾値を設定して,対数尤度関数の変化分がその閾値を下回ったときに繰り返しを終了するように設定しています。更新式は,上でお伝えした通りです。thr)print("Iteration:"+str(n_iter))print("log-likelihood:"+str(log_list))実際に,上で定義したEMアルゴリズムを実行しています。ここでは,何回目の繰り返しなのか,対数尤度関数の値はいくらなのかを表示するようにしています。閾値は,適当に0.01としました。np.allcloseを利用しても良いかもしれません。クラスタの個数は何回か試したところ4個がよさそうでしたので,あらかじめ定めてしまいました。局所解に陥らなければ,出力はおおよそ以下のようになります。Pi.ravel()])今回は,負担率とパラメータを出力するのでした。負担率はすでにまとまった配列として定義されているため,今回新しくパラメータをまとめて配列化しておきます。その際,ravel()で一次元化しておかなければ,ファイルとして出力する際にエラーを起こしてしまいます。0)負担率の大きさから,データ点が属するクラスタを判別しています。クラスタとは,どの正規分布から生成したものとみなせるかを表しているものです。azim=45)plt.show()上で求めたindexをもとに,色分けして10000個のデータをプロットしています。for文で回していますが,もっと他に良い方法があれば教えていただけると助かります。writer.writerows(params)最後に,負担率とパラメータをファイルに出力して終了です。●PRML「パターン認識と機械学習<第9章>」pythonあわせて読みたいアカデミック【これなら分かる!】変分ベイズ詳解&Python実装。最尤推定/MAP推定との比較まで。Beginaidアカデミック【サーベイまとめ】自己回帰モデルとは?係数の求め方まで出来る限り分かりやすく簡潔に説明します。コメントをキャンセルメールアドレスが公開されることはありません。付いている欄は必須項目です過去問解説販売中!zuka京都大学で機械学習を学んでいます。iBTの出題形式が大幅変更!新テストの注意点まとめ。情報学のおすすめ参考書HOMEアカデミックプログラミング【かゆい所に手が届く】EMアルゴリズム解説とPythonによるGMM(混合ガウス分布)への実装。