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

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

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