パラメータ推定(多変量正規分布)

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

変量の数は「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関数を使って、最尤値となるパラメータを推定してみる(初期パラメータは適当に選択)
ここで、最適化手法を変えて試してみた。

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

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

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


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