/*************************************************************/
/*                   主座標分析 : PCOORD                     */
/*                                                           */
/* DATE: 1997/07/19                                          */
/* CODE: Tokuhisa SUZUKI                                     */
/*************************************************************/

proc iml ;

/*-------------------------------------------------------------
                類似度行列 E （対称行列）の指定
-------------------------------------------------------------*/

    e = {  0  -1  -5 -17 -20 -25 -13  -9 ,
          -1   0  -4 -16 -25 -32 -20 -16 ,
          -5  -4   0  -4 -13 -20 -16 -20 ,
         -17 -16  -4   0  -9 -16 -20 -32 ,
         -20 -25 -13  -9   0  -1  -5 -17 ,
         -25 -32 -20 -16  -1   0  -4 -16 ,
         -13 -20 -16 -20  -5  -4   0  -4 ,
         -9  -16 -20 -32 -17 -16  -4   0  
        }
    ;

/*-----------------------------------------------------------*/

    n = ncol( e )   ;                          /*  対象の数  */
    t_mean = e[ : ] ;                          /*  全体平均  */
    c_mean = e[ :,] ;                          /*  列の平均  */
    r_mean = e[,: ] ;                          /*  行の平均  */

                      /* 行列 E を２重中心化した行列 A 作成  */
    a = e - repeat(r_mean, 1, n) - repeat(c_mean, n, 1) + t_mean ;

    call eigen( eigen, evec, a )    ; /* 行列 A の固有値問題 */
    evec = evec` # sqrt( abs(eigen)); /* ノルム１を固有値に  */
    x    = t( evec )                ; /* 固有ベクトルの転置  */
    proport = eigen / trace( a )    ; /* 寄与率              */
    cumulat = cusum( proport  )     ; /* 累積寄与率          */

    reset noname ;
    print   "類似度行列：Ｅ"     ,, e
          / "２重中心化行列：Ａ" ,, a ;
    print / eigen  [ format = 10.4 colname = '固有値'     ]
            proport[ format =  9.4 colname = '寄与率'     ]
            cumulat[ format =  9.4 colname = '累積寄与率' ] ;
    print / '座標行列：Ｘ',,  x[ format = 9.4 ] ;

    create pcoord from x ; /* 座標行列を出力 */
    append from x        ;
    close pcoord         ;

    /*------------------------ 性質１ -----------------------*/
    x2 = x[ , 1:2 ] ; /* 座標行列の第1，2列だけを抽出 */
    xx = x2 * x2`   ;
    print / '性質１：２次元・座標行列から２重中心化行列の再生',,
            '座標行列（２次元）:X2' x2[ format = 9.4 ] ,
            '再生行列( x2 * x2` )'  xx[ format = 9.0 ] ;

    /*------------------------ 性質２ -----------------------*/
    u = j( n, 1, 1 ) ;
    d = diag( xx ) * u*u`  -  2 # xx  +  u*u` * diag( xx ) ;
    e2 = -2 # e ;
    print / '性質２：主座標空間の点間平方距離は類似度行列の−２倍',,
            '各サンプルの（２次元）点間平方距離行列' d ,
            '類似度行列の各要素の−２倍 ' e2 ;

    /*------------------------ 問１ -----------------------------*/
    mean = x[ :, ] ;
    print / '問１：主座標空間は各次元につき平均０に中心化',,
            '座標行列の各列の平均ベクトル' mean[ format = 4.0 ] ;

    /*------------------------ 問２ -----------------------------*/
    call eigen( eigen_e, evec_e, e ) ;    /* 行列 E の固有値問題 */
    evec_e = evec_e` # sqrt( abs(eigen_e)); /* ノルムを固有値に  */
    y    = t( evec_e ) ;                  /* 固有ベクトルの転置  */
    y2 = y[ , 1:2 ]    ;         /* 座標行列の第1，2列だけを抽出 */
    yy = y2 * y2`      ;
    dy= diag( yy ) * u*u`  -  2 # yy  +  u*u` * diag( yy ) ;
    print / '問２：Ｅを直接スペクトル分解しても点間距離は同じ',,
            'Ｅをスペクトル分解した座標行列' y[ format = 9.4 ] ,
            'Ｅをスペクトル分解した（２次元）点間平方距離行列' dy ;

quit ;
    /*------------------------ 問３ -----------------------------*/
proc princomp data = pcoord out = prin n = 2 ;
    var col1 - col2 ;
run ;
data prin ;
    set prin( keep = col1 col2 prin1 prin2 ) ;
    ratio1 = col1 / prin1 ;
    ratio2 = col2 / prin2 ;
run ;
title '問３：座標行列Ｘを主成分分析すると主成分得点はＸに比例' ;
proc print data = prin ;
    var col1 col2 prin1 prin2 ratio1 ratio2 ;
run ;
title ;

/*---------- 主座標分析による第１軸×第２軸の布置 -----------*/
data label ;
    set pcoord ;
    text     = left( put( _n_, 2. ) ) ;
    x        = col1       ;
    y        = col2       ;
    xsys     = '2'        ;
    ysys     = '2'        ;
    size     = 1          ;
    style    = 'SWISS' ;
    function = 'LABEL'    ;
run ;
/*******
goptions hsize   = 17 cm  vsize   = 20 cm
         horigin =  2 cm  vorigin =  3 cm
         device = winprtm
;
******/


proc gplot data = pcoord ;
    plot col2 * col1 /
         anno   = label frame
         vref   = 0     href   = 0
         lvref  = 3     lhref  = 3
         vaxis  = axis1 haxis  = axis1
         vminor = 1     hminor = 1
    ;
    axis1  length = 6 cm
           offset = ( 2       )
           value  = ( h = 0.5 f = swiss )
           label  = ( h = 0.5 f = swiss )
    ;
    symbol v = none ;
    footnote1 f = mincho h = 1 '主座標分析' ;
run ;
quit ;
footnote1 ;


/**********************************************************
 model-2

    e = {
           0    -1.00 -2.25 -4.10 -4.45 -5.00 -3.60 -3.00 ,
          -1.00  0    -2.00 -4.00 -5.00 -5.65 -4.45 -4.00 ,
          -2.25 -2.00  0    -2.00 -3.60 -4.45 -4.00 -4.45 ,
          -4.10 -4.00 -2.00  0    -3.00 -4.00 -4.45 -5.65 ,
          -4.45 -5.00 -3.60 -3.00  0    -1.00 -2.25 -4.10 ,
          -5.00 -5.65 -4.45 -4.00 -1.00  0    -2.00 -4.00 ,
          -3.60 -4.45 -4.00 -4.45 -2.25 -2.00  0    -2.00 ,
          -3.00 -4.00 -4.45 -5.65 -4.10 -4.00 -2.00  0
        }
    ;


 model-3

    e = {
          0      0.950  0.445 -0.315 -1.000 -0.950 -0.445  0.315,
          0.950  0      0.705  0.000 -0.950 -1.000 -0.705  0.000,
          0.445  0.705  0      0.705 -0.445 -0.705 -1.000 -0.705,
         -0.315  0.000  0.705  0      0.315  0.000 -0.705 -1.000,
         -1.000 -0.950 -0.445  0.315  0      0.950  0.445 -0.315,
         -0.950 -1.000 -0.705  0.000  0.950  0      0.705  0.000,
         -0.445 -0.705 -1.000 -0.705  0.445  0.705  0      0.705,
          0.315  0.000 -0.705 -1.000 -0.315  0.000  0.705  0
        }
    ;

 model-4

    e = {
           0 -1 -3 -5 -6 -7 -5 -3,
          -1  0 -2 -4 -7 -8 -6 -4,
          -3 -2  0 -2 -5 -6 -4 -6,
          -5 -4 -2  0 -3 -4 -6 -8,
          -6 -7 -5 -3  0 -1 -3 -5,
          -7 -8 -6 -4 -1  0 -2 -4,
          -5 -6 -4 -6 -3 -2  0 -2,
          -3 -4 -6 -8 -5 -4 -2  0
        }
    ;
***********************************************************/