「統計学関連なんでもあり」の過去ログ---028

★ Re: 一般化固有値問題の解と正準相関分析 ★

2434.　Re: 一般化固有値問題の解と正準相関分析　tes　2004/02/22 (日) 05:20
├2453.　Re^2: 一般化固有値問題の解と正準相関分析　tes　2004/02/24 (火) 00:56
│├2472.　Re^3: 固有値問題の解　uchan　2004/02/24 (火) 20:39
│└2463.　Re^3: 一般化固有値問題の解と正準相関分析　青木繁伸　2004/02/24 (火) 14:40
├2452.　Re^2: 一般化固有値問題　青木繁伸　2004/02/23 (月) 22:42
└2443.　Re^2: 一般化固有値問題の解と正準相関分析　uchan　2004/02/23 (月) 10:38

2434.　Re: 一般化固有値問題の解と正準相関分析　tes　　2004/02/22 (日) 05:20

こんばんは。正準相関分析のプログラムを作っています。
私はfortranしか使えないので，残念ながらuchanさんや青木さんに教えてもらったCのプログラムは読めませんでした。そこで数値計算の本を参考にfortranでパワー法のプログラムを作りました。n次元の相関行列からn個の固有値を計算する事ができました。精度もなかなか高いです。しかし，固有値に対する固有ベクトルが相関変量と一致しませんでした。固有値に対する固有ベクトル(a1,a2,a3・・・an)がそのまま相関変量u=b1*x1+b2*x2・・・bn*xn の(b1,b2,b3・・・bn)になるのではないのでしょうか？

　　　　　[このページのトップへ]

2453.　Re^2: 一般化固有値問題の解と正準相関分析　tes　　2004/02/24 (火) 00:56

データの標準化は済んでおります。標準偏差で割る事により，正規化した相関行列が得られました。
http://gucchi24.hp.infoseek.co.jp/SEIJUNEX.htm
このページを参考にプログラムを作ったのですが，最大固有値は0.936となりこのページと一致しました。しかし固有値から固有ベクトルを求めるアルゴリズムがわかりません。従来の固有値の問題では固有ベクトルUの条件がU'*U=1　(U'はUの転置行列)　であるのに対して，正準相関行列の固有値ベクトルはUの条件は
U'*R11*U=1
であるところがポイントだと思います。私の作ったパワー法のプログラムはU'*U=1なので，結果として違う固有ベクトルが出てくるのだと思います。

ps,QR分解のプログラムも考えてみたのですが，ヘッセンベルグ行列への変換でつまづいています。

　　　　　[このページのトップへ]

2472.　Re^3: 固有値問題の解　uchan　　2004/02/24 (火) 20:39

固有値をgとしたとき

A= R11(-1)*R12*R22(-1)*R21
B= R22(-1)*R21/√g
の関係がありますから

重みベクトルを v1=(a,b) v2=(p,q)とすると

A*v1=g*v1
v2=B*v1
の関係が成り立っているので固有値gが求まっていれば上の関係から重みベクトルv1,v2が求められるはずです

FORTRANなら汎用計算機用の数値計算ライブラリーが通常用意されているのではないかと思うのでソースコードが公開されていたらプログラム自作の参考に出来るかも知れません ^^;

　　　　　[このページのトップへ]

2463.　Re^3: 一般化固有値問題の解と正準相関分析　青木繁伸　　2004/02/24 (火) 14:40

何事によらず，書籍の情報というのは役に立つものです。
村越勝弘「FORTRANのための数値計算法」日本コンピュータ協会，1972年
という古い本を本棚から引っ張り出して，再確認（だいぶ前に見たときの書き込みもあった(^_^;)）。

269ページからの記述とプログラムを参考にアルゴリズムを追いかけ，だんだん思い出してきました。
今回は R という強い見方があるので，ちゃんと理解できているか段階を追って確認できました。R は，こんな風に線形代数の勉強にも使えるんですね(^_^;)

あれこれ解説するのも大変だけど，R で記述すると非常に簡単なので，プログラムを示して解説に代えます。
Ax=λBx, A:実対称行列，B:実対称正値行列，λ:スカラー，x:列行列
一般化固有値問題を解くプログラム
geneig <- function(A, B) { ans <-eigen(B) T <- ans$vectors A2 <- t(T) %*% A %*% T B <- diag(1/sqrt(ans$values)) B2 <- B %*% A2 %*% B ans2 <- eigen(B2) list(values=ans2$values, vectors=T %*% B %*% ans2$vectors) }
戻り値は，固有値と固有ベクトル

これを利用する正準相関分析のプログラム
x <- read.table("cancor.data2") #データ行列を読む r <- cor(x) #相関係数行列を計算する gr1=1:7 #変数群の定義 gr2=8:12 #変数群の定義 S11 <- r[gr1, gr1] # 行列の切り出し S22 <- r[gr2, gr2] S12 <- r[gr1, gr2] geneig(S12 %*% solve(S22) %*% t(S12), S11) # 片方の解 geneig(t(S12) %*% solve(S11) %*% S12, S22) # もう一方の解
階数の小さい方の解を求めて，それを変形して階数の大きい方の解を求めて...なんて書いてあるけど，面倒くさいので上のようにした。

　　　　　[このページのトップへ]

2452.　Re^2: 一般化固有値問題　青木繁伸　　2004/02/23 (月) 22:42

変なタイトルを付けないようにしてくださいね。教えてくださいというような意味のないタイトルより悪い(^_^;)。

R11(-1)・R12・R22(-1)・R21 の固有値は目的の固有値と同じになるようですが，固有ベクトルは違うようですよ。どう違うか，どうすれば正しい固有ベクトルが求められれるか，フォローする余裕がありませんが。

　　　　　[このページのトップへ]

2443.　Re^2: 一般化固有値問題の解と正準相関分析　uchan　　2004/02/23 (月) 10:38

>固有値に対する固有ベクトルが相関変量と一致しませんでした。固有値に対する固有ベクトル(a1,a2,a3・・・an)がそのまま相関変量u=b1*x1+b2*x2・・・bn*xn の(b1,b2,b3・・・bn)になるのではないのでしょうか？

正準相関分析はあまり詳しくはないのですが(大部分の多変量解析の基礎のベースとも考えられるので重要なんでしょうけど)

正準相関分析で重みベクトル(固有ベクトルのこと)を求める前にデータを標準化する必要がありますがしてますでしょうか

はずしていたらごめんなさい

青木先生がフォローしてくれるかも知れません
PS.
パワー法より精度がよい方法としてQR分解を使うという手があります。fortranではポインターに替わるものなかった?ような気がしますのでコーディングが少し手間かも参考まで

　　　　　[このページのトップへ]