/*---------------------------------------------*/ /* */ /* SAS/IML で「正準判別分析」の計算方法を確認 */ /* */ /* */ /* AUTHOR: 鈴木督久 */ /* DATE: 1999.3.1 */ /*---------------------------------------------*/ options nosource pageno = 1 ; %macro candisc ( data = , /* 分析データセット名 */ class = , /* 分類変数名 */ var = /* 独立変数名：省略形は不可 */ ) ; /* 平均値データ _mean_ 作成 -------------------------*/ proc summary data = &data ; var &var ; output out = _tmean_ mean = ; run ; proc summary data = &data nway ; class &class ; var &var ; output out = _mean_ mean = ; run ; /* 群数データ groupnum と &groupnum 作成 --------------------------------------*/ data groupnum ; keep groupnum ; set _mean_( obs = 1 ) nobs = n ; groupnum = n ; output ; call symput( 'groupnum', trim(left( put( groupnum, 3. ))) ) ; stop; run ; /* 群別データ __g1 __g2 ... 作成 -------------------------------------*/ proc sort data = &data out = _totdat_ ; by &class ; run ; data __g0 %do i = 1 %to &groupnum ; __g&i %end ; ; keep &class &var ; set _totdat_ ; by &class ; if first.&class then _i_ + 1 ; %do i = 1 %to &groupnum ; if _i_ = &i then output __g&i ; %end ; output __g0 ; /* 全体データ */ run ; /*----------------------------------------------- 以下は proc IML （行列演算言語）による計算 -----------------------------------------------*/ proc iml ; /*-------- 汎用的なモジュールを準備 -----------*/ start csscp( x ) ; /* 偏差平方和積和行列を返す */ n = nrow( x ) ; sum = x[ + , ] ; xpx = t(x) * x - t(sum) * sum / n ; return( xpx ) ; finish csscp ; start cov( x ) ; /* 不偏分散共分散行列を返す */ n = nrow( x ) ; sum = x[ + , ] ; xpx = t(x) * x - t(sum) * sum / n ; df = n - 1 ; cov = ( 1/df ) * xpx ; return( cov ) ; finish cov ; start corr( x ) ; /* 相関行列を返す */ n = nrow( x ) ; sum = x[ + , ] ; xpx = t(x) * x - t(sum) * sum / n ; df = n - 1 ; cov = ( 1/df ) * xpx ; s = diag( 1 / sqrt( vecdiag( xpx ) ) ) ; corr = s * xpx * s ; return( corr ) ; finish corr ; start std( x ) ; /* 標準化データ行列を返す */ n = nrow( x ) ; df = n - 1 ; mean = x[ + , ] / n ; x = x - repeat( mean, n, 1 ) ; ss = x[ ## , ] ; std = sqrt( ss / df ) ; z = x * diag( 1 / std ) ; /* 標準得点の積和の平均は相関係数 */ corr = ( 1 / df ) * z` * z ; return( z ) ; finish std ; /*------ ＳＡＳデータセットを行列に変換 -------*/ %do i = 0 %to &groupnum ; /* 生データ行列:全体(g0)と群別(g1..) */ use __g&i ; read all var{ &var } into g&i ; close __g&i ; %end ; use _tmean_ ; read all var{ &var } into tmean ;/* 全体の平均：行ベクトル */ close _tmean_ ; use _mean_ ; read all var{ &var } into wmean ;/* 群別平均行列：群×変数 */ close _mean_ ; use _mean_ ; read all var{ _freq_ } into freq ;/* 各群個体数：列ベクトル*/ close _mean_ ; nvar = ncol( wmean ) ; /* 独立変数の数 ：スカラー */ ngrp = nrow( wmean ) ; /* 群の数 ：スカラー */ nobs = freq[ + , ] ; /* 全体の個体数 ：スカラー */ vname = { &var } ; /* 独立変数名の文字列：横ベクトル */ ndim = min( ngrp - 1 , nvar ) ; /* 次元数 */ print vname, nvar ngrp freq ndim , '全体平均',tmean, '群別平均',wmean ; /*------------------------------------------------------------- 偏差平方和積和行列からの計算 田中・垂水(1995)統計解析ハンドブック・多変量解析．共立出版 --------------------------------------------------------------*/ freqx = repeat( freq, 1, nvar ) ; tmeanz = wmean` - repeat(tmean`,1,ngrp); /* 群平均−全平均 */ Bsscp = (freqx` # tmeanz) * tmeanz`; /* 群間平方和積和行列 */ Wsscp = j( nvar , nvar, 0 ) ; /* プールした群内平方和積和行列 */ Tsscp = j( nvar , nvar, 0 ) ; /* 全体平方和積和行列 */ %do i = 1 %to &groupnum ; zg&i = g&i - repeat(wmean[&i,],freq[&i,],1); /* 群別平均偏差 */ wg&i = zg&i.` * zg&i ; /* 郡別の平方和積和行列 */ wsscp = wsscp + wg&i ; tg&i = g&i - repeat(tmean, freq[&i,], 1);/* 全平均との差 */ tg&i = tg&i.` * tg&i ; /* 郡平方和積和行列 */ tsscp = tsscp + tg&i ; %end ; WaddB = WSSCP + BSSCP ; title '田中・垂水(1995)統計解析ハンドブック・多変量解析．共立出版' ; print / Wsscp, Bsscp, Tsscp,, '平方和の加法定理(W+B)の確認', WaddB ; /* （Ｂ−λＷ）ａ＝０の一般固有値問題 */ call geneig( Lamda, a, Bsscp, Wsscp ) ; Lamda = Lamda[ 1:ndim ] ; a = a[ , 1:ndim ] ; /* 相関比: η2 = Sb / St の最大化問題 */ call geneig( eta2, e, bsscp, tsscp ) ; eta2 = eta2[ 1:ndim ] ; e = e[ , 1:ndim ] ; print / Lamda[format=14.10] a[format=14.10] , eta2 [format=14.10] e[format=14.10] ; /*--- 固有ベクトルの基準化 ---*/ a1 = a[ , 1 ] ; a1 = a1 / sqrt( (a1` * Wsscp * a1) / (nobs - ngrp) ) ; e1 = e[ , 1 ] ; e1 = e1 / sqrt( (e1` * Wsscp * e1) / (nobs - ngrp) ) ; print '分散比：固有ベクトルａの基準化' a1[ format = 14.10 ] ; print '相関比：固有ベクトルｅの基準化' e1[ format = 14.10 ] ; /*-------------------------------------------------------------- 分散共分散行列による計算 豊田(1996)非線形多変量解析．朝倉書店 --------------------------------------------------------------*/ Y = t( g0 ) ; /* 全ての群を込みにしたデータ行列：変数×個体 */ Yt = repeat( t( tmean ),1, nobs );/* 全平均の行列：変数×個体 */ twmean = t( wmean ) ; %do i = 1 %to &groupnum ; twmean&i = twmean[ ,&i ] ; freq&i = freq[ &i, ] ; wmean&i = repeat( twmean&i , 1, freq&i ) ; Yb = Yb || wmean&i ; /* 各群の平均値を個体数並べる */ %end ; St = (1/nobs ) * (Y - Yt) * (Y - Yt)` ;/* 全体分散共分散行列 */ Sw = (1/nobs ) * (Y - Yb) * (Y - Yb)` ;/* 群内分散共分散行列 */ Sb = (1/nobs ) * (Yb - Yt) * (Yb - Yt)` ;/* 群間分散共分散行列 */ BaddW = Sw + Sb ; call geneig( eta2, w, Sb, St ) ; call geneig( lamda, a, Sb, Sw ) ; Lamda = Lamda[ 1:ndim ] ; a = a[ , 1:ndim ] ; eta2 = eta2[ 1:ndim ] ; w = w[ , 1:ndim ] ; title '豊田(1996)非線形多変量解析．朝倉書店' ; print / '標本分散共分散行列', st, sw, sb , baddw / lamda[format=14.10] a[format=14.10] , eta2 [format=14.10] w[format=14.10] ; /*--- 固有ベクトルの基準化 ---*/ w1 = w[ , 1 ] ; w1 = w1 / sqrt( w1` * sw * w1 ) ; a1 = a[ , 1 ] ; a1 = a1 / sqrt( a1` * sw * a1 ) ; print '分散比：固有ベクトルａの基準化' a1[format=14.10] ; print '相関比：固有ベクトルｗの基準化' w1[format=14.10] ; /*-------------------------------------------------------------- Gnanadesikan(1977)Method for Statistical Data Analysis of Multivariate Observations. John Wiley & Sons. 丘本正・磯貝恭史訳(1979)統計的多変量データ解析．日科技連．p.75 --------------------------------------------------------------*/ title 'ニャナデシカン（丘本・磯貝訳(1979)）統計的多変量データ解析' ; /*--- ｐ×ｐのプールした群内（不偏）分散共分散行列 ---*/ wdf = nobs - ngrp ; wss = j( nvar, nvar, 0 ) ; %do i = 1 %to &groupnum ; wss&i = csscp( g&i ) ; wss = wss + wss&i ; %end ; wcov = ( 1/wdf ) * wss ; /*--- ｐ×ｐの群間（不偏）分散共分散行列 ---*/ bdf = ngrp - 1 ; bss = (freqx` # tmeanz) * tmeanz`;/* 群間平行和積和行列 */ bcov = ( 1/bdf ) * bss ; /*--- ｐ×ｐの全体（不偏）分散共分散行列 ---*/ tcov = cov( g0 ) ; call geneig( lamda, a, bcov, wcov ) ; Lamda = Lamda[ 1:ndim ] ; a = a[ , 1:ndim ] ; print / '不偏分散共分散行列', wcov, bcov, tcov ; print / lamda[format=18.10] a[format=18.10] ; t = half( wcov ) ; /* 群内のコレスキー分解 */ Bt = t(inv(t)) * bcov * inv(t) ; /* 変換後の群間分散行列 */ call eigen( Lamda, ev , bt ) ; /* 対称行列の固有値問題 */ Lamda = Lamda[ 1:ndim ] ; ev = ev[ , 1:ndim ] ; ev = inv( t ) * ev ; print t , bt , Lamda[format=18.10] ev[format=18.10] ; title '不偏分散による計算' ; bdf = ( nobs * ( ngrp - 1 ) ) / ngrp ; bcov = ( 1/bdf ) * bss ; %* bcov = tcov - wcov ; call geneig( lamda, a, bcov, wcov ) ; Lamda = Lamda[ 1:ndim ] ; a = a[ , 1:ndim ] ; print / '不偏分散共分散行列', wcov, bcov, tcov ; print / lamda[format=18.10] a[format=18.10] ; /*-------------------------------------------------------------- 松原望(1992)統計的決定．放送大学．p.119 --------------------------------------------------------------*/ title '松原望(1992)統計的決定' ; /*--- ｐ×ｐのプールした群内（不偏）分散共分散行列 ---*/ df = nobs - 1 ; wss = j( nvar, nvar, 0 ) ; %do i = 1 %to &groupnum ; wss&i = csscp( g&i ) ; wss = wss + wss&i ; %end ; wcov = ( 1/df ) * wss ; /*--- ｐ×ｐの群間（不偏）分散共分散行列 ---*/ bss = (freqx` # tmeanz) * tmeanz`;/* 群間平行和積和行列 */ bcov = ( 1/df ) * bss ; /*--- ｐ×ｐの全体（不偏）分散共分散行列 ---*/ tcov = cov( g0 ) ; call geneig( lamda, a, bcov, wcov ) ; Lamda = Lamda[ 1:ndim ] ; a = a[ , 1:ndim ] ; print / '松原(1992)', wcov, bcov, tcov ; print / lamda[format=18.10] a[format=18.10] ; title ; run ; quit ; %mend candisc ; options source ; %candisc( data = dat.irisx, /* Fisher の菖蒲データ */ class = spec, var = v1 v2 v3 v4 ) title ' PROC candisc による結果' ; proc candisc data = dat.irisx all out = canout ; class spec ; var v1 - v4 ; run ; title '正準判別変量の確認' ; proc candisc data = canout( keep = can1 can2 spec ) all ; class spec ; var can: ; run ; title 'PROC glm による多変量分散分析の結果' ; proc glm data = dat.irisx ; class spec ; model v1 v2 v3 v4 = spec ; manova h = spec / printe printh ; run ; quit ; title ;