/*=======================================================*/ /* 数量化V類[コレスポンデンス分析] */ /* SAS/IMLによるアルゴリズムの確認 */ /* Greenacre (1984) 参照 */ /*=======================================================*/ proc iml ; X = { 4 2 3 2 , 4 3 7 4 , 25 10 12 4 , 18 24 33 13 , 10 6 7 2 } ; /* データ行列(頻度行列)X */ n = X[ + ] ; /* Xの要素総計・スカラーn */ P = ( 1/n ) * X ; /* コレスポンデンス行列P */ print p ; nr = nrow( X ) ; /* 行の要素の数をnr */ nc = ncol( X ) ; /* 列の要素の数をnc */ k = min( nr, nc ) ; /* 行列Xのサイズ k */ r = P[ ,+ ] ; /* 行和(周辺比率)r */ c = P[ +, ] ; /* 列和(周辺比率)c */ print r c ; Dr = diag( r ) ; /* rの対角行列Dr */ Dc = diag( c ) ; /* cの対角行列Dc */ Rp = inv( Dr ) * P ; /* 行プロファイル */ Cp = inv( Dc ) * P` ; /* 列プロファイル */ *--------------- <Pの一般化特異値分解> --------------; * P と Ω(Dr) および Φ(Dc) が与えられている場合 ; * 以下の特異値分解をすることができる. ; * Y = Ω1/2 * P * Φ1/2 = UΔV` ただし U`U = V`V = I ; * この時,F = Ω-1/2 * U, G = Φ-1/2 * V と置いて,上式 ; * をPについて解くと, ; * P = FΔG` (ただし F`ΩF = G`ΦG = I ) となる.これ ; * をPについての一般化特異値分解という.すなわち,Pを ; * Yのように変換した時の,Yの通常の特異値分解になる. ; *-------------------------------------------------------; Y = sqrt( Dr ) * P * sqrt( Dc ) ; /* PをYに変換 */ call svd( U, L, V, Y ) ; /* Yの特異値分解 */ print y ; F = inv( sqrt( Dr ) ) * U ; /* 左一般化特異ベクトル */ G = inv( sqrt( Dc ) ) * V ; /* 右一般化特異ベクトル */ Fiden = F` * Dr * F ; /* 単位行列になる確認 */ Giden = G` * Dc * G ; /* 単位行列になる確認 */ Z = F * diag( L ) * G` ; /* ZがPに一致の確認 */ *-------------------------------------------------------; * Pの一般化特異値分解は,Pを変換したYの通常の特異値 ; * 分解である.数量化V類(コレスポンデンス分析)では, ; * Y=Ω-1/2 * P * Φ-1/2 と変換する.Yは独立性からの ; * 偏差行列である.PをP-rc で置換えると,最後の特異値,; * 特異ベクトルが自明解(trivial solutions)となり好都合. ; * A = Ω1/2 * U, B = Φ1/2 * V と置くと, ; * P = AΔB` (ただし,A`Ω-1A = B`Φ-1B = I)となる. ; *-------------------------------------------------------; Y = inv(sqrt( Dr )) * ( P - r*c ) * inv(sqrt( Dc )) ; call svd( U, L, V, Y ) ; /* 変換したYの特異値分解 */ print y ; A = sqrt( Dr ) * U ; /* 左一般化特異ベクトル */ B = sqrt( Dc ) * V ; /* 右一般化特異ベクトル */ Aiden = A` * inv( Dr ) * A ; /* 単位行列になる確認 */ Biden = B` * inv( Dc ) * B ; /* 単位行列になる確認 */ Z = A * diag( L ) * B` + r*c ; /* ZがPに一致の確認 */ Du = diag( L ) ; /* 特異値を要素に持つ対角行列Du */ row = inv( Dr ) * A * Du ; /* 行座標を基準化する */ col = inv( Dc ) * B * Du ; /* 列座標を基準化する */ print / L [ colname = '特異値' format = 6.3 ] row [ colname = '行座標' format = 6.3 ] col [ colname = '列座標' format = 6.3 ] ; run ;