[Excel][VBA]EMアルゴリズムと混合正規分布でデータを分類する1

マサムネのイメージ図 機械学習

のっぴきならぬ事情により、python などで計算を行う事が出来ない人もいるかもしれません。そのような人の為に、エクセルだけで完結するコンテンツを作ります。
今回の記事では、混合正規分布によって、データを分類します。混合正規分布のパラメーターの推定には、EMアルゴリズムを使います。
今回の記事は、VBAでEstepまで実装します。

問題設定

データを混合正規分布で近似する目的は複数あります。

  1. データをクラスタ分けする
  2. データに含まれるクラスタの比率を推定する
  3. クラスタ毎の平均と分散を推定する

上記3つが主な目的です。やりようによっては、クラスタの数Kも推定する事が出来ますが、今回の記事では、Kは人間が与える事にします。
この記事の手法を使うのは、以下のような状況です。

  • 複数の山を持つデータを取り合えずクラスタ分けしたい
  • データが持つそれぞれの山の平均値や分散が知りたい
  • 異常検知手法として使い、異常なデータがどんな振る舞いをするか調べたい

混合正規分布

混合正規分布については、以下の記事で詳しく説明しています。

EMアルゴリズム

どんな確率モデルを使うのかだけ解説しておきます。
混合正規分布では、データ\( X= \{x_1 , \cdots , x_N \} \)が次のような過程に沿って生成されていると考えます。

  1. データをクラスタ\(1,\cdots , K \)を割り当てる確率分布\( p( \pi ) \)があり、混合比率\( ( \pi _1 , \cdots , \pi _K ) \)が決まる。
  2. 各kに対して、パラメータ\( \theta _k = (\mu _k , \Sigma _k )\)が事前分布\( p ( \theta _k ) \)から生成される。1
  3. データ\(x_n \) に対応するクラスタの割り当て\(s_n \)が比率\( \pi \)から選ばれる。
  4. \( s_n \)によって選ばれたk が、\( x_n \)のクラスタ番号になる。また、データ\( x_n \)が\( p(x_n |\theta _k ) \)で生成される。

モデルの気持ちは上記の通りですが、実際に使う式とモデルの図を列挙しておきます。

混合正モデルの図
混合正モデルの図

\(x_n \)のクラスタkが決まっている時、正規分布から生成されていると考えるのが、混合正規分布モデルです。2
$$\begin{eqnarray}
p(x_n | \theta _k ) = \mathcal{N} (x_n | \mu _k , \Sigma _k )
\end{eqnarray}$$
クラスタは\(s_n \)によって与えられます。\( \pi \)をパラメータにしたカテゴリー分布を使います。
$$\begin{eqnarray}
p(s_n | \pi ) = { \rm Cat } (s_n |\pi ) =\prod _{i=1} ^{K} \pi _i ^{s_{n, i}}
\end{eqnarray}$$
ただし、\(s_n \)はK次元のベクトルで、一つの成分だけが1で、それ以外の成分が0になっています。また、\(\sum _{i=1}^{K}\pi _i =1 \)です。
\(s_n \)の情報がある時、\(x_n \)を生成する確率分布は以下のようになります。
$$\begin{eqnarray}
p(x_n | s_n , \Theta ) = \prod_{k=1}^{K} p(x_n | \theta _k ) ^{s_{n,k} }
\end{eqnarray}$$
本来は、\( \pi \)やパラメーター\( \Theta \)も確率分布から生成されていると考えますが、今回は簡単の為に、定数として扱い、最尤法で推定されるべきパラメーターと考えます。
データ\(x \)自身が従う確率分布は、今までの仮定を使って以下のように計算されます。
$$\begin{eqnarray}
p(x_n ) &=& \sum _{s_n } p(s_n ) p(x_n | s_n , \Theta ) \\
&=&\sum _{k=1} ^{K} \pi _{k} \mathcal{N}(x_n | \mu _k , \Sigma _k )
\end{eqnarray}$$
ただし、\( \sum _{s_n } \)は、取り得る\(s_n \)すべてについて和を取るという意味です。3
\(p(x_n ) \)達を尤度だと思って、最尤法によりパラメータ\( \mu _k , \Sigma _k , \pi _k \)を決定していきます。
$$\begin{eqnarray}
L= \prod_{n=1}^{N} p(x_n)
\end{eqnarray}$$
その為に、EMアルゴリズムと呼ばれる手法を使います。

EMアルゴリズム

尤度に対して、最尤法を使うと、補助変数\( \gamma (s_{n,k} ) \)を使って以下のように推定値が求められます。
$$\begin{eqnarray}
\mu _k &=& \frac{ \sum_{n} \gamma (s_{nk}) x_n }{N_k} \\
\Sigma _k &=&
\frac{1}{N_k} \sum_n \gamma (s_{nk} ) (x_n -\mu_k ) (x_n – \mu _k ) ^{T} \\
\pi _k &=& \frac{N_k }{N} \\
\gamma(s_{nk}) &=&
\frac{ \pi_k \mathcal{N} \left( x_n| \mu _k , \Sigma _k \right) }{ \sum_j \pi_j \mathcal{N} \left( x_n| \mu _j , \Sigma _j \right) } \\
N_k &=& \sum_n \gamma (s_{nk} )
\end{eqnarray}$$
\( \mu ,\Sigma , \pi \)は一応式の形が分かりますが、\( \gamma (s_{n,k } ) \)達に依存しているので、解けたわけではありません。EMアルゴリズムでは、次のように段階的にパラメーターを求めます。

[E ステップ]
\( \mu _k , \Sigma _k , \pi _k \)を適当に与え、\( \gamma (s_{n,k} ) , N_k \)を計算する。

[Mステップ]
\( \mu _k , \Sigma _k , \pi _k \) を計算する。

Eステップ、Mステップの順に繰り返し計算を行い4 、尤度が収束したら計算を終了するのがEMアルゴリズムです。

VBAで実装する

以上の計算を、VBAで実装します。
dataというシートにデータが有るという前提で、以下の手順で計算させます。諸々の計算は、result シートを使います。5

  1. データの形を取得する(データの数、次元)
  2. クラスタ数Kを聞くメッセージボックスを表示する
  3. 1,2の情報から,\( \mu _k , \Sigma _k , \pi _k \)を適当に与える。
  4. EMアルゴリズムに沿って計算する

今回は、4.の途中である、Estep の計算までを実装します。6

配列の値を初期化するのに、成分毎で値を設定したので、凄く見づらいコードになってしまいました。7
大事そうな場所だけ解説します。基本的に、定義通りに計算するだけなのですが、python のように配列を作って、配列の一部を取り出そうとしたりすることが(多分)出来ません。
その為に、平均、分散、元データの一部を取り出す関数を作っています。

また、多変量正規分布を得る関数が無いので、自力で計算させています。

mu_l, Sigma_l はクラスタ毎の正規分布の平均と分散です。エクセルについている行列積の関数や、行列の逆行列を求める関数を駆使して計算させています。

以下のようなデータを与えてマクロを実行できます。

データ
emアルゴリズムで使うデータ

Eステップまでのコードを実行すると、result シートに平均や分散、Γの値が記録されます。

resultシート
result シート

まとめ

  • 混合正規分布モデルの説明をした
  • EMアルゴリズムの説明をした
  • Eステップまでvbaで実装した
  1. 平均値と分散の事です。
  2. 上の過程4です。
  3. \(s_{n,1} =1 \)となる\(s_n \)、\(s_{n,2} =1 \)となる\(s_n \) 、という具合です。
  4. 2回目以降のEステップでは、直前のMステップで計算された \( \mu _k , \Sigma _k , \pi _k \) を使います。
  5. 使ったエクセルファイルはgithubにあげておきます。
  6. vbaらしくないコードだと思うので、簡潔に書けるところがあれば教えてください。
  7. もしかしたら、result シートに直接出力するようにするかもしれません。
タイトルとURLをコピーしました