ichou1のブログ

主に音声認識、時々、データ分析のことを書く

多変量正規分布(その2)

前回は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」であった。
改めて確率密度の式を見てみる。


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

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

# 相関係数が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変量正規分布の確率密度の積になる。
これは同時に分析せずとも、それぞれの変量を別個に調べればよいことを示している。

多変量正規分布

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

以下の視点で考えた時、

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

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

  1. (1変量)正規分布
  2. (1変量)混合正規分布
  3. 多変量正規分布
  4. 多変量混合正規分布

今回は上記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)」ケース。
「相関がある」ケースではどうか、次回に続く。

パラメータ推定(混合正規分布)

混合分布を使うことで、より柔軟に分布をモデル化できる。

例えば、ピークが2つあるような分布。
f:id:ichou1:20181028085746p:plain

これは以下2つの分布を混合したもの。
分布1 : 正規分布N(2, 5)、混合比「0.6」
分布2 : 正規分布N(10, 5)、混合比「0.4」

この混合分布は音声認識において使われていて、隠れマルコフモデル(HMM)の各状態の出力確率が混合正規分布(混合ガウス分布)に従うとしたモデルを「GMM-HMM」と呼ぶ。

具体的にどのようにモデル化しているかは以下の記事に書いた。
work-in-progress.hatenablog.com

今回は混合正規分布最尤推定に関するメモ。

混合数をKとおく。

混合パラメータ


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

観測データがN個のときの対数尤度関数


\begin{eqnarray}
l(\pi_k, \mu_k, \sigma_k^2 ) &=& \log L(\pi_k, \mu_k, \sigma_k^2 ) \\
&=& \log \displaystyle \prod_{ n = 1 }^N p(x_n \ | \ \pi_k, \mu_k, \sigma_k^2 ) \\
&=& \displaystyle \sum_{ n = 1 }^N \log p(x_n \ | \ \pi_k, \mu_k, \sigma_k^2 ) \\
&=& \displaystyle \sum_{ n = 1 }^N \log \left( \displaystyle \sum_{ k = 1 }^K \pi_k \ f(x_n \ | \ \mu_k, \sigma_k^2 ) \right)
\end{eqnarray}

xが与えられたとき、それに対応する確率密度を各分布ごとに分解した割合を仮置きするために、潜在変数zを導入する。
潜在変数zは、K個の要素からなり、どれか1つだけ1を取り、他は全て0となる。
xについてn番目の値が与えられたとき、k番目のzが1となる確率を負担率と呼び、以下の式で計算する。



\begin{eqnarray}
\gamma(z_{nk}) &=& P(z_k =1 | x_n) \\
\\
&=& \frac{\pi_k N(x_n|\mu_k, \sigma_k^2)}{\displaystyle \sum_{k=1}^K \pi_k N(x_n|\mu_k, \sigma_k^2)}
\end{eqnarray}

分母は、混合分布における、n番目のxの確率密度。
分子は、k番目の分布における、n番目のxの確率密度となるので、負担率は



0 \leq \gamma(z_{nk}) \leq 1 \\
\begin{eqnarray}
\displaystyle \sum_{k=1}^K \gamma(z_{nk}) &=& 1
\end{eqnarray}

この辺は以下のサイトを参考にさせていただいた。
混合ガウスモデルとEM - 人工知能に関する断創録

負担率はN個のxのそれぞれについて求まる。
この負担率を使って、パラメータ(平均、分散、混合比)の最尤値を計算する。



\begin{eqnarray}
\mu_{k \ new} &=& \frac{ \displaystyle \sum_{ n = 1 }^N \left( \ \gamma(z_{nk}) \ x_n \ \right) }{ \displaystyle \sum_{ n = 1 }^N \gamma(z_{nk}) } \\
\\
\sigma_{k \ new}^2 &=& \frac{ \displaystyle \sum_{ n = 1 }^N  \left( \ \gamma(z_{nk}) \ (x_n -\mu_k)^2 \ \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}


負担率を求めるためにはパラメータが必要で、パラメータの最尤値を計算するためには負担率が必要と、循環してしまうので、以下の2ステップに分けて計算する。

  • Eステップ:現在の混合比、平均、分散に対応する負担率を計算する
  • Mステップ:Eステップで計算した負担率をもとに、混合比、平均、分散の最尤値を求める


混合ガウス分布の数が「2」、各分布が以下であるサンプルデータ「500」個用意し、推定の過程を追ってみる。
分布1 : 正規分布、平均「2」、分散「5」、混合比「0.6」
分布2 : 正規分布、平均「10」、分散「5」、混合比「0.4」

rnorm関数に与えるseedは「123」を指定。

 [1]  0.7467384  1.4853075  5.4853777  2.1576616  2.2890962  9.7872438
 [7]  3.0306400 -0.8287629  0.4641503  1.0034695  9.7822757  2.8045682
[13]  2.8961522 10.4597300  8.6216076  5.9956592  9.7052800 -2.3974896
(snip)
[486]  3.741583972  4.045495218 11.760084269  5.637816986 13.701085639
[491] 11.511050698  5.139551842  4.893659583  9.523030728 10.700400827
[496]  9.798809343 12.393746629  2.368595386  8.831393539  3.234662183

密度プロット(赤い線はcurve関数で書いた真の確率密度関数
f:id:ichou1:20181030194850p:plain

「R」には「mixtool」ライブラリで「normalmixEM」関数が用意されているので、これを使えば最尤値を求められる。

検算の意味で最初に値を求めておく。

mixtool_val <- normalmixEM(data, k=2)
最尤値での対数尤度
mixtool_val$all.loglik[length(mixtool_val$all.loglik)]
-1388.84
最尤値(混合比、平均、分散)
mixtool_val$lambda
[1] 0.589741 0.410259    # 混合比 

mixtool_val$mu
[1]  1.969816 10.129432  # 平均

mixtool_val$sigma^2
[1] 4.498699 4.566088    # 分散


先ほどの式で最尤値を100回、計算してみた結果。
初期値は、以下の値を与えた。
分布1 : 混合比0.2、平均1、分散2
分布2 : 混合比0.8、平均5、分散2

f:id:ichou1:20181031163918p:plain

負担率の推移。
f:id:ichou1:20181031163929p:plain

パラメータ推定(ニュートン・ラフソン法)

前回の続き。

ニュートン・ラフソン法」を使って、2変数パラメータの最尤推定値を求めてみる。
以下のサイトを参考にさせていただいた。
http://www012.upp.so-net.ne.jp/doi/sas/numerical/Newton_Raphson/Newton_Raphson.pdf

正規分布の対数尤度関数を再掲。

対数尤度関数(正規分布



\begin{eqnarray}
l(\mu, \sigma^2 ) &=& \log L(\mu, \sigma^2 ) \\
 &=& \log \displaystyle \prod_{ i = 1 }^n f(x_i | \mu, \sigma^2 ) \\
 &= & - \frac{n}{2} \log 2\pi - \frac{n}{2} \log \sigma^2 - \frac{1}{2\sigma^2} \displaystyle \sum_{ i = 1 }^{ n } (x_i - \mu)^2
\end{eqnarray}

対数尤度関数をパラメータで1階偏微分する。


\begin{eqnarray}
\nabla l(\mu, \sigma^2) &=&
\left(
  \begin{array}{c}
    \frac{ \partial l(\mu, \sigma^2 ) }{ \partial \mu } \\
    \frac{ \partial l(\mu, \sigma^2 ) }{ \partial \sigma^2 }
  \end{array}
\right) \\
 &=&
\left(
  \begin{array}{c}
    \frac{1}{\sigma^2} \displaystyle \sum_{ i = 1 }^{ n } (x_i - \mu) \\
    - \frac{n}{2\sigma^2} + \frac{1}{2(\sigma^2)^2} \displaystyle \sum_{ i = 1 }^{ n } (x_i - \mu)^2
  \end{array}
\right)
\end{eqnarray}

2階偏微分を求める。



\begin{eqnarray}
\nabla^2 l(\mu, \sigma^2) &=&
\left(
  \begin{array}{cc}
    \frac{ \partial }{ \partial \mu } \frac{ \partial l(\mu, \sigma^2 )}{ \partial \mu }
    &&
    \frac{ \partial }{ \partial \sigma^2 } \frac{ \partial l(\mu, \sigma^2 )}{ \partial \mu }
    \\
    \frac{ \partial }{ \partial \mu } \frac{ \partial l(\mu, \sigma^2 )}{ \partial \sigma^2 }
    &&
    \frac{ \partial }{ \partial \sigma^2 } \frac{ \partial l(\mu, \sigma^2 )}{ \partial \sigma^2 }
  \end{array}
\right)
\\
&=&
\left(
  \begin{array}{cc}
    - \frac{n}{\sigma^2}
    &&
    - \frac{1}{(\sigma^2)^2} \displaystyle \sum_{ i = 1 }^{ n } (x_i - \mu)
    \\
    - \frac{1}{(\sigma^2)^2} \displaystyle \sum_{ i = 1 }^{ n } (x_i - \mu)
    &&
    \frac{n}{2(\sigma^2)^2 } - \frac {1} {(\sigma^2)^3} \displaystyle \sum_{ i = 1 }^{ n } (x_i - \mu)^2
  \end{array}
\right)
\end{eqnarray}

繰り返し計算の回数をkとおくと、漸化式は以下のとおり。



\begin{eqnarray}
\left(
  \begin{array}{c}
    \mu_{k+1}
    \\
    \sigma_{k+1}^2
  \end{array}
\right)
&=&
\left(
  \begin{array}{c}
  	\mu_k
    \\
    \sigma^2_k
  \end{array}
\right)
-
\left(
\nabla^2 l(\mu_k, \sigma_k^2 )
\right)^{-1}
\nabla l(\mu_k, \sigma_k^2 )
\\
&=&
\left(
  \begin{array}{c}
  	\mu_k
    \\
    \sigma^2_k
  \end{array}
\right)
-
  \frac{\nabla l(\mu_k, \sigma_k^2 )}{\nabla^2 l(\mu_k, \sigma_k^2 ) }
\end{eqnarray}

検算。

データ
# 平均2、分散5の正規分布に従うサンプルを200個
set.seed(123); data <- rnorm(200, mean=2, sd=sqrt(5))

  [1]  0.746738355  1.485307487  5.485377748  2.157661556  2.289096164
  [6]  5.835001897  3.030639969 -0.828762916  0.464150333  1.003469540
 [11]  4.737130109  2.804568177  2.896152207  2.247494077  0.757101438
(snip)
[186]  1.559101297  4.481857217  2.189478346  3.686115522  0.883549109
[191]  2.479514290  1.273980231  2.211495199 -0.002093333 -0.931041334
[196]  6.465904894  3.343225764 -0.797927823  0.633391465 -0.650814055

初期パラメータは、平均1、分散1とおいて、先ほどの式を計算した結果。

f:id:ichou1:20181029185002p:plain
最尤値となるパラメータは、平均「1.980836」、分散「4.425514」


検算。Rの関数で計算した結果。

# 準ニュートン法による最適解
> optim(c(1, 1), fn=log_likelihood_norm, obs=data, method="BFGS", control=list(fnscale=-1), hessian=TRUE)
結果
$par
[1] 1.980839 4.425508

$value
[1] -432.5263

$counts
function gradient 
      31       10 

$convergence
[1] 0

$message
NULL

$hessian
              [,1]          [,2]
[1,] -4.519255e+01  3.119993e-05
[2,]  3.119993e-05 -5.105931e+00

小数点第五位まで一致していることが確認できた。

また、k回目の最尤値でのヘッセ行列を固有値分解してみる。

H_0 <- matrix(c(-200.000000 , -196.167181 , -196.167181 , -977.510617  ), nrow=2)
H_1 <- matrix(c(-152.232794 ,  -35.661895 ,  -35.661895 , -340.743870  ), nrow=2)
H_2 <- matrix(c(-108.715893 ,   -7.280127 ,   -7.280127 , -113.101423  ), nrow=2)
H_3 <- matrix(c( -79.549506 ,   -1.429134 ,   -1.429134 ,  -39.900363  ), nrow=2)
H_4 <- matrix(c( -61.130669 ,   -0.254284 ,   -0.254284 ,  -15.933080  ), nrow=2)
H_5 <- matrix(c( -50.656386 ,   -0.036104 ,   -0.036104 ,   -7.966423  ), nrow=2)
H_6 <- matrix(c( -46.162098 ,   -0.002919 ,   -0.002919 ,   -5.555946  ), nrow=2)
H_7 <- matrix(c( -45.231581 ,   -0.000058 ,   -0.000058 ,   -5.123588  ), nrow=2)
H_8 <- matrix(c( -45.192558 ,   -0.000000 ,   -0.000000 ,   -5.105933  ), nrow=2)
H_9 <- matrix(c( -45.192490 ,   -0.000000 ,   -0.000000 ,   -5.105903  ), nrow=2)

固有値分解の結果
f:id:ichou1:20181029193007p:plain

最尤値をとるパラメータを与えたとき、ヘッセ行列の固有ベクトルは標準基底になっている。
また、各パラメータの2階偏微分の値が固有値になっている。


推定するパラメータが2つの場合、フィッシャー情報行列は以下として定義される。



\begin{eqnarray}
I(\mu, \sigma^2)
&=& - E(\nabla^2 (\mu, \sigma^2 ))
\end{eqnarray}

今回は1セットのデータのみ用いており、パラメータ推定値も1セットのみであるから、フィッシャー情報行列は、対数尤度関数の2階偏微分にマイナスを掛けたものになる。



\begin{eqnarray}
I(\mu, \sigma^2)
&=& - \nabla^2 (\mu, \sigma^2 )
\end{eqnarray}



今回、各繰り返しにおけるフィッシャー情報量と対数尤度は以下のとおりであった。
f:id:ichou1:20181029205949p:plain

平均を固定して、分散を変えてみると以下のとおりであった。
f:id:ichou1:20181029212531p:plain

フィッシャー情報行列は、対数尤度関数のスコア関数(1階偏微分)の分散として定義されるから、以下のことがいえる。

  • スコア関数の分散が大きい(=フィッシャー情報量が大きい)ほど、相対的に尤度は下がる
  • スコア関数の分散が小さい(=フィッシャー情報量が小さい)ほど、相対的に尤度は上がる

パラメータ推定(最尤推定法)



今、手元にn個の観測値(標本)が得られているとする。

統計モデルとして表現するには、

  1. 観測値が従うであろう確率分布(確率密度関数)の形
  2. 確率密度関数のパラメータ

が必要となる。

1番目の確率分布の形が分かったとして、2番目のパラメータを推定する方法の1つである「最尤推定法」に関するメモ。

ある1つの観測値が、他の観測値に影響を及ぼさない(確率変数は独立している、という)場合、それぞれの観測値の確率密度を掛けあわせたもの(同時確率)が尤度関数となる。

観測値つまり確率変数をx、未知パラメータをθとおくと、尤度関数は以下の形で表される。
\[
L(\theta) = \displaystyle \prod_{ i = 1 }^n f(x_i | \theta)
\]

尤度関数すなわち同時確率を最大にするパラメータθの値が「最尤推定値」になる。
小数の乗算はアンダーフローを引き起こすので、尤度関数の対数をとった「対数尤度関数」を最大化するパラメータθを求める。

具体的に見てみる。
確率分布の形が「正規分布」だとすると、推定したいパラメータθ は、
 \theta = \left( \mu , \sigma^2 \right)
となる。
μは平均、σ標準偏差を表す)

正規分布確率密度関数


f(x | \mu, \sigma^2 )=\frac{1}{\sqrt{2\pi\sigma^{2}}} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2} \right)

を対数尤度関数の式にあてはめて、式変形する。


\begin{eqnarray}
\log L(\mu, \sigma^2 ) &=& \log \displaystyle \prod_{ i = 1 }^n f(x_i | \mu, \sigma^2 ) \\
 &= & - \frac{n}{2} \log 2\pi - \frac{n}{2} \log \sigma^2 - \frac{1}{2\sigma^2} \displaystyle \sum_{ i = 1 }^{ n } (x_i - \mu)^2
\end{eqnarray}

導出過程は以下のサイトを参照。
統計学的画像再構成法である

これをRの関数として定義する。
第1項は定数なので、最大化を考えるにあたっては無視できるが、今回は確認用に残しておく。

対数尤度関数(正規分布
# arg1 : vector(size 2), [1]mean, [2]variance  
# arg2 : vector(size N), observation
log_likelihood_norm <- function(param, obs){
  mu <- param[1];
  sigma2 <- param[2];
  term1 <- (-1) * length(obs) * 0.5 * log(2*pi)
  term2 <- (-1) * length(obs) * 0.5 * log(sigma2);
  term3 <- (-1) * 0.5 / sigma2 * sum((obs-mu)^2);

  return ( term1 + term2 + term3 );
}

検算。

平均「2」、分散「5」の正規分布に従う乱数を10個用意する。

> data <- rnorm(10, mean=2, sd=sqrt(5))
> data
 [1] -0.01731340  1.79861872  2.99598784  1.78515021  2.42412228 -1.61575147
 [7]  1.74058535  3.84369888  1.93563703  0.03983286
Rのdnorm関数を使って求めた結果
value <- 0
for (i in 1:length(data)) {
    value <- value + log(dnorm(data[i], mean=2, sd=sqrt(5)))
}                             
value
# -19.80805
 
value <- 1
for (i in 1:length(data)) {
    value <- value * dnorm(data[i], mean=2, sd=sqrt(5))
}                             
log(value)
# -19.80805
作成した関数での結果
log_likelihood_norm(c(2, 5), data)
# -19.80805

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

値が最大になるパラメータをRのoptim関数を使って探してみる。

最適解(サンプル数10)
optim(c(0.5, 1), fn=log_likelihood_norm, obs=data, control=list(fnscale=-1), hessian=FALSE)
# $par
# [1] 1.493105 2.313280

# $value
# [1] -18.38531

# $counts
#  function gradient 
        73       NA 

# $convergence
# [1] 0

出てきた最適値は、平均が「1.493105」、分散が「2.313280」

統計量(サンプル数10)
mean(data)
# [1] 1.493057

# 不偏分散
var(data)
# [1] 2.571645

# 標本分散
variance <- function(x) { var(x) * (length(x)-1) / length(x) }
variance(data)
# [1] 2.31448

サンプル数を100に増やしてみる。

最適解(サンプル数100)
data <- rnorm(100, mean=2, sd=sqrt(5))
optim(c(0.5, 1), fn=log_likelihood_norm, obs=data, control=list(fnscale=-1), hessian=FALSE)
# $par
# [1] 1.918190 4.851344

# $value
# [1] -220.8375

# $counts
# function gradient 
        59       NA 

# $convergence
# [1] 0

出てきた最適値は、平均が「1.918190」、分散が「4.851344」

統計量(サンプル数100)
mean(data)
# [1] 1.917897
var(data)
# [1] 4.898471
variance(data)
# [1] 4.849486

さらにサンプル数を増やしてみる。

最適解(サンプル数1000)
data <- rnorm(1000, mean=2, sd=sqrt(5))
optim(c(0.5, 1), fn=log_likelihood_norm, obs=data, control=list(fnscale=-1), hessian=FALSE)
# $par
# [1] 1.891001 4.752660

# $value
# [1] -2198.32

# $counts
# function gradient 
        63       NA 

# $convergence
# [1] 0
統計量(サンプル数1000)
mean(data)
# [1] 1.891252
var(data)
# [1] 4.7577
variance(data)
# [1] 4.752943

サンプル数が多くなるほど、標本分散と不偏分散の差が少なくなっている。

対数尤度関数が返す値の3次元散布図を書いてみる。
meanが「2の前後3、-1から5まで1刻みの7通り」、sigma2は「5の前後3、2から8まで1刻みの7通り」の49パターンとした。

f:id:ichou1:20181017164527p:plain

一番高くなるところのパラメータが最尤推定値となる。

最尤推定値を求めるにあたって、Rのoptim関数では以下の最適化アルゴリズムが使える模様。

method="Nelder-Mead"  # シンプレックス法、関数が微分できなくても可
method="BFGS"         # 準ニュートン法
method="L-BFGS-B"     # 準ニュートン法の変種(推定するパラメータのlower/upperを指定)
method="CG"           # 共役勾配法
method="SANN"         # シミュレーテッド・アニーリング法、関数が微分できなくても可

関数の最大・最小化 - RjpWiki

methodパラメータの省略時デフォルトは「Nelder-Mead」(シンプレックス法
optim function | R Documentation

The default method is an implementation of that of Nelder and Mead (1965), that uses only function values and is robust but relatively slow.
It will work reasonably well for non-differentiable functions.

(参考)多変量正規分布のパラメータ推定で、optim関数のmethodを変えて比較してみた結果
パラメータ推定(多変量正規分布) - ichou1のブログ

ヘッセ行列の固有値による凸性の判定

最適化を考える時、どの方向に向かって修正していくかを決める必要がある。

扱う対象が関数化できるなら、その最大値(あるいは最小値)を求める問題に帰着するので、関数の凸性が大事な概念になってくる。

今回は2変数関数を対象に、2階偏微分であるヘッセ行列を使って凸性を調べる方法を見てみる。

以下のサイトを参考にさせていただいた。
http://www2.kaiyodai.ac.jp/~yoshi-s/Lectures/Optimization/2013/lecture_1.pdf

偏微分の順番で微分結果は変わらないので、ヘッセ行列は対称行列になる。

3パターンの関数で見てみる。

[パターン1] 狭義凸関数

ヘッセ行列が正定値(固有値が全て0より大きい)
グラフが全ての方向で丸く突き出ている

[パターン2] 凸関数(狭義凸関数ではない)

ヘッセ行列が半正定値(固有値が全て0以上)
2点を結ぶ線分とグラフが重なる部分がある

[パターン3] 凸面と凹面が混在する関数

ヘッセ行列が不定固有値が正と負の両方を持つ)


[パターン1] 狭義凸関数

f:id:ichou1:20181014084615p:plain

勾配ベクトル


\begin{eqnarray}
\nabla f(x,y) = \begin{pmatrix} 2x \\ 2y \end{pmatrix}
\end{eqnarray}


ヘッセ行列


\begin{eqnarray}
\nabla^2 f(x,y) = \begin{pmatrix} 2 & 0 \\ 0 & 2 \end{pmatrix}
\end{eqnarray}


ヘッセ行列を固有値分解する。
固有値は「2」であることを確認

> eigen(matrix(c(2, 0, 0, 2), nrow=2, ncol=2, byrow=TRUE))
eigen() decomposition
$values
[1] 2 2

$vectors
     [,1] [,2]
[1,]    0   -1
[2,]    1    0

角度を変えてみる
y軸が奥行きと重なるように見たところ

f:id:ichou1:20181014084715p:plain
x軸が奥行きと重なるように見たところ

f:id:ichou1:20181014084730p:plain

全ての方向で丸く突き出た形をしている。

[パターン2] 凸関数(狭義凸関数ではない)

f:id:ichou1:20181014085017p:plain

勾配ベクトル


\begin{eqnarray}
\nabla f(x,y) = \begin{pmatrix} 1 \\ 2y \end{pmatrix}
\end{eqnarray}


ヘッセ行列


\begin{eqnarray}
\nabla^2 f(x,y) = \begin{pmatrix} 0 & 0 \\ 0 & 2 \end{pmatrix}
\end{eqnarray}


ヘッセ行列を固有値分解する。
固有値は「2」と「0」であることを確認

> eigen(matrix(c(0, 0, 0, 2), nrow=2, ncol=2, byrow=TRUE))
eigen() decomposition
$values
[1] 2 0

$vectors
     [,1] [,2]
[1,]    0   -1
[2,]    1    0

角度を変えてみる。
y軸が奥行きと重なるように見たところ

グラフの底が、線分の形をしている。
f:id:ichou1:20181014085026p:plain
x軸が奥行きと重なるように見たところ

f:id:ichou1:20181014085039p:plain

[パターン3] 凸面と凹面が混在する関数

f:id:ichou1:20181014085301p:plain

勾配ベクトル


\begin{eqnarray}
\nabla f(x,y) = \begin{pmatrix} 2x+2y \\ 2x-2y \end{pmatrix}
\end{eqnarray}


ヘッセ行列


\begin{eqnarray}
\nabla^2 f(x,y) = \begin{pmatrix} 2 & 2 \\ 2 & -2 \end{pmatrix}
\end{eqnarray}


ヘッセ行列を固有値分解する。
固有値は「2.828427」と「-2.828427」であることを確認

> eigen(matrix(c(2, 2, 2, -2), nrow=2, ncol=2, byrow=TRUE))
eigen() decomposition
$values
[1]  2.828427 -2.828427

$vectors
           [,1]       [,2]
[1,] -0.9238795  0.3826834
[2,] -0.3826834 -0.9238795

角度を変えてみる
y軸が奥行きと重なるように見たところ

f:id:ichou1:20181014085310p:plain
x軸が奥行きと重なるように見たところ

f:id:ichou1:20181014085326p:plain

凸面と凹面が混在している。

行列の低ランク近似と近似度

特異値分解を使った「行列の低ランク近似」に関するメモ。

まずは、視覚的に確認しやすい画像データで試してみる。

高さ200px、幅320px、インデックスカラーのRAWデータを用意する。
f:id:ichou1:20181010153239p:plain

1pxを1バイトで表現、0〜255の範囲の値をとる200行 x 320列の行列ができる。

以下は左上を拡大した画像。
先頭6バイト分のみ値を添えてみた。値が小さいほど黒に近く、値が大きいほど白に近くなる。
f:id:ichou1:20181010160014p:plain

CSVに落としたものをRに読み込み、

A <- data.matrix(read.table("200_320_gray_scale.csv", header=F, sep=","))

ランク数を「r」とおいて、特異値分解する。

A.svd <- svd(A)

u <- A.svd$u[,1:r, drop=F]
d <- diag(A.svdA$d)[1:r, 1:r, drop=F]
v <- A.svd$v[,1:r, drop=F]

そして、再度、csvに書き戻す。

write.table(round(u %*% d %*% t(v)),
                 str_c("200_320_gray_scale_low_rank_approximate_", r), ".csv"),
                 sep=",", row.names=FALSE, col.names=FALSE)
r = 1

f:id:ichou1:20181010161502p:plain

r = 1 (拡大)

f:id:ichou1:20181010162151p:plain

r = 10

f:id:ichou1:20181010162446p:plain

r = 10 (拡大)

f:id:ichou1:20181010163421p:plain

r = 20

f:id:ichou1:20181010162541p:plain

r = 20 (拡大)

f:id:ichou1:20181010163239p:plain


続けて、フロベニウスノルム(全成分の二乗和の平方根)を確認してみる。

フロベニウスノルム(元々の行列)
norm(A, "F")  # 26607.34 

続けて、近似した行列のフロベニウスノルム

フロベニウスノルム(ランク1、5、10、15、20、25、50、75、100で近似した行列)
norm(A_rank1,   "F")  # 25696.45
norm(A_rank5,   "F")  # 26186.34
norm(A_rank10,  "F")  # 26333.71
norm(A_rank15,  "F")  # 26412.39
norm(A_rank20,  "F")  # 26460.39
norm(A_rank25,  "F")  # 26493.76
norm(A_rank50,  "F")  # 26570.16
norm(A_rank75,  "F")  # 26594.95
norm(A_rank100, "F")  # 26603.22

近似に使うランク数を大きくするほど、フロベニウスノルムも大きくなっていることが確認できる。
尚、特異値の上位20は以下のとおりであった。

> A.svd$d[1:20]
 [1] 25695.8853  3754.6842  2374.1272  1782.2953  1588.0870  1466.8724
 [7]  1308.7928  1248.0294  1112.3381  1042.8967   983.9771   963.9730
[13]   916.1656   841.3105   837.3214   794.9816   735.2430   722.9311
[19]   685.1815   634.1618

下位20は以下のとおり。

> A.svd$d[180:200]
 [1] 15.898742 15.417982 14.987853 14.435300 14.174012 13.721447 13.289164
 [8] 12.747505 12.436280 11.677728 11.473563 11.297093 10.792612 10.623457
[15]  9.714841  9.523776  9.005887  8.843195  8.200154  8.091183  7.484766
各ランクの特異値(ランク1〜200)

f:id:ichou1:20181012210530p:plain

各ランクの特異値(ランク11〜200)

f:id:ichou1:20181012213726p:plain


(参考)特異値を全部足すと以下のとおり。

> sum(A.svd$d)
[1] 75484.84
近似に使うランク数とフロベニウスノルム

f:id:ichou1:20181012210506p:plain


フロベニウスノルムは以下のやり方でも求まる。

列方向または行方向の内積を求めて(分散共分散行列になる)、

AT_A <- t(A) %*% A # column-wise 320row x 320col
A_AT <- A %*% t(A) # row-wise    200row x 200col    

トレース(対角成分の和)の平方根を求める。

> sqrt(sum(diag(AT_A)))
[1] 26607.34
> sqrt(sum(diag(A_AT)))
[1] 26607.34