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

多変量混合正規分布のパラメータ推定に入る前に、どのようなモデルを扱おうとしているのか整理しておく。

以下の視点で考えた時、

変量が「１」or「多」
「混合分布である」or「混合分布ではない」

以下の4パターンに分類される。

(1変量)正規分布
(1変量)混合正規分布
多変量正規分布
多変量混合正規分布

今回は上記3の「多変量正規分布」について見てみる。

変量の数が「2」、「変量x」と「変量y」で表されるとすると、確率密度関数は以下のとおり。

$\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}$

以下のサイトを参考にさせていただいた。
http://cse.naro.affrc.go.jp/minaka/R/R-binormal.html

具体的に値を求めてみる。

変量x : 平均「1」、分散「2」の正規分布に従う
変量y : 平均「3」、分散「4」の正規分布に従う
変量x と変量y には相関がない(共分散が0)

であるとすると、分散共分散行列は以下のとおりとなる。

$\begin{eqnarray} \varSigma_{xy} &=& \begin{pmatrix} 2 & 0 \\ 0 & 4 \end{pmatrix} \end{eqnarray}$

このデータはどんな形をしているか。
Rの「mvtnorm」パッケージで、多変量正規分布の乱数を作れる「rmvnorm」関数があるので、これを使う。

# 各変量の平均
x_mu <- 1
y_mu <- 3

# 各変量の分散と共分散
x_sigma2 <- 2 
y_sigma2 <- 4
xy_cov   <- 0

# 分散共分散行列
Sigma.zero <- matrix(c(x_sigma2, xy_cov, xy_cov, y_sigma2), nrow=2, byrow=TRUE)

# 乱数生成
N <- 10000
set.seed(123); data <- rmvnorm(n=N, mean=c(x_mu, y_mu), sigma=Sigma.zero)
          [,1]     [,2]
[1,] 0.2073677 2.539645
[2,] 3.2043464 3.141017
[3,] 1.1828405 6.430130
(snip)
 [9998,] -0.7927958106  2.2253545108
 [9999,]  0.9042444874  4.8101121211
[10000,]  3.6597604578  4.8591241077

チェック。
標本分散、標本共分散については自作関数を使用。
work-in-progress.hatenablog.com

# xの標本平均
mean(data[,1])
[1] 0.9998627

# yの標本平均
mean(data[,2])
[1] 2.977238

# xの標本分散（自作関数を使用）
variance(data[,1])
[1] 2.011234

# yの標本分散（自作関数を使用）
variance(data[,2])
[1] 3.978216

# xとyの標本共分散（自作関数を使用）
covariance(data[,1], data[,2])
[1] 0.002540545

# xとyの相関係数
> cor(data[,1], data[,2])
[1] 0.0008981552

2次元散布図
f:id:ichou1:20181104115856p:plain

3次元散布図
f:id:ichou1:20181104120150p:plain

確率密度の計算。

先ほどの式をRで実装してみる。

f_2var_norm <- function(x, y, xmean, ymean, Sigma){

  # x と y をセンタリングしたcolumn vector（縦ベクトル）
  A <- matrix(c((x-xmean), (y-ymean)), nrow=2, byrow=TRUE)

  # 指数関数のべき乗部分
  power <- -0.5 * (t(A) %*% solve(Sigma) %*% A)
  
  denominator <- 2 * pi * sqrt(det(Sigma))
  
  return ( 1 / denominator  * exp(power) )
}

x=0.5、y=2.5の確率密度を計算してみる。

x_mu <- 1
y_mu <- 3
f_2var_norm(0.5, 2.5, x_mu, y_mu, Sigma.zero)
           [,1]
[1,] 0.05123421

検算。
「mvtnorm」パッケージの「dmvnorm」関数の値。

x_mu <- 1
y_mu <- 3
dmvnorm(matrix(c(0.5, 2.5), ncol=2), mean=c(x_mu, y_mu), sigma=Sigma.zero)
[1] 0.05123421

一致していることを確認。

各変量（xとy）ごとの形状を見るため、横からの角度でプロットしてみる。
f:id:ichou1:20181106204016p:plain
左図は変量xに着目した分布、右図は変量yに着目した分布

比較用として、1変量の確率密度を出力してみる。
f:id:ichou1:20181106204226p:plain

1変量正規分布における、平均値での確率密度は

# N(1, 2)における、平均値(=1)での確率密度
> dnorm(1, mean=1, sd=sqrt(2))
[1] 0.2820948

# N(3, 4)における、平均値(=3)での確率密度
> dnorm(3, mean=3, sd=sqrt(4))
[1] 0.1994711

一方、2変量正規分布における、変量xと変量yの平均値での確率密度は

> x_mu <- 1
> y_mu <- 3
> dmvnorm(matrix(c(1, 3), ncol=2), mean=c(x_mu, y_mu), sigma=Sigma.zero)
[1] 0.05626977

これは、1変量正規分布における平均値での確率密度を掛け合わせたものと等しい。

> dnorm(1, mean=1, sd=sqrt(2)) * dnorm(3, mean=3, sd=sqrt(4))
[1] 0.05626977

これは、2変量に「相関がない（＝相関係数が0）」ケース。
「相関がある」ケースではどうか、次回に続く。