/*=======================================================*/
/*          数量化�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 }     ; /* データ行列（頻度行列）Ｘ */

    n  = X[ + ]            ; /* Ｘの要素総計・スカラーｎ */
    P  = ( 1/n ) * X       ; /* コレスポンデンス行列Ｐ   */
print p ;
    nr = nrow( X )         ; /* 行の要素の数をｎr        */
    nc = ncol( X )         ; /* 列の要素の数をｎc        */
    k  = min( nr, nc )     ; /* 行列Ｘのサイズ k         */

    r  = P[ ,+ ]           ; /* 行和（周辺比率）ｒ       */
    c  = P[ +, ]           ; /* 列和（周辺比率）ｃ       */
print r c ;
    Dr = diag( r )         ; /* ｒの対角行列Ｄｒ         */
    Dc = diag( c )         ; /* ｃの対角行列Ｄｃ         */

    Rp = inv( Dr ) * P     ; /* 行プロファイル           */
    Cp = inv( Dc ) * P`    ; /* 列プロファイル           */

  *--------------- ＜Ｐの一般化特異値分解＞ --------------;
  * Ｐ と Ω（Ｄr） および Φ（Ｄc） が与えられている場合 ;
  * 以下の特異値分解をすることができる．                  ;
  * Ｙ = Ω1/2 * Ｐ * Φ1/2 = UΔV`  ただし U`U = V`V = I ;
  * この時，F = Ω-1/2 * U, G = Φ-1/2 * V と置いて，上式 ;
  * をＰについて解くと，                                  ;
  *  P = FΔG`  （ただし F`ΩF = G`ΦG = I ） となる.これ ;
  * をＰについての一般化特異値分解という．すなわち，Ｐを  ;
  * Ｙのように変換した時の，Ｙの通常の特異値分解になる．  ;
  *-------------------------------------------------------;

    Y = sqrt( Dr ) * P * sqrt( Dc ) ; /* ＰをＹに変換    */
    call svd( U, L, V, 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`    ; /* ＺがＰに一致の確認    */

  *-------------------------------------------------------;
  * Ｐの一般化特異値分解は，Ｐを変換したＹの通常の特異値  ;
  * 分解である．数量化�V類（コレスポンデンス分析）では，  ;
  * Ｙ＝Ω-1/2 * Ｐ * Φ-1/2 と変換する．Ｙは独立性からの ;
  * 偏差行列である．ＰをＰ-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 ) ;   /* 変換したＹの特異値分解 */
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 ; /* ＺがＰに一致の確認 */

    Du  = diag( L ) ;   /* 特異値を要素に持つ対角行列Ｄu */
    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 ;