2018-11-17

確率分布の推定

R 音声認識

これまで、モデルを表現するための要素である「確率分布の形」（確率密度関数）が分かっているものとして考えてきた。

確率分布の形が分からないケースを考えてみる。
（現実世界で推定しようとする対象は、こっちの方が多いと思われる）

以前、多変量混合正規分布のパラメータ推定の際に使った3000個の2変量データがあるとする。
f:id:ichou1:20181117082202p:plain

以下のことが分かっているので、対数尤度関数を決められた。

混合分布であり、混合数は2
各分布は正規分布に従う

ここで、（混合分布ではない）多変量正規分布と仮定してたらどうなるか。

3000個の標本の基本統計量を求めると以下のとおり。

# 平均
mean(data[,1])  # 変量x
[1] 3.827617
mean(data[,2])  # 変量y
[1] 5.80732

# 標本分散
variance(data[,1])  # 変量x
[1] 5.8918
variance(data[,2])  # 変量y
[1] 8.044347

# 標本共分散
covariance(data[,1],data[,2])  # 変量xと変量y
[1] 2.895581

これを最尤パラメータとして、標本の対数尤度を求めてみると（前回作った自作関数を使用）、対数尤度は「-14009.41」

log_likelihood_Bivar_norm(c(3.827617, 5.80732, 5.8918, 8.044347, 2.895581), data)
          [,1]
[1,] -14009.41

一方、「2変量混合正規分布（混合数2）」としての標本の対数尤度は「13616.7956」
f:id:ichou1:20181117083137p:plain

標本からモデルを決めるときは、対数尤度が最大になるような分布を選べばよい。

分布を推定するには大きく３つの方法に分けられる模様。

パラメトリックモデルを用いる手法（モデルを関数型と仮定して、関数のパラメータを調整）
ノンパラメトリックモデルを用いる手法（モデルを関数型と仮定しない）
セミパラメトリックな手法（モデルを関数型と仮定、パラメータの数を系統的に増やせるようにする）

混合分布モデルや階層型ニューラルネットワークは、3つの手法に該当。

（参考サイト）
確率密度分布の推定

音声認識で言うと、出力確率を前者に当てはめたのが「GMM-HMM」、後者に当てはめたのが「DNN-HMM」

続いて、モデルが決まった後の話。
ある点についての確率密度を求めたいとする。

左図が混合分布モデル、右図は混合分布ではないモデル。
f:id:ichou1:20181117095457p:plain

2つのモデルに対して、任意の変量xとyを与えたときの確率密度。
f:id:ichou1:20181117110126p:plain

x=3、y=5で誤差が大きくなっている（混合分布の2つの山に挟まれる部分）
混合分布モデルでは確率密度が「0.0133」になる点が、混合分布ではないモデルにおいては「0.0238」と確率が高くなっており、想定とは違う方向に「尤もらしい」を押し上げるように作用している。

2018-11-17

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

音声認識 R

前回、多変量混合正規分布のパラメータ推定について書いたが、（混合分布ではない）多変量正規分布のパラメータ推定についても書いておく。

変量の数は「2」、「変量x」と「変量y」で表されるとすると、2変量正規分布の確率密度関数は

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

$\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} \varSigma_{xy} &=& \begin{pmatrix} \sigma_{xx} & \sigma_{xy} \\ \sigma_{yx} & \sigma_{yy} \end{pmatrix} \end{eqnarray}$

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

対数尤度関数

$\begin{eqnarray} l(\mu_x, \mu_y, \varSigma_{xy}) &=& \log L(\mu_x, \mu_y, \varSigma_{xy} ) \\ &=& \log \displaystyle \prod_{ n = 1 }^N f(x_n,y_n \ | \ \mu_x, \mu_y, \varSigma_{xy} ) \\ &=& \displaystyle \sum_{ n = 1 }^N \log f(x_n,y_n \ | \ \mu_x, \mu_y, \varSigma_{xy} ) \end{eqnarray}$

この対数尤度関数を最大化するパラメータを求めることになる。

1変量正規分布の最尤値が標本平均、標本分散であったように、多変量正規分布についても、最尤値は標本平均、標本分散、標本共分散として解析的に求まる。

$\begin{eqnarray} \mu_x &=& \frac{1}{N} \displaystyle \sum_{ n = 1 }^N x_n &=& \bar{x} \\ \\ \mu_y &=& \frac{1}{N} \displaystyle \sum_{ n = 1 }^N y_n &=& \bar{y} \\ \end{eqnarray}$

$\begin{eqnarray} \varSigma_{xy} &=& \frac{1}{N} \displaystyle \sum_{ n = 1 }^N \left( \begin{array}{c} x_n - \mu_{x} \\ y_n - \mu_{y} \\ \end{array} \right) \left( \begin{array}{c} x_n - \mu_{x} \\ y_n - \mu_{y} \\ \end{array} \right)^{ \mathrm{ T } } \end{eqnarray}$

（参考サイト）
https://home.hiroshima-u.ac.jp/tkurita/lecture/prnn/node7.html

これは、EMアルゴリズムのMステップで、負担率を「1」とおいたものに等しい。

「mvtnorm」パッケージの「rmvnorm」関数を使い、2変量正規分布に従う乱数を作って、実際に計算してみる。
乱数生成時のパラメータは以下のとおりとした（seedは"123"を指定）
f:id:ichou1:20181116190100p:plain

data
               [,1]        [,2]
  [1,]  0.207367739  2.53964502
  [2,]  3.204346438  3.14101678
  [3,]  1.182840469  6.43012997
(snip)
[498,]  1.442972377  2.82004961
[499,]  2.513938298  0.29779923
[500,]  0.260908379  2.50161864

基本統計量の確認。これが最尤値を与えるパラメータとなる。

# 平均
mean(data[,1])
[1] 1.109889

mean(data[,2])
[1] 2.909105

# 標本分散
variance(data[,1])
[1] 1.891138

variance(data[,2])
[1] 4.047196

# 標本共分散
covariance(data[,1], data[,2])
[1] -0.1962958

対数尤度を計算してみる。

対数尤度を計算する関数

log_likelihood_Bivar_norm <- function(param, obs){

  xmean    <- param[1]
  ymean    <- param[2]
  x_sigma2 <- param[3]
  y_sigma2 <- param[4]
  xy_cov   <- param[5]

  # 標本数
  n <- length(obs[,1])

  # 標本共分散
  Sigma <- matrix(c(x_sigma2, xy_cov, xy_cov, y_sigma2), nrow=2, byrow=TRUE)

  # 対数尤度(Log-Likelihood)
  LL <- 0
  for (i in 1:n) {
     LL <- LL + log(dmvnorm(c(obs[i,1], obs[i,2]), mean=c(xmean, ymean), sigma=Sigma))
  }
  return (LL)
}

この関数を使って、最尤値となるパラメータでの対数尤度を計算。

log_likelihood_Bivar_norm(c(1.109889, 2.909105, 1.8911375, 4.0471964, -0.1962958), data)
[1] -1926.477

検算として、optim関数を使って、最尤値となるパラメータを推定してみる（初期パラメータは適当に選択）
ここで、最適化手法（method）を変えて試してみた。

optim(par=c(1,1,2,3,1), fn=log_likelihood_Bivar_norm, obs=data, control=list(fnscale=-1), method=" ??? ")

各methodの概要
パラメータ推定（最尤推定法） - ichou1のブログ

結果。
CG（共役勾配法）以外は収束。
f:id:ichou1:20181116212216p:plain

パラメータ推定値と対数尤度。
１行目は、標本平均、標本分散、標本共分散と、その対数尤度。
収束はしなかったが、CG（共役勾配法）の推定値が一番近い結果となった。
f:id:ichou1:20181116213716p:plain

SANN（シミュレーテッド・アニーリング法）に関しては、なかなか結果が返ってこなかった。
関数が呼ばれた回数を見てみると以下のとおり。
f:id:ichou1:20181116211310p:plain

2018-11-14

統計学における変量と変数の違い

「RとS-PLUSによる多変量解析」という本を読んだ時のメモ。

以下、重回帰分析に関するくだりを引用。
書籍では、二行目は「この手法は、従属変数と複数の説明変数の間の関係を探るために使われる」という記載であるが、「従属変数」を「目的変数」に書き換えた

重回帰分析は、単線形回帰モデルを一般化し、説明変数を二つ以上に拡張したモデルである。
この手法は、目的変数と複数の説明変数の間の関係を探るために使われる。
（中略）
ここで特に注意が必要なのは、説明変数は厳密に言うとまったくランダムな変数とは見なされないことであり、従って重回帰分析は、本質的には目的変数のみをランダム変数とする単変量の解析技術であることである。
そこでこの技術を、本来の多変量の分析手順と区別するため、多変数解析と呼ぶことがある。

「変量」というのは、「ランダムな変動を伴う変数」のこと。

統計学は、観測値からモデルを導き出すものであった。
その観測値が不確定な偶然現象（ランダムな現象）によるものであれば、それは「変量」と呼ぶ。
「変量」でなければ何なのか。さらに引用。

説明変数は厳密には定数と仮定されている。
すなわち、それらはランダム変数と見なされない。
ただし現実には、その仮定が満たされるのはまれであり、重回帰分析の結果は、説明変数の観測値に条件付けられていると解釈される。

厳密には「定数」としての扱いになるようだが、実際には、「変数」の扱いでよさそう。
「重回帰分析」は「単変量多変数解析」にあたるということになる。

「変数」同士の関係性についても、意識して使い分けるようにする必要がある。

条件付けられる（因果関係がある）
相関がある

もう一度、同じ箇所を引用。

重回帰分析は、単線形回帰モデルを一般化し、説明変数を二つ以上に拡張したモデルである。
この手法は、目的変数と複数の説明変数の間の関係を探るために使われる。
（中略）
重回帰分析の結果は、説明変数の観測値に条件付けられていると解釈される。

最後に英語での表記も書いておく。

変量（variate）
変数（variable）

「variate」と出てきたら、意図して使い分けているのか注意が必要。

2018-11-11

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

音声認識 R

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

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

先ほどの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パッケージ関数の値と一致していることを確認。

2018-11-08

多変量混合正規分布

R 音声認識

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

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

ここで、∑は「変量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

2018-11-06

多変量正規分布（その2）

R 音声認識

前回は2変量に相関がないケースを見たが、今回は相関があるケースを見てみる。

各変量の平均と分散は前回と同じ。相関係数だけ変えてみる。

変量x : 平均「1」、分散「2」の正規分布に従う
変量y : 平均「3」、分散「4」の正規分布に従う
変量x と変量y には正の相関がある（＝共分散が「1」であるとして相関係数は「0.35」）

分散共分散行列は以下のとおりとなる。

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

> Sigma.positive_cor
     [,1] [,2]
[1,]    2    1
[2,]    1    4

set.seed(123); data <- rmvnorm(n=N, mean=c(x_mu, y_mu), sigma=Sigma.positive_cor)
             [,1]       [,2]
 [1,]  0.15661631  2.3779756
 [2,]  3.17597129  3.6033091
 [3,]  1.68913945  6.4304202
(snip)
 [9998,] -0.8679178275  1.856718741
 [9999,]  1.1757277988  4.769811627
[10000,]  3.8768372227  5.398123198

チェック。

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

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

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

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

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

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

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

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

変量xと変量yの平均値での確率密度を求めてみる。

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

相関がないケースでの確率密度は「0.05626977」であった。
改めて確率密度の式を見てみる。

分母の行列式の部分を比べてみる。

# 相関係数が0.35
> sqrt(det(Sigma.positive_cor))
[1] 2.645751

# 相関係数が0
> sqrt(det(Sigma.zero))
[1] 2.828427

相関がある方が、行列式の値が小さくなっている。
これは、確率密度が描く山を上に押し上げている要素だと思われる。

比率で見てみる。

# 確率密度での比較（相関係数が0の場合を基準にして何倍か）
>  0.06015491 / 0.05626977
[1] 1.069045

# 行列式の平方根での比較（相関係数が0の場合を基準にして何倍か）
> 2.828427 / 2.645751
[1] 1.069045

また、指数関数の中の、分散共分散行列の逆行列を比べてみる。

# 相関係数が0.35
> solve(Sigma.positive_cor)
           [,1]       [,2]
[1,]  0.5714286 -0.1428571
[2,] -0.1428571  0.2857143

# 相関係数が0
> solve(Sigma.zero)
     [,1] [,2]
[1,]  0.5 0.00
[2,]  0.0 0.25

これは、確率密度が描く山を楕円から扁平に押しつぶしている要素だと思われる。

視点を変えてみる。変量yに着目した形状。
f:id:ichou1:20181109080326p:plain
左図が相関あり、右図が相関なし

上から見下ろした形状。
f:id:ichou1:20181109080343p:plain
左図が相関あり、右図が相関なし

各変量が互いに関連していれば、同時に分析することで、データの特徴が分かる。

一方、各変量に相関がない場合は、前回見たように、2変量の確率密度は1変量正規分布の確率密度の積になる。
これは同時に分析せずとも、それぞれの変量を別個に調べればよいことを示している。

2018-11-04

多変量正規分布

R 音声認識

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

以下の視点で考えた時、

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

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

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

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

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

ここで、∑は「変量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）」ケース。
「相関がある」ケースではどうか、次回に続く。