様々な多変量正規分布の最尤推定量の導出方法をみていく。
平均ベクトルと共分散行列で解説したコレスキー分解を用いた方法以外にも、スペクトル分解、帰納法、行列の微分を用いた導出方法をみていく。
多変量正規分布に関する最尤推定については、平均ベクトルと共分散行列の最尤推定量を参照されたい。
対数尤度関数の最大化
多変量正規分布の対数尤度関数は次のように表現することができる。
\begin{align}\label{eq1}f(\boldsymbol{G}) = -N\log |\boldsymbol{G}| \mathrm{tr}(\boldsymbol{G}^{-1}\boldsymbol{D}).\tag{1}\end{align}
\eqref{eq1}の構成に関する詳細は平均ベクトルと共分散行列の最尤推定量を参照されたい。
共分散行列の最尤推定量を得るためには、\eqref{eq1}を\(\boldsymbol{G}\)に関して最大化させる必要がある。また、この最大化させる\(\boldsymbol{G}\)が最尤推定量となる。
ここでは、様々な最大化の方法について解説する。
スペクトル分解を用いた導出
まず、スペクトル分解を用いた導出方法をみていく。コレスキー分解を用いた導出方法と同様に、\(\boldsymbol{D} = \boldsymbol{EE}^T\)、\(\boldsymbol{E}^T\boldsymbol{G}^{-1}\boldsymbol{E} = \boldsymbol{H}\)とおく。このとき、\eqref{eq1}は次のように書き換えられる。
\begin{align}\label{eq2} f = -N\log |\boldsymbol{D}| + N\log |\boldsymbol{H}| -\mathrm{tr}(\boldsymbol{H}).\tag{2}\end{align}
スペクトル分解を用いると、行列\(\boldsymbol{H}\)について以下がいえる。
\begin{align}|\boldsymbol{H}| &= \prod_{i=1}^p\lambda_i,\\ \mathrm{tr}(\boldsymbol{H}) &= \sum_{i=1}^p\lambda_i.\end{align}
したがって、\eqref{eq2}は次となる。
\begin{align}f &= -N\log |\boldsymbol{D}| + N\log \prod_{i=1}^p \lambda_i -\sum_{i=1}^p \lambda_i\\ &= -N\log |\boldsymbol{D}| + \sum_{i=1}^p\left(N\log\lambda_i - \lambda_i\right).\end{align}
上式を\(\lambda_i,\ i=1, \ldots, p\)に関して最大化させることを考える。
\begin{align}\cfrac{\partial f}{\partial \lambda_i} &=\cfrac{\partial}{\partial \lambda_i}\left[ -N\log |\boldsymbol{D}| + \sum_{j=1}^p\left(N\log\lambda_j - \lambda_j\right)\right]\\&= N\cfrac{1}{\lambda_i} -1 \end{align}
であることから
\begin{align}&N\cfrac{1}{\lambda_i} -1 = 0\\ &\Leftrightarrow \lambda_i = N.\end{align}
故に、\(\lambda_i=N, \ i=1, \ldots, p\)から固有方程式に関して\( (N - \lambda)^p = 0 \Leftrightarrow |N\boldsymbol{I} - \lambda\boldsymbol{I}| = 0 \)が成り立つ。よって、\eqref{eq2}は\(|\boldsymbol{H}| = N\boldsymbol{I}\)のときに最大化される。また、その最大値は次となる。
\begin{align}-N\log|\boldsymbol{D}| + \sum_{i=1}^p(N\log N -N) = -N\log|\boldsymbol{D} +pN\log N pN.\end{align}
帰納法を用いた導出
次に帰納法を用いた導出方法を紹介する。正方行列\(\boldsymbol{A}\)を次のように分割する。
\begin{align}\boldsymbol{A} = \begin{pmatrix}\boldsymbol{A}_{11} & \boldsymbol{A}_{12} \\ \boldsymbol{A}_{21} & \boldsymbol{A}_{22} \end{pmatrix}.\end{align}
ブロック行列の行列で解説したように次が成り立つ。
\begin{align}\label{eq3} |\boldsymbol{A}| = | \boldsymbol{A}_{11} - \boldsymbol{A}_{12}\boldsymbol{A}_{22}^{-1}\boldsymbol{A}_{21}| \cdot | \boldsymbol{A}_{22}|.\tag{3}\end{align}
ここで、\eqref{eq2}の\(\boldsymbol{H}\)を次のようにブロック行列で表現する。
\begin{align}\boldsymbol{H} = \begin{pmatrix}\boldsymbol{H}_{i-1} & \boldsymbol{h}_{(i)} \\ \boldsymbol{h}_{(i)}^T & h_{ii} \end{pmatrix}, \ \ \ \ i = 1, \ldots, p.\end{align}
このように、\(i=1, \ldots, p\)について帰納法により示していく。
\(i=1\)のとき、
\begin{align}f &= -N\log|\boldsymbol{D}| +N\log|h_{11}| - \mathrm{tr}(h_{11})\\&= -N\log | \boldsymbol{D}| +N\log h_{11} - h_11. \end{align}
したがって、スペクトル分解を用いた導出方法と同様に、\(h_{11} = N\)のときに最大化される。
次に\(i=k-1\)のとき、\(\boldsymbol{H}_{k-1} = N\boldsymbol{I}_{k-1}\)のときに、\eqref{eq2}が最大化されると仮定し、\(i=k\)のとき、
\begin{align}\boldsymbol{H}_k = \begin{pmatrix}\boldsymbol{H}_{k-1} & \boldsymbol{h}_{(k)}\\ \boldsymbol{h}_{(k)}^T & h_{kk}\end{pmatrix},\ \ \ \ \boldsymbol{H}_{k-1} : (k-1)\times (k-1)\end{align}
について、\eqref{eq2}は\(\boldsymbol{H}_k = N\boldsymbol{I}_k\)のときに最大化されることを示す。ここで、\(\boldsymbol{H}_k\)は正定値行列であることより、そのブロック行列である\(\boldsymbol{H}_{k-1}\)も正定値行列である（ブロック行列の正定値性より）。また正定値性より\(\boldsymbol{H}_{22\cdot1} = h_{kk} - \boldsymbol{h}_{(k)}^T\boldsymbol{H}_{k-1}\boldsymbol{h}_{(k)} > 0\)である。よって、\eqref{eq3}に\(\boldsymbol{H}_k\)を適用すると次を得る。
\begin{align} |\boldsymbol{H}_k| &= |\boldsymbol{H}_{22\cdot 1}|\cdot| \boldsymbol{H}_{k-1}| \\&= \boldsymbol{H}_{22\cdot1} |\boldsymbol{H}_{k-1}|.\end{align}
したがって、\eqref{eq2}は次のように書き換えられる。
\begin{align}f &= N\log|\boldsymbol{D}| + N\log \boldsymbol{H}_{22\cdot1} |\boldsymbol{H}_{k-1}| - \left[\mathrm{tr}(\boldsymbol{H}_{k-1}) + h_{kk}\right]\\ &= -N\log |\boldsymbol{D}| + \left[N\log|\boldsymbol{H}_{k-1}| - \mathrm{tr}(\boldsymbol{H}_{k-1})\right] +N\log(h_{kk} - \boldsymbol{h}_{(k)}^T \boldsymbol{H}_{k-1}^{-1}\boldsymbol{h}_{(k)}) - h_{kk}.\end{align}
今、\(\boldsymbol{H}_{k-1}\)と\(h_{kk}\)を固定したとき、\(\log\)が単調増加関数であることから、上の関数の\(\boldsymbol{h}_{(k)}\)に関する最大化は\(\boldsymbol{h}_{(k)} = \boldsymbol{0}\)のときにおこる。このとき、上式は
\begin{align} f &= -N\log|\boldsymbol{D}| + \left[N\log |\boldsymbol{H}_{k-1} -\mathrm{tr}(\boldsymbol{H}_{k-1})|\right] + (N\log h_{kk} - h_{kk})\end{align}
となる。まず、この式の\(\boldsymbol{H}_{k-1}\)に関する最大化を考える。この式の第2項は\(i=k-1\)のときの\eqref{eq2}の第2項、第3項と一致する。そのため、\(\boldsymbol{H}_{k-1} = N\boldsymbol{I}_{k-1}\)のときに最大化される。また、スペクトル分解を用いた導出方法と同様に、\(h_{kk} =N\)のときに最大化される。したがって、\(i=k\)のとき\(\boldsymbol{H}_k = N\boldsymbol{I}_k\)によって\eqref{eq2}は最大化される。故に、帰納法により、\(i=1, \ldots, p\)について、\eqref{eq2}は\(\boldsymbol{H} = N\boldsymbol{I}\)のときに最大である。また、その最大値は
\begin{align}f -N\log |\boldsymbol{D} + pN \log N -pN\end{align}
である。
行列の微分を用いた導出
\eqref{eq1}の第2項と第3項を次の式で表す。
\begin{align}f(\boldsymbol{C}; \boldsymbol{D}) = N\log|\boldsymbol{C}| -\mathrm{tr}(\boldsymbol{CD}).\end{align}
\eqref{eq1}から、\(\boldsymbol{C} = \boldsymbol{G}^{-1}\)という関係があることがわかる。関数\(f(\boldsymbol{C}; \boldsymbol{D})\)を最大化させるために、次の行列微分を用いる。
\begin{align}\cfrac{\partial \mathrm{tr}(\boldsymbol{CD}) }{\partial cc_{ij}} &= d_{ji},\\ \cfrac{\partial |\boldsymbol{C}|}{\partial c_{ij}} &= C_{ij}\end{align}
ここに、\(C_{ij}\)は\(c_{ij}\)の余因子である。\(\boldsymbol{C} = \boldsymbol{\Sigma}^{-1}\)の各要素についての微分を考える。\(\mathrm{tr}(\boldsymbol{CD})\)の線形性と\(\log|\boldsymbol{C}|\)が凹関数であることより、\(f(\boldsymbol{C}; \boldsymbol{D})\)は\(\boldsymbol{C}\)は\(\boldsymbol{C}\)について凹関数である。\(\boldsymbol{C}\)が正定値行列ではなくなるにつれて、\(f(\boldsymbol{C}; \boldsymbol{D}) \to \infty\)となる。すなわち、\(\boldsymbol{C}\)の最小固有値が\(0\)に近づくにつれて\(\log|\boldsymbol{C}|\to \infty\)である。したがって\(\boldsymbol{C}\)は正定値行列であるので、\(\boldsymbol{C}\)に関して\(f(\boldsymbol{C}; \boldsymbol{D})\)の最大値が一意に存在することがいえる。
\(f(\boldsymbol{C}; \boldsymbol{D})\)の\(\boldsymbol{C} = \boldsymbol{\Sigma}^{-1}\)に関する微分は
\begin{align}\cfrac{\partial f(\boldsymbol{C}; \boldsymbol{D})}{\partial c_{ij}} &= \cfrac{\partial }{\partial c_{ij}}\left[ N\log|\boldsymbol{C}| - \mathrm{tr}(\boldsymbol{CD})\right]\\&= N|\boldsymbol{C}|^{-1} \cfrac{\partial |\boldsymbol{C}|}{dc_{ij}} -d_{ji}\\&= N|\boldsymbol{C}|^{-1} C_{ij}-d_{ji}\end{align}
である。これを\(0\)とおくと次を得る。
\begin{align}&N|\boldsymbol{C}|^{-1} C_{ij}-d_{ji} =0\\ &\Leftrightarrow N|\cfrac{C_{ij}}{|\boldsymbol{C}|} = d_{ji}\\ &\Leftrightarrow c^{ij} = \cfrac{d_{ij}}{N},\end{align}
ここに、\(c^{ij}\)は\(\boldsymbol{C}^{-1}\)の\(i, j\)成分である。よって、各\(i, j\)成分について上式のときに関数\(f(\boldsymbol{C}; \boldsymbol{D})\)は最大値を取る。したがって
\begin{align}&\boldsymbol{C}^{-1} = \cfrac{1}{N}\boldsymbol{D}\\ &\Leftrightarrow \boldsymbol{C} = N\boldsymbol{D} \end{align}
のときに、\(f(\boldsymbol{C}; \boldsymbol{D})\)は最大となる。また、その最大値は
\begin{align}f(N\boldsymbol{D}^{-1}; \boldsymbol{D}) &= N\log|N\boldsymbol{D}^{-1}| -\mathrm{tr}(\boldsymbol{N\boldsymbol{D}^{-1}|D})\\ &= - N\log|\boldsymbol{D}| + N\log |N\boldsymbol{I}| - pN\\&= - N\log|\boldsymbol{D}| + pN\log N- pN\end{align}
である。
