多変量混合正規分布

前回の多変量正規分布の続き。
「多変量混合正規分布」について見てみる。

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

混合パラメータ


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

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}

ここで、∑は「変量x」と「変量y」の分散共分散行列を表す。


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


混合比は各クラスタの確率密度を押し下げるように作用する。
任意のxとyの確率密度は、各クラスタでの確率密度の和になる。


具体的な値をあてはめてみる。
パラメータは下表のとおりとした。
f:id:ichou1:20181114210802p:plain

分布1
分布2
2変量混合正規分布確率密度関数(Rで実装)
# parameter
K1_x_Norm <- c(1,2)
K1_y_Norm <- c(3,4)
K2_x_Norm <- c(5,3)
K2_y_Norm <- c(7,5)
K1_cov <- 0.9
K2_cov <- -0.8
K1_w <- 0.3
K2_w <- 1 - K1_w

f_Bivar_GMM <- function(x, y){

  K1_x_mu <- K1_x_Norm[1]
  K1_y_mu <- K1_y_Norm[1]
  K2_x_mu <- K2_x_Norm[1]
  K2_y_mu <- K2_y_Norm[1]
  K1_Sigma <- matrix(c(K1_x_Norm[2], K1_cov, K1_cov, K1_y_Norm[2]), nrow=2, byrow=TRUE)
  K2_Sigma <- matrix(c(K2_x_Norm[2], K2_cov, K2_cov, K2_y_Norm[2]), nrow=2, byrow=TRUE)

  return ((K1_w * dmvnorm(matrix(c(x, y), ncol=2), mean=c(K1_x_mu, K1_y_mu), sigma=K1_Sigma)) + 
  (K2_w * dmvnorm(matrix(c(x, y), ncol=2), mean=c(K2_x_mu, K2_y_mu), sigma=K2_Sigma)) )

}


確率密度関数が描く形状を視覚的に確認してみる。
左図 : 2変量正規分布を2つ(それぞれ、上記の分布1と分布2に従う)
右図 : 2変量混合正規分布(混合分布の要素は上記の分布1と分布2)

f:id:ichou1:20181108214809p:plain

x軸が [ -7 : 13 ]、y軸が[ -5 : 15 ]の範囲で上から見た図
f:id:ichou1:20181108215940p:plain

この形状が正しいことを積分で確認してみる。
確率密度関数積分結果は確率を表すので、右図の積分結果は(ほぼ)「1」になるはずである。
(左図は多変量正規分布2つ分なので「2」)

「cubature」パッケージで重積分を求める関数が用意されているので、これを使う。
引数はベクトルで渡す必要があるので、ラッパー関数を用意する。

# "x" is vector
f_Bivar_GMM_for_integrate <- function(x){ f_Bivar_GMM(x[1], x[2])}
library(cubature)

x_range_lower <- -7
x_range_upper <- 13
y_range_lower <- -5
y_range_upper <- 15
adaptIntegrate(f_Bivar_GMM_for_integrate, lowerLimit=c(x_range_lower, y_range_lower), upperLimit=c(x_range_upper, y_range_upper))

結果。積分値が「0.9998677」

$integral
[1] 0.9998677

$error
[1] 9.838417e-06

$functionEvaluations
[1] 5355

$returnCode
[1] 0

積分範囲を広げてみると、

x_range_lower <- -100
x_range_upper <- 100
y_range_lower <- -100
y_range_upper <- 100
adaptIntegrate(f_Bivar_GMM_for_integrate, lowerLimit=c(x_range_lower, y_range_lower), upperLimit=c(x_range_upper, y_range_upper))

ほぼ「1」であることを確認。

$integral
[1] 1.000001

$error
[1] 9.929906e-06

$functionEvaluations
[1] 6613

$returnCode
[1] 0