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

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

変量の数は「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、２で、それぞれ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