tensorflowメモ(手書き文字認識その2)

前回生成したモデルを掘り下げてみる。
重み付けを行うパラメータが「final_w」に入っているとする。

final_w.shape
(64, 10)

# 最小値
numpy.min(final_w)
-1.68

# 最大値
numpy.max(final_w)
1.55

「final_w」をグレースケール変換して視覚化してみる。
変換手順は以下のとおり。

  • 各値に、「1.68」を加算する(結果は 0 から 3.36(= 1.68 x 2)の範囲になる)
  • 各値を「75.893」倍(= 255 / 3.36)する(結果は 0 から 255の範囲になる、元々のゼロは「127.5」になる)
  • 各値の端数を切り捨てる

各クラス(0から9)ごとの重みパラメータの変換結果は以下のとおり。

重みパラメータを可視化

f:id:ichou1:20181211112830p:plain

中間値が灰色だとすると、白に近い部分は正、黒に近い部分は負に該当する。

これをゼロとラベル付けした画像に適用してみる。
f:id:ichou1:20181211113542p:plain

重みパラメータの白に近い部分は正の値、黒の部分は負の値として作用する。
outputとして出力された値を足し上げたものを使って各クラスの尤度を判定する。

入力画像をモデルに通した後の各クラスの値
0:  2.6651433
1: -1.6024852
2: -0.6885104
3: -0.2711072 
4: -0.0937684 
5:  0.1937362 
6: -0.3147157 
7: -0.5869334 
8:  0.0727276
9:  0.6259142

クラス"0"の値が一番大きく、次いで、クラス"9"の値が大きくなっている。

クラス"9"の重みを見てみると、左下部分に入力があるとマイナスに作用するようになっている。
クラス"9"のoutputでも左下部分がマイナスとして出力されており、結果、インプット画像がクラス"9"であることの尤度を押し下げている。
(クラス"0"とクラス"9"の特徴分けをしているのは左下部分と言える)

tensorflowメモ(手書き文字認識その1)

ライブラリ「scikit-learn」の手書き文字データセット「digits」を使って、手書き文字認識を試してみる。

入力データ

データセットの構成は以下のとおり。

  • 8 x 8ピクセルの画像
  • 「0」から「9」までの手書き文字
  • 明暗を0から16までの値で表現
  • 画像数は1797

1番めの画像のデータ表現。
64個(=8x8)の値が1列に並んで格納されている。

array([ 0,  0,  5, 13,  9,  1,  0,  0,
        0,  0, 13, 15, 10, 15,  5,  0,
        0,  3, 15,  2,  0, 11,  8,  0,
        0,  4, 12,  0,  0,  8,  8,  0,
        0,  5,  8,  0,  0,  9,  8,  0,
        0,  4, 11,  0,  1, 12,  7,  0,
        0,  2, 14,  5, 10, 12,  0,  0,
        0,  0,  6, 13, 10,  0,  0,  0])

画像にしてみる。
f:id:ichou1:20181208085216p:plain

明度を反転。
f:id:ichou1:20181208085521p:plain

これは数字のゼロにあたるものとして、ラベル付けをしている。

学習する際は、0から1の範囲になるようスケーリングする(16で割る)。

array([0.    , 0.    , 0.3125, 0.8125, 0.5625, 0.0625, 0.    , 0.    ,
       0.    , 0.    , 0.8125, 0.9375, 0.625 , 0.9375, 0.3125, 0.    ,
       0.    , 0.1875, 0.9375, 0.125 , 0.    , 0.6875, 0.5   , 0.    ,
       0.    , 0.25  , 0.75  , 0.    , 0.    , 0.5   , 0.5   , 0.    ,
       0.    , 0.3125, 0.5   , 0.    , 0.    , 0.5625, 0.5   , 0.    ,
       0.    , 0.25  , 0.6875, 0.    , 0.0625, 0.75  , 0.4375, 0.    ,
       0.    , 0.125 , 0.875 , 0.3125, 0.625 , 0.75  , 0.    , 0.    ,
       0.    , 0.    , 0.375 , 0.8125, 0.625 , 0.    , 0.    , 0.    ])

全データ(1797個)のうち2割を検証用とし、残りの1437個をトレーニング用とする。

モデル

モデルを以下のとおり表現する。

trainY = X * W + b
  • trainY : 計算結果を格納する行列(minibatch-size x 10)
  • X : 入力データを格納する行列(minibatch-size x 64)
  • W : 重みを格納する行列(64 x 10)
  • b : バイアスを格納する1次元配列(要素数10)

tensorflowとしては以下のとおり定義する。

  • X : tf.placeholder(tf.float32, [None, 64])
  • W : tf.Variable(tf.zeros([64, 10]))
  • b : tf.zeros([10])

「W」と「b」については、オール0で初期化しているが、これでも最適解は求まる。

学習

損失関数として、計算結果と正解ラベルとのクロスエントロピーを計算する。

loss = - tf.reduce_sum(tf.log(tf.nn.softmax(trainY)) * "正解ラベル(one-hot表現)", axis=1)

この値を最小化する問題として、オプティマイザに渡す(ここでは「勾配降下法」を使う)

tf.train.GradientDescentOptimizer("learning-rate").minimize(loss)

ミニバッチサイズを「128」、学習率を「0.001」、エポックを「100」回に設定してトレーニングを実施。
各エポックにおける、テストデータに対する認識精度は以下のとおりであった。

学習 1回目完了のモデル :  accuracy = 80.6(%)
学習11回目完了のモデル :  accuracy = 88.6(%)
学習21回目完了のモデル :  accuracy = 90.8(%)
学習31回目完了のモデル :  accuracy = 91.7(%)
学習41回目完了のモデル :  accuracy = 93.1(%)
学習51回目完了のモデル :  accuracy = 93.6(%)
学習61回目完了のモデル :  accuracy = 93.9(%)
学習71回目完了のモデル :  accuracy = 94.7(%)
学習81回目完了のモデル :  accuracy = 95.0(%)
学習91回目完了のモデル :  accuracy = 95.0(%)

モデルの初期パラメータをオールゼロに設定したが、1回の学習で80%、最大95%の認識精度が出ている。
バイアスはオールゼロのままで、更新されることはなかった。


損失計算の過程を追ってみる。
11回目の学習が完了したモデルに、冒頭の画像データ("ゼロ"とラベリングされたもの)を渡すとする。

入力画像1つをモデルに通した後の結果
data_0 = np.dot(np.reshape(X[0], (1,64)), train_W) + train_b
[[ 2.6651433 -1.6024852 -0.6885104 -0.27110717 -0.09376842 0.19373618 -0.31471574 -0.5869334 0.07272757 0.62591416  ]]
softmax関数を通した結果
data_0_softmax = sess.run(tf.nn.softmax(data_0))
[[ 0.64753246 0.00907515 0.02263541 0.03436087 0.0410281 0.05469443 0.03289465 0.02505547 0.04846071 0.08426274 ]]
対数をとる
data_0_softmax_log = sess.run(tf.log(data_0_softmax)
[[ -0.43458635 -4.702215 -3.78824 -3.3708367 -3.193498 -2.9059935 -3.4144452 -3.6866632 -3.0270019 -2.4738154 ]]
入力画像に紐付けられる正解ラベル(one-hot表現)
Y[0]
[1. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
計算結果と正解ラベルとのクロスエントロピーを計算
data_0_softmax_log_xent = data_0_softmax_log * np.reshape(Y[0], (1,10))
[[ -0.43458635  -0.  -0.  -0.  -0.  -0.  -0.  -0.  -0.  -0. ]]
2次元目のレベルで足し合わせ
data_0_softmax_log_xent_sum = sess.run(tf.reduce_sum(data_0_softmax_log_xent, axis=1))
[ -0.43458635 ]

ここで尤度計算は、確率同士の掛け算を扱うことになるが、アンダーフローを引き起こすので、対数尤度に置き換えて最適化を考えている。

正解に近いほど値は大きくなり(=絶対値が小さくなる)、正解から乖離するほど値は小さくなる(絶対値が大きくなる)。

マイナスを掛けて最小化問題に置き換える
[ 0.43458635 ]

ミニバッチサイズ分、これと同様の値が得られる。

学習においては、この値が小さくなるよう、モデルのパラメータをチューニングする。

確率分布の推定

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

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

以前、多変量混合正規分布のパラメータ推定の際に使った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つの方法に分けられる模様。

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

混合分布モデルや階層型ニューラルネットワークは、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」と確率が高くなっており、想定とは違う方向に「尤もらしい」を押し上げるように作用している。

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

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

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

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

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

以下、重回帰分析に関するくだりを引用。

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

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

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

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

  • 変量(variate)
  • 変数(variable)


同じように、変量の関係性についても、意識して使い分けるようにする。

  • 相関がある
  • 条件付けられる(因果関係がある)

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

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

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

変量の数は「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}
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}
\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、2で、それぞれ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パッケージ関数の値と一致していることを確認。

多変量混合正規分布

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

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

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


混合比は各クラスタの確率密度を押し下げるように作用する。
任意のxとyの確率密度は、各クラスタでの確率密度の和になる。


具体的な値をあてはめてみる。
パラメータは下表のとおりとした。
f:id:ichou1:20181114210802p:plain

分布1
分布2
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