多変量混合正規分布 - ichou1のブログ

変量の数は「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

混合比は「0.3」
xは正規分布N(1, 2)、yは正規分布N(3, 4)に従う
xとyの共分散は「0.9」、相関係数は「0.32」

分布2

混合比は「0.7」
xは正規分布N(5, 3)、yは正規分布N(7, 5)に従う
xとyの共分散は「-0.8」、相関係数は「-0.21」

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