パラメータ推定(多変量混合正規分布)

前回の多変量混合正規分布の続き。
パラメータ推定に関するメモ。

変量の数は「2」、「変量x」と「変量y」で表されるとし、混合数をKとおくと、2変量混合正規分布確率密度関数

2変量混合正規分布確率密度関数


\begin{eqnarray}
p(x,y \ | \ \pi_k, \mu_{x_k}, \mu_{y_k}, \varSigma_{xy_k} ) &=& \displaystyle \sum_{ k = 1 }^K \pi_k \ \ f(x, y \ | \ \mu_{x_k}, \mu_{y_k}, \varSigma_{xy_k} ) 
\end{eqnarray}

ここで、


\begin{eqnarray}
f(x, y \ | \ \mu_x, \mu_y, \varSigma_{xy} ) &=& \frac{1}{2 \pi \sqrt {\vert \varSigma_{xy} \vert} } \exp \left( - \frac{1}{2} 
\left(
  \begin{array}{c}
    x - \mu_x \\
    y - \mu_y \\
  \end{array}
\right)^{ \mathrm{ T } }
 {\varSigma_{xy}}^{-1}
  \left(
    \begin{array}{c}
      x - \mu_x \\
      y - \mu_y \\
    \end{array}
  \right)
\right) 
\end{eqnarray}

\begin{eqnarray}
\displaystyle \sum_{ k = 1 }^K \pi_k &=& 1 \quad \left( 0 \leq \pi_k \leq 1 \right) 
\end{eqnarray}

\begin{eqnarray}
\varSigma_{xy} &=& \begin{pmatrix} \sigma_{xx} & \sigma_{xy} \\ \sigma_{yx} & \sigma_{yy} \end{pmatrix}
\end{eqnarray}


観測データがN個であるとき、対数尤度関数は

対数尤度関数


\begin{eqnarray}
l(\pi_k, \mu_{x_k}, \mu_{y_k}, \varSigma_{xy_k}) &=& \log L(\pi_k, \mu_{x_k}, \mu_{y_k}, \varSigma_{xy_k} ) \\
&=& \log \displaystyle \prod_{ n = 1 }^N p(x_n,y_n \ | \ \pi_k, \mu_{x_k}, \mu_{y_k}, \varSigma_{xy_k} ) \\
&=& \displaystyle \sum_{ n = 1 }^N \log p(x_n,y_n \ | \ \pi_k, \mu_{x_k}, \mu_{y_k}, \varSigma_{xy_k} ) \\
&=& \displaystyle \sum_{ n = 1 }^N \log \left( \displaystyle \sum_{ k = 1 }^K \pi_k \ f(x_n, y_n \ | \ \mu_x, \mu_y, \varSigma_{xy} ) \right)
\end{eqnarray}


この対数尤度関数を最大化するパラメータを求める。
1変量混合正規分布のパラメータ推定と同様、EMアルゴリズムを使う。

Eステップ(負担率を求める)は、



\begin{eqnarray}
\gamma(z_{nk}) &=& P(z_k =1 | \ x_n,y_n) \\ \\
&=& \frac{\pi_k N(x_n,y_n| \ \mu_{x_k}, \mu_{y_k}, \varSigma_{xy_k})}{\displaystyle \sum_{k=1}^K \pi_k N(x_n,y_n| \ \mu_{x_k}, \mu_{y_k}, \varSigma_{xy_k} )}
\end{eqnarray}

分母は、分布全体における、サンプルn番目の(x,y)の確率密度。
分子は、k番目の分布における、サンプルn番目の(x, y)の確率密度。

ここで、



0 \leq \gamma(z_{nk}) \leq 1 \\
\begin{eqnarray}
\displaystyle \sum_{k=1}^K \gamma(z_{nk}) &=& 1 \\
\displaystyle \sum_{n=1}^N \displaystyle \sum_{k=1}^K \gamma(z_{nk}) &=& N
\end{eqnarray}


Mステップ(平均、分散、混合比を求める)は、



\begin{eqnarray}
\mu_{x_{k \ new}} &=& \frac{\displaystyle \sum_{ n = 1 }^N \left( \ \gamma(z_{nk}) \ x_n \ \right) }{\displaystyle \sum_{ n = 1 }^N  \gamma(z_{nk}) }  \\
\mu_{y_{k \ new}} &=& \frac{\displaystyle \sum_{ n = 1 }^N \left( \ \gamma(z_{nk}) \ y_n \ \right) }{\displaystyle \sum_{ n = 1 }^N  \gamma(z_{nk}) }  \\
\varSigma_{xy_{k \ new}} &=& \frac{
\displaystyle \sum_{ n = 1 }^N \left( \ \gamma(z_{nk}) \
\left(
  \begin{array}{c}
    x_n - \mu_{x_k} \\
    y_n - \mu_{y_k} \\
  \end{array}
\right)
\left(
    \begin{array}{c}
      x_n - \mu_{x_k} \\
      y_n - \mu_{y_k} \\
    \end{array}
\right)^{ \mathrm{ T } } \
\right)
}{\displaystyle \sum_{ n = 1 }^N  \gamma(z_{nk}) }  
\\
\pi_{k \ new} &=& \frac{ \displaystyle \sum_{ n = 1 }^N \gamma(z_{nk}) }{ \displaystyle \sum_{ k = 1 }^K \displaystyle \sum_{ n = 1 }^N \gamma(z_{nk}) }
\end{eqnarray}



サンプルを作って、パラメータ推定の過程を追ってみる。

サンプル数は3000個、rmvnorm関数を使って乱数を生成する。

乱数生成時のパラメータは以下のとおりとした(seedは「123」を指定)
f:id:ichou1:20181114211157p:plain

クラスタ1、2で、それぞれ900個、2100個の乱数を生成して、行結合。

                 [,1]        [,2]
   [1,]  0.1601659285  2.39412269  # start of cluster1 data
   [2,]  3.1835310020  3.55591811
   [3,]  1.6374620473  6.43394682
 (snip)
 [898,]  1.2690186956  3.34487192
 [899,]  1.1127917066  2.61540279
 [900,]  1.3495507885  1.07626005  # end of cluster1 data
 [901,]  4.0825508528  6.60102620  # start of cluster2 data
 [902,]  7.6669223995  6.84108593
 [903,]  4.8747766727 10.79301010
 (snip)
[2998,]  3.9620380727 10.26483654
[2999,]  5.1786302856  5.54080566
[3000,]  5.3288312182  4.77818261  # end of cluster2 data
散布図

f:id:ichou1:20181111075131p:plain


「mixtool」ライブラリの「mvnormalmixEM」関数を使えば、多変量混合正規分布のパラメータ最尤値を求められる。

検算の意味で最初に値を求めておく。
結果が返ってくるまで少々時間がかかっており、topコマンドを見てみるとCPUをフルに使っていた。

mixtool_val <- mvnormalmixEM(data)
number of iterations= 115 
最尤値での対数尤度
mixtool_val$all.loglik[length(mixtool_val$all.loglik)]
[1] -13616.8
最尤値(混合比)
mixtool_val$lambda
[1] 0.3008678 0.6991322
最尤値(平均)
mixtool_val$mu
[[1]]
[1] 1.047571 3.067607

[[2]]
[1] 5.023995 6.986341
最尤値(分散共分散行列)
# cluster 1
mixtool_val$sigma[[1]]
          [,1]      [,2]
[1,] 1.7966296 0.8329409
[2,] 0.8329409 4.1450229

# cluster 2
mixtool_val$sigma[[2]]
           [,1]       [,2]
[1,]  2.8968299 -0.9050588
[2,] -0.9050588  5.1021322

先ほどのEMアルゴリズムで最尤値を計算してみた結果。
初期パラメータは以下のとおりとした。
f:id:ichou1:20181114211821p:plain

対数尤度の遷移は以下のとおり。

対数尤度(初期パラメータ)

f:id:ichou1:20181114211804p:plain

対数尤度(EMステップ1回目終了後)

f:id:ichou1:20181114212001p:plain

対数尤度(EMステップ2回目終了後)

f:id:ichou1:20181114212933p:plain

対数尤度(EMステップ10回目終了後)

f:id:ichou1:20181114212516p:plain

対数尤度(EMステップ115回目終了後)

f:id:ichou1:20181114212205p:plain

mixtoolsパッケージ関数の値と一致していることを確認。