共分散行列の求め方

数個のデータから共分散行列を生成するプログラムを実装しようと
考えています。そこで、次のようなプログラムを書きましたが、
結果があっているかがわかりません。共分散行列の生成方法に
詳しい方がおられましたら、どうか教えて頂きたいと思います。
共分散行列生成に用いた式は
Σ = E[(X-E[X])(X-E[X])^T] 　です。
(Σ：共分散行列　X：データ E[X]：Xの平均　^T：転置行列）


double xave = 0.0; //式のE[x]にあたる部分
double work [ ]; //データを入れる一次元配列（式のXにあたる）
double a[][]; //共分散行列を入れる配列

work[0] = 3; //5つのデータを格納
work[1] = -2;
work[2] = 5;
work[3] = -3;
work[4] = 2;

for(i = 0; i < N; i++){　　 //E[x]を求めるためにxaveに全要素合計を入れる
xave += work[i];
}
xave /= N; 　　　　　　　//全要素合計をデータ数で割る
for(i = 0; i < N; i++){ 　　//要素Xから平均を引く(X-E[x]にあたる）
work[i] -= xave;
}
for(i = 0; i < N; i++){
for(j = 0; j < N; j++){
a[i][j] = work[i]*work[j]/N; 　//a[][]に値を格納 Nで割っているのは式の一番外側のEにあたる
}
}

＜実行結果＞

0.0 | 0.0 | 0.0 | 0.0 | 0.0
0.0 | 4.0 | -6.0 | 8.0 | -8.0
0.0 | -6.0 | 9.0 | -12.0 | 12.0
0.0 | 8.0 | -12.0 | 16.0 | -16.0
0.0 | -8.0 | 12.0 | -16.0 | 16.0

No.2 ベストアンサー 回答者： quaestio 回答日時： 2011/10/16 15:39

> ・言語は一応Javaでやっています。

Javaでしたか、ずいぶん前にちょっと使ってみただけでしたので思いつきませんでした。

> Σ = E[(X-E[X])(X-E[X])^T] の式から、(X-E[X])を一列の列ベクトル、
> (X-E[X])^Tを一行の行ベクトルと考え、二つの行列の積を考えるとちょうど
> 正方行列になるので、これでいいかと思っていました。しかし、この考え方
> では間違っているとういことでしょうか。

はい、そうなります。
E[X]は列ベクトルとなりますが、Xの要素が異なる変数なので、それぞれの期待値が等しくなるとは限りません。
なので、

> for(i = 0; i < N; i++){ 　　//要素Xから平均を引く(X-E[x]にあたる）
> work[i] -= xave;
> }

とするのはおかしいのです。
引くなら、work[i]の値を引くべきです。
（データ数が1個なので平均はwork[i]そのもの。その結果共分散行列は5×5の零行列になる）


> それならば、共分散行列を生成するには、二次元以上のデータが必要ということでしょうか。

例えば、下記のようなデータなら計算できます。

番号　　変数１　　変数２　　変数３　　変数４　　変数５
　１　　　　５　　　　　３　　　　　４　　　　　６　　　　１
　２　　　　…
…

上のような形式のデータなら変数ごとの分散と変数間の共分散を求めて、共分散行列を得ることが出来ます。

この回答へのお礼

なるほど、quaestioさんの書いて下さった説明を読み、納得できました。
確かに一次元のデータのみでは共分散行列が求められないですね。
何度も親切に答えて頂き感謝します。ありがとうございました。

お礼日時：2011/10/16 16:08

No.1 回答者： quaestio 回答日時： 2011/10/16 14:06

いろいろと問題がありますが、一番の問題は一次元のデータから共分散行列は計算できません。

aは対角要素がデータの平均からの偏差の二乗で、非対角要素がデータ間の平均からの偏差の積であるだけで、これは共分散行列とは言いません。

その他の問題としては、
・プログラムの質問なのに使用した言語の記載が無い（C/C++だとは想像できますが）
・使用している変数の説明がない（Nの値は？）
・結果表示の関数が記載されていないが、多分間違っている（何故0.0という結果がでるのですか？）

この回答への補足

ご回答ありがとうございます。いろいろと抜け落ちていたり、
間違って記述していました。すみませんでした。

・言語は一応Javaでやっています。
・変数Nは配列の要素数（データの数）です。

実行結果ですが、配列の表示範囲を間違えていました。実際は
4.0 -6.0 8.0 -8.0 2.0
-6.0 9.0 -12.0 12.0 -3.0
8.0 -12.0 16.0 -16.0 4.0
-8.0 12.0 -16.0 16.0 -4.0
2.0 -3.0 4.0 -4.0 1.0
でした。

また、プログラムの最後のa[i][j] = work[i]*work[j]/N;の部分も
正しくはa[i][j] = work[i]*work[j];で /N　はなかったです。

一次元のデータから共分散行列は計算できないということですが、
Σ = E[(X-E[X])(X-E[X])^T] の式から、(X-E[X])を一列の列ベクトル、
(X-E[X])^Tを一行の行ベクトルと考え、二つの行列の積を考えるとちょうど
正方行列になるので、これでいいかと思っていました。しかし、この考え方
では間違っているとういことでしょうか。
それならば、共分散行列を生成するには、二次元以上のデータが必要ということ
でしょうか。

補足日時：2011/10/16 14:33