多変量正規分布のパラメータである平均ベクトルと共分散行列の最尤推定量を導出する。
単変量のときと同様に、平均ベクトルの最尤推定量は標本平均ベクトルとなり、共分散行列の差烏有推定量は標本共分散行列となることを示す。
導出する中で、単変量の標本平均と標本分散を多変量に拡張したものであることもみていく。
平均ベクトルと共分散行列の最尤推定量
\(p\)変量正規分布からの観測ベクトルが与えられたとき、平均ベクトル\(\boldsymbol{\mu}\)と共分散行列\(\boldsymbol{\Sigma}\)の推定量が知りたい。そこで最尤推定量を用いる。
多変量正規分布の下で、様々な推定や検定問題をするときに最尤推定量が非常に有用である。最尤推定量は最適な性質を持つことが知られている。
\(N(\boldsymbol{\mu}, \boldsymbol{\Sigma})\)に従う確率ベクトル\(\boldsymbol{X}\)の\(N\)個の観測ベクトルを\(\boldsymbol{x}_1, \ldots, \boldsymbol{x}_N\)とする。ここに\(N>p\)とする。このとき、尤度関数は
\begin{align}L &= \prod_{\alpha = 1}^N\cfrac{1}{(2\pi)^{\frac{1}{2}p} |\boldsymbol{\Sigma}|^{\frac{1}{2}}} \exp\left[-\cfrac{1}{2}(\boldsymbol{x}_{\alpha} - \boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}_{\alpha} - \boldsymbol{\mu})\right]\\ \label{eq1}&= \cfrac{1}{(2\pi)^{\frac{1}{2}pN} |\boldsymbol{\Sigma}|^{\frac{1}{2}}} \exp\left[-\cfrac{1}{2}\sum_{\alpha=1}^N (\boldsymbol{x}_{\alpha} - \boldsymbol{\mu})^T\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}_{\alpha} - \boldsymbol{\mu})\right]\tag{1}\end{align}
である。尤度関数はベクトル\(\boldsymbol{x}_1, \ldots, \boldsymbol{x}_N\)が与えられたときの\(\boldsymbol{\mu}\)と\(\boldsymbol{\Sigma}\)の関数である。\(\boldsymbol{\mu}\)と\(\boldsymbol{\Sigma}\)の関数であることを強調するために、次のように\(\boldsymbol{\mu}\)と\(\boldsymbol{\Sigma}\)を\(\boldsymbol{\mu}^*\)と\(\boldsymbol{\Sigma}^*\)で置き換える。
\begin{align}\label{eq2} L= \cfrac{1}{(2\pi)^{\frac{1}{2}pN} |\boldsymbol{\Sigma}^*|^{\frac{1}{2}}} \exp\left[-\cfrac{1}{2}\sum_{\alpha=1}^N (\boldsymbol{x}_{\alpha} - \boldsymbol{\mu}^*)^T\boldsymbol{\Sigma}^{*-1}(\boldsymbol{x}_{\alpha} - \boldsymbol{\mu}^*)\right]\tag{2}\end{align}
また、対数尤度関数は
\begin{align}\log L &= \log \left\{\cfrac{1}{(2\pi)^{\frac{1}{2}pN} |\boldsymbol{\Sigma}^*|^{\frac{1}{2}}} \exp\left[-\cfrac{1}{2}\sum_{\alpha=1}^N (\boldsymbol{x}_{\alpha} - \boldsymbol{\mu}^*)^T\boldsymbol{\Sigma}^{*-1}(\boldsymbol{x}_{\alpha} - \boldsymbol{\mu}^*)\right]\right\}\\ \label{eq3} &= -\cfrac{1}{2}pN\log 2\pi - \cfrac{1}{2}N\log |\boldsymbol{\Sigma}^*| - \cfrac{1}{2}\sum_{\alpha=1}^N (\boldsymbol{x}_{\alpha} - \boldsymbol{\mu}^*)^T\boldsymbol{\Sigma}^{*-1}(\boldsymbol{x}_{\alpha} - \boldsymbol{\mu}^*)\tag{3} \end{align}
となる。\(\log L \)は\(L\)の単調増加関数であるため、その最大値は\(\boldsymbol{\mu}^*\)、\(\boldsymbol{\Sigma}^*\)の空間における最大値と\(L\)の最大値は同じ点で起きる。\(\boldsymbol{\mu}\)、\(\boldsymbol{\Sigma}\)の最尤推定量は\(\log L\)を最大化させるベクトル\(\boldsymbol{\mu}^*\)、正定値行列\(\boldsymbol{\Sigma}^*\)であることが分かる。
ここで、標本平均ベクトルを次で定義する。
\begin{align}\label{eq4} \bar{\boldsymbol{x}} &= \cfrac{1}{N}\sum_{\alpha=1}^N \boldsymbol{x}_{\alpha} = \begin{pmatrix}\cfrac{1}{N}\sum_{\alpha=1}^Nx_{1\alpha}\\ \vdots \\\cfrac{1}{N}\sum_{\alpha=1}^Nx_{p\alpha}\end{pmatrix} = \begin{pmatrix}\bar{x}_{1}\\\vdots\\\bar{x}_p\end{pmatrix},\tag{4} \end{align}
ここに、\(\boldsymbol{x}_{\alpha} = (x_{1\alpha}, \ldots, x_{p\alpha})^T\)、\(\bar{x}_i = \sum_{\alpha=1}^Nx_{i\alpha}/N\)である。平均との差の2乗とクロス積から成る行列\(\boldsymbol{A}\)を次で定義する。
\begin{align}\sum_{\alpha =1}^N(\boldsymbol{x}_{\alpha} -\bar{\boldsymbol{x}})(\boldsymbol{x}_{\alpha} -\bar{\boldsymbol{x}})^T &= \sum_{\alpha= 1}^N\begin{pmatrix}x_{1\alpha} - \bar{x}_1\\ \vdots \\ x_{p\alpha} - \bar{x}_p\end{pmatrix} \begin{pmatrix} x_{1\alpha} - \bar{x}_1 &\cdots & x_{p\alpha} - \bar{x}_p\end{pmatrix}\\ &= \begin{pmatrix} \sum_{\alpha=1}^N(x_{1\alpha} - \bar{x}_1)^2 &\cdots & \sum_{\alpha=1}^N(x_{1\alpha} - \bar{x}_1)(x_{p\alpha} - \bar{x}_p) \\ \vdots & & \vdots\\ \sum_{\alpha=1}^N(x_{p\alpha} - \bar{x}_p)(x_{1\alpha} - \bar{x}_1) & \cdots & \sum_{\alpha=1}^N(x_{p\alpha} - \bar{x}_p)^2\end{pmatrix} \\&=\label{eq5} \left[\sum_{\alpha = 1}^N(x_{i\alpha} - \bar{x}_i)(x_{j\alpha} - \bar{x}_j)\right] ,\ \ \ \ i, j = 1,\ldots, p.\tag{5}\end{align}
この行列\(\boldsymbol{A}\)を変形するために、次の補題を用いる。
補題１
\(\boldsymbol{x}_1, \ldots, \boldsymbol{x}_N\)を\(N\)個のベクトルとし、\(\bar{\boldsymbol{x}}\)は\eqref{eq3}で定義したベクトルである。このとき任意のベクトル\(\boldsymbol{b}\)に対して、次が言える。
\begin{align}\label{eq6} \sum_{\alpha = 1}^N(\boldsymbol{x}_{\alpha} - \boldsymbol{b})(\boldsymbol{x}_{\alpha} - \boldsymbol{b})^T = \sum_{\alpha=1}^N(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})^T + (\bar{\boldsymbol{x}} - \boldsymbol{b})(\bar{\boldsymbol{x}} - \boldsymbol{b})^T.\tag{6}\end{align}
証明
\begin{align} \sum_{\alpha = 1}^N(\boldsymbol{x}_{\alpha} - \boldsymbol{b})(\boldsymbol{x}_{\alpha} - \boldsymbol{b})^T&= \sum_{\alpha=1}^N \left[(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})+ (\bar{\boldsymbol{x}} - \boldsymbol{b})\right]\left[(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})+ (\bar{\boldsymbol{x}} - \boldsymbol{b})\right]^T\\ &= \sum_{\alpha = 1}^N \Bigl[(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})^T + (\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})(\bar{\boldsymbol{x}} - \boldsymbol{b})^T\\&\ \ \ \ + (\bar{\boldsymbol{x}} - \boldsymbol{b})(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})^T+(\bar{\boldsymbol{x}} - \boldsymbol{b})(\bar{\boldsymbol{x}} - \boldsymbol{b})^T\Bigr]\\ &= \sum_{\alpha=1}^N(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})^T + N(\bar{\boldsymbol{x}} - \bar{\boldsymbol{x}})(\bar{\boldsymbol{x}}-\boldsymbol{b})^T\\ &\ \ \ \ + N(\bar{\boldsymbol{x}} - \boldsymbol{b})(\bar{\boldsymbol{x}} - \bar{\boldsymbol{x}}) + N(\bar{\boldsymbol{x}} - \boldsymbol{b})(\bar{\boldsymbol{x}} - \boldsymbol{b})^T\\&= \sum_{\alpha=1}^N(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})^T + (\bar{\boldsymbol{x}} - \boldsymbol{b})(\bar{\boldsymbol{x}} - \boldsymbol{b})^T.\ \ \ \ □ \end{align}
補題1の\(\boldsymbol{b}\)に\(\boldsymbol{\mu}^*\)を適用することで、次を得る。
\begin{align}\sum_{\alpha=1}^N(\boldsymbol{x}_{\alpha} - \boldsymbol{\mu}^*)(\boldsymbol{x}_{\alpha} - \boldsymbol{\mu}^*)^T &= \sum_{\alpha=1}^N(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})^T + N(\bar{\boldsymbol{x}} - \boldsymbol{\mu}^*)(\bar{\boldsymbol{x}} - \boldsymbol{\mu}^*)^T\\ \label{eq7} &= \boldsymbol{A} + N(\bar{\boldsymbol{x}} - \boldsymbol{\mu}^*)(\bar{\boldsymbol{x}} - \boldsymbol{\mu}^*)^T.\tag{7}\end{align}
\eqref{eq7}と二次形式をトレースで表現することによって、次の関係を得る。
\begin{align}\sum_{\alpha=1}^N(\boldsymbol{x}_{\alpha} - \boldsymbol{\mu}^*)^T\boldsymbol{\Sigma}^{*-1}(\boldsymbol{x}_{\alpha} - \boldsymbol{\mu}^*) &= \mathrm{tr}\sum_{\alpha=1}^N\boldsymbol{\Sigma}^{*-1}(\boldsymbol{x}_{\alpha} - \boldsymbol{\mu}^*)(\boldsymbol{x}_{\alpha} - \boldsymbol{\mu}^*)^T\\ &= \mathrm{tr}\boldsymbol{\Sigma}^{*-1}\sum_{\alpha=1}^N(\boldsymbol{x}_{\alpha} - \boldsymbol{\mu}^*)(\boldsymbol{x}_{\alpha} - \boldsymbol{\mu}^*)^T\\ &= \mathrm{tr}\boldsymbol{\Sigma}^{*-1}\left[\boldsymbol{A} + N(\bar{\boldsymbol{x}} - \boldsymbol{\mu}^*)(\bar{\boldsymbol{x}} - \boldsymbol{\mu}^*)^T\right]\\ &= \mathrm{tr}\boldsymbol{\Sigma}^{*-1}\boldsymbol{A} + \mathrm{tr}\boldsymbol{\Sigma}^{*-1}N(\bar{\boldsymbol{x}} - \boldsymbol{\mu}^*)(\bar{\boldsymbol{x}} - \boldsymbol{\mu}^*)^T \\ \label{eq8}&= \mathrm{tr}\boldsymbol{\Sigma}^{*-1}\boldsymbol{A} + N(\bar{\boldsymbol{x}} - \boldsymbol{\mu}^*)^T\boldsymbol{\Sigma}^{*-1}(\bar{\boldsymbol{x}} - \boldsymbol{\mu}^*).\tag{8} \end{align}
このことから、\eqref{eq2}の対数尤度関数は次の書き換えることができる。
\begin{align}\label{eq9} \log L &= -\cfrac{1}{2}pN \log 2\pi - \cfrac{1}{2}N\log |\boldsymbol{\Sigma}^*| -\cfrac{1}{2}\mathrm{tr}\boldsymbol{\Sigma}^{*-1}\boldsymbol{A} - \cfrac{1}{2}N(\bar{\boldsymbol{x}} - \boldsymbol{\mu}^*)^T\boldsymbol{\Sigma}^{*-1}(\bar{\boldsymbol{x}} - \boldsymbol{\mu}^*).\tag{9}\end{align}
\(\boldsymbol{\Sigma}^*\)は正定値行列であるため、\(\boldsymbol{\Sigma}^{*-1}\)も正定値行列である。さらに、\(N(\bar{\boldsymbol{x}} - \boldsymbol{\mu}^*)^T\boldsymbol{\Sigma}^{*-1}(\bar{\boldsymbol{x}} - \boldsymbol{\mu}^*) \geq 0\)であり、\(\boldsymbol{\mu}^* = \bar{\boldsymbol{x}}\)のときのみ\(0\)をとる。\(\boldsymbol{\Sigma}^*\)の正定値性から\eqref{eq9}は負であるため、\(\boldsymbol{\mu}^*\)に関して\(\boldsymbol{\mu}^* = \bar{\boldsymbol{x}}\)のとき、\eqref{eq9}は最大となる。つまり、\(\boldsymbol{\mu}\)の最尤推定量は\(\bar{\boldsymbol{x}}\)である。次に、\(\boldsymbol{\Sigma}^*\)に関しての最大化を考える。\eqref{eq9}の第2、第3項を最大化させるために、次の補題を用いる。
補題２
\(\boldsymbol{D}\)が\(p\)次の正定値行列であるとき
\begin{align}\label{eq10} f(\boldsymbol{G}) =-N\log |\boldsymbol{G}| - \mathrm{tr}\boldsymbol{F}^{-1}\boldsymbol{D}\tag{10}\end{align}
の最大値は、正定値行列\(\boldsymbol{G}\)に関して存在し、\(\boldsymbol{G} = (1/N)\boldsymbol{D}\)のとき最大となり、次の最大値をもつ。
\begin{align}\label{eq11} f\bigl((1/N)\boldsymbol{D}\bigr) = pN\log N -N\log |\boldsymbol{D}| -pN.\tag{11}\end{align}
証明 \(\boldsymbol{D} = \boldsymbol{EE}^T\)、\(\boldsymbol{E}^T\boldsymbol{G}^{-1}\boldsymbol{E} = \boldsymbol{H}\)とする。このとき、
\begin{align}& \boldsymbol{E}^T\boldsymbol{G}^{-1}\boldsymbol{E} = \boldsymbol{H}\\ &\Leftrightarrow \boldsymbol{G}^{-1} = (\boldsymbol{E}^T)^{-1}\boldsymbol{HE}^{-1}\\ &\Leftrightarrow \boldsymbol{G} = \boldsymbol{EH}^{-1}\boldsymbol{E}^T\end{align}
である。このことから
\begin{align} |\boldsymbol{G} &= \boldsymbol{E}| \cdot |\boldsymbol{H}^{-1}| \cdot |\boldsymbol{E}^T| \\ &= \Leftrightarrow |\boldsymbol{EE}^T|\cdot|\boldsymbol{H}|^{-1} \\ &= \cfrac{|\boldsymbol{D}|}{|\boldsymbol{H}|}\end{align}
であり
\begin{align}\mathrm{tr}(\boldsymbol{G}^{-1}\boldsymbol{D}) &= \mathrm{tr}(\boldsymbol{G}^{-1}\boldsymbol{EE}^T)\\&=\mathrm{tr}(\boldsymbol{E}^T\boldsymbol{G}^{-1}\boldsymbol{E}) \\&= \mathrm{tr}(\boldsymbol{H}) \end{align}
である。\(\boldsymbol{D} = \boldsymbol{EE}^T\)の構成法はスペクトル分解などによって可能である。よって、\eqref{eq10}は次のように表される。
\begin{align}f &= -N\log \cfrac{|\boldsymbol{D}|}{|\boldsymbol{H}|} -\mathrm{tr}(\boldsymbol{H})\\\label{eq12} &= -N\log |\boldsymbol{D}| + N\log |\boldsymbol{H}| -\mathrm{tr}(\boldsymbol{H}).\tag{12}\end{align}
\(\boldsymbol{H} = \boldsymbol{TT}^T\)とする。ここに、\(\boldsymbol{T}\)は下三角行列である。この分解に関しては、コレスキー分解を参照されたい。\(\boldsymbol{T}\)について次がいえる。
\begin{align}|\boldsymbol{T}|^2 &= \begin{vmatrix}t_{11} & 0 & 0 &\cdots & 0\\ t_{21} & t_{22} & 0 & \cdots & 0\\ \vdots & \vdots & \vdots & & \vdots \\t_{p1} & t_{p2} & t_{3p} &\cdots & t_{pp}\end{vmatrix}^2\\&=\left(\prod_{i=1}^pt_{ii}\right)^2\\\label{eq13}&= \prod_{i=1}^pt_{ii}^2.\tag{13} \end{align}
また
\begin{align}\mathrm{tr}(\boldsymbol{TT}^T) &= \mathrm{tr}\left[\begin{matrix}t_{11}^2 & t_{11}t_{21} & \cdots & t_{11}t_{p1} \\ t_{21}t_{11} & t_{21}^2 + t_{22}^2 &\cdots & t_{21}t_{p1} + t_{22}t_{p2} \\ \vdots & \vdots & & \vdots\\ t_{p1}t_{11} & t_{p1}t_{21} + t_{p2}t_{22} & \cdots & t_{p1}^2 + \cdots + t_{pp}^2\end{matrix}\right]\\ &= \sum_{i=1}^pt_{ii}^2 +\sum_{i>j}t_{ij}^2\end{align}
であることから
\begin{align}f &= N\log |\boldsymbol{D}| + N\log|\boldsymbol{TT}^T| -\mathrm{tr}(\boldsymbol{TT}^T)\\&= -N\log|\boldsymbol{D}| + N\log |\boldsymbol{T}|^2 -\mathrm{tr}(\boldsymbol{TT}^T)\\&= -N\log|\boldsymbol{D}| +N\log \prod_{i=1}^pt_{ii}^2 -\left(\sum_{i=1}^pt_{ii}^2 + \sum_{i>j}t_{ij}^2\right) \\&= -N\log|\boldsymbol{D}| + \sum_{i=1}^p(N\log t_{ii}^2 - t_{ii}^2) -\sum_{i>j}t_{ij}^2\end{align}
の最大値は\(t_{ii} = \sqrt{N}\)、\(t_{ij} = 0,\ i\neq j\)で起きる。これは
\begin{align}&\cfrac{\partial f}{\partial t_{ii}} = 0\\ &\Leftrightarrow \cfrac{\partial}{\partial t_{ii}} \left[\sum_{j=1}^p (N\log t_{jj}^2 - t_{jj}^2)\right] \\&\Leftrightarrow N \cfrac{2}{t_{ii}}- 2t_{ii} = 0\\&\Leftrightarrow t_{ii}^2 = N\\&\Rightarrow t_{ii} = \sqrt{N}\end{align}
であることを用いた。\(t_{ii} = \sqrt{N}\)、\(t_{ij} = 0,\ i\neq j\)のとき、\(\boldsymbol{H}\)は次のようになる。
\begin{align}\boldsymbol{H} &= \boldsymbol{TT}^T\\&= \begin{pmatrix}t_{11}^2 & t_{11}t_{21} & \cdots & t_{11}t_{p1} \\ t_{21}t_{11} & t_{21}^2 + t_{22}^2 &\cdots & t_{21}t_{p1} + t_{22}t_{p2} \\ \vdots & \vdots & & \vdots\\ t_{p1}t_{11} & t_{p1}t_{21} + t_{p2}t_{22} & \cdots & t_{p1}^2 + \cdots + t_{pp}^2\end{pmatrix} \\&= \begin{pmatrix}N & 0 & \cdots & 0\\ 0 & N &\cdots & 0\\\vdots &\vdots & & \vdots\\ 0 & 0 &\cdots & N\end{pmatrix}\\ &= N \boldsymbol{I}. \end{align}
よって
\begin{align}\boldsymbol{G} &= \boldsymbol{EH}^{-1}\boldsymbol{E}^T \\&= \cfrac{1}{N}\boldsymbol{EE}^T \\&= \cfrac{1}{N}\boldsymbol{D}\end{align}
のときに最大値をとり、その最大値は
\begin{align} f &= -N\log | \boldsymbol{D}| +N\log |N\boldsymbol{I}| -\mathrm{tr}(N\boldsymbol{I}) \\&= -N\log|\boldsymbol{D}| + N\log N^p -pN\\&= -N\log|\boldsymbol{D}| + pN\log N -pN \end{align}
である。□
ここでは、コレスキー分解を対数尤度関数の最大化に用いたが、スペクトル分解、帰納法、行列・ベクトル微分を用いて最大化することも可能である。これは多変量正規分布の最尤推定量の様々な導出方法で紹介している。
定理１ 平均ベクトルと共分散行列の最尤推定量
平均ベクトルと共分散行列の最尤推定量
\(\boldsymbol{x}_1, \ldots, \boldsymbol{x}_N\)が\(N(\boldsymbol{\mu}, \boldsymbol{\Sigma})\)からの\(p\)観測ベクトルであるとき、\(\boldsymbol{\mu}\)と\(\boldsymbol{\Sigma}\)の最尤推定量\(\hat{\boldsymbol{\mu}}\)と\(\hat{\boldsymbol{\Sigma}}\)はそれぞれ次で与えられる。
\begin{align}\hat{\boldsymbol{\mu}} &= \bar{\boldsymbol{x}} = \cfrac{1}{N}\sum_{\alpha=1}^N \boldsymbol{x}_{\alpha}, \\ \hat{\boldsymbol{\Sigma}} &= \cfrac{1}{N}\sum_{\alpha=1}^N(\boldsymbol{x}_{\alpha} \bar{\boldsymbol{x}})(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})^T.\end{align}
また、補題1に\(\boldsymbol{b} = \boldsymbol{0}\)を適用することで次の関係を得る。
\begin{align}&\sum_{\alpha = 1}^N\boldsymbol{x}_{\alpha}\boldsymbol{x}_{\alpha}^T = \sum_{\alpha = 1}^N (\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})^T + N\bar{\boldsymbol{x}}\bar{\boldsymbol{x}}^T\\&\Leftrightarrow \sum_{\alpha = 1}^N (\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})^T = \sum_{\alpha=1}^N\boldsymbol{x}_{\alpha}\boldsymbol{x}_{\alpha}-N\bar{\boldsymbol{x}}\bar{\boldsymbol{x}}^T\end{align}
したがって、共分散行列\(\boldsymbol{\Sigma}\)の最尤推定量である標本共分散行列\(\hat{\boldsymbol{\Sigma}} = (1/N)\sum_{\alpha = 1}^N(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})(\boldsymbol{x}_{\alpha} - \bar{\boldsymbol{x}})^T\)は次のように表現できる。
\begin{align}\hat{\boldsymbol{\Sigma}} &= \cfrac{1}{N}\sum_{\alpha=1}^N \boldsymbol{x}_{\alpha}\boldsymbol{x}_{\alpha}^T - \bar{\boldsymbol{x}}\bar{\boldsymbol{x}}^T .\end{align}
これは、単変量のときの公式
\((1/N)\sum_{\alpha=1}^N(x_{\alpha} - \bar{x})^2 = (1/N)\sum_{\alpha=1}^N x_{\alpha}^2 - \bar{x}^2\)
を\(p\)次元に拡張したものであることが分かる。
