/*************************************************************/ /* 主座標分析 : 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 を2重中心化した行列 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)); /* ノルム1を固有値に */ x = t( evec ) ; /* 固有ベクトルの転置 */ proport = eigen / trace( a ) ; /* 寄与率 */ cumulat = cusum( proport ) ; /* 累積寄与率 */ reset noname ; print "類似度行列:E" ,, e / "2重中心化行列:A" ,, a ; print / eigen [ format = 10.4 colname = '固有値' ] proport[ format = 9.4 colname = '寄与率' ] cumulat[ format = 9.4 colname = '累積寄与率' ] ; print / '座標行列:X',, x[ format = 9.4 ] ; create pcoord from x ; /* 座標行列を出力 */ append from x ; close pcoord ; /*------------------------ 性質1 -----------------------*/ x2 = x[ , 1:2 ] ; /* 座標行列の第1,2列だけを抽出 */ xx = x2 * x2` ; print / '性質1:2次元・座標行列から2重中心化行列の再生',, '座標行列(2次元):X2' x2[ format = 9.4 ] , '再生行列( x2 * x2` )' xx[ format = 9.0 ] ; /*------------------------ 性質2 -----------------------*/ u = j( n, 1, 1 ) ; d = diag( xx ) * u*u` - 2 # xx + u*u` * diag( xx ) ; e2 = -2 # e ; print / '性質2:主座標空間の点間平方距離は類似度行列の−2倍',, '各サンプルの(2次元)点間平方距離行列' d , '類似度行列の各要素の−2倍 ' e2 ; /*------------------------ 問1 -----------------------------*/ mean = x[ :, ] ; print / '問1:主座標空間は各次元につき平均0に中心化',, '座標行列の各列の平均ベクトル' mean[ format = 4.0 ] ; /*------------------------ 問2 -----------------------------*/ 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 / '問2:Eを直接スペクトル分解しても点間距離は同じ',, 'Eをスペクトル分解した座標行列' y[ format = 9.4 ] , 'Eをスペクトル分解した(2次元)点間平方距離行列' dy ; quit ; /*------------------------ 問3 -----------------------------*/ 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 '問3:座標行列Xを主成分分析すると主成分得点はXに比例' ; proc print data = prin ; var col1 col2 prin1 prin2 ratio1 ratio2 ; run ; title ; /*---------- 主座標分析による第1軸×第2軸の布置 -----------*/ 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 } ; ***********************************************************/