分散共分散行列の計算法

$p$次元のデータが$n$個ある。これを$n \times p$の行列$X$で表す。このとき分散共分散行列は以下に計算できる。$1_n$は全ての成分が$1$の$n$次元縦ベクトルである。

\begin{align}
&S = (1/(n-1))X'J_n X \tag{1} \\
&J_n = I_n - (1/n)1_n (1_n)' \tag{2}
\end{align}

ここで以下の性質がある。(3)式は平均偏差にほかならない。

\begin{align}
&(X_{.j}- 1_n({1/n})\sum_i X_{ij})=J_n X_{.j} \tag{3}\\
&(J_n)^2=J_n  \tag{4}\\
&J_n' = J_n  \tag{5}
\end{align}

以上を、matlabで確認した例を以下に示す。　

N = 5;
P = 3;
X = randn( N, P );
one = ones( N, 1 );
J = eye(N) - 1/N*(one*one');
fprintf( ' X^T * J * X ');
(X'* J * X)/(N-1)
fprintf( ' cov( X ) ')
cov( X )

fprintf( ' Jを二乗しても値は変わらない　J - J*J = 0 \n');
J - J*J

以下は計算例です。

X^T * J * X 
ans =

    0.8465   -0.2367    0.3679
   -0.2367    0.2478   -0.1620
    0.3679   -0.1620    1.0746

 cov( X ) 
ans =

    0.8465   -0.2367    0.3679
   -0.2367    0.2478   -0.1620
    0.3679   -0.1620    1.0746

 Jを二乗しても値は変わらない　J - J*J = 0 

ans =

   1.0e-15 *

   -0.2220         0         0    0.0278    0.0278
         0   -0.2220         0    0.0278    0.0278
         0         0   -0.2220    0.0278    0.0278
    0.0278    0.0278    0.0278   -0.1110    0.0278
    0.0278    0.0278    0.0278    0.0278   -0.1110