探索的因子分析（EFA）の推定法の１つである
重みなし最小２乗法（ULS）のSAS/IMLによる
プログラム（ Minres 法）．
 Harman(1976) Modern Factor Analysis.
に基づいて書き始め，まだデバグ中です．
時間ができたら再開します．
（鈴木督久, Sep. 1999）


/*----------------------------------------------*/
/*    探索的因子分析の最小２乗法（ミンレス法）  */
/*    Unweighted Least Squares Method[ ULS ]    */
/*----------------------------------------------*/
/*                                              */
/*  R    : 相関行列                             */
/*  p    : 変数の個数                           */
/*  q    : 因子の個数                           */
/*  h    : 共通性                               */
/*  m    : 固有値                               */
/*  e    : 固有ベクトル                         */
/*  F    : 因子負荷行列                         */
/*  conv : 収束基準                             */
/*  iter : 反復回数                             */
/*                                              */
/*  Author: Tokuhisa SUZUKI                     */
/*----------------------------------------------*/
proc iml ;

start uls ;

/* 初期設定 */
	q       = 0         ; /* 明示指定も可 */
	conv    = 0.001     ; /* 明示指定も可 */
	maxiter = 20        ; /* 明示指定も可 */

	call eigen( val, vec, R )        ;
	if q = 0 then q = sum( val > 1 ) ;
	iq   = 1:q ;

/* 事前共通性と初期因子負荷 */
	p   = ncol( r ) ;
	smc = i(p) - diag(1 / vecdiag( inv(r) )); /* 芝(1979), p62 */
	one = i(p) ;
	h2  = smc ;                    /* 事前共通性 h2 の選択指定 */
	G   = ( R - diag(R) ) + h2  ;  /* 事前推定値の代入した行列 */
	call eigen(m, e, G ) ;
	eq  = e[ ,iq]        ;         /* 最初のｑ個の固有ベクトル */
	md  = diag( sqrt( m[iq, ] )) ; /* 最初のｑ個の固有値対角化 */
	F2  = eq * md ;                /* 初期因子負荷行列         */

	eigenval = t( val )          ; /* Ｒの固有値ベクトルを横に */
	preeigen = t( m   )          ; /* 事前固有値ベクトルを横に */
	priors   = t( vecdiag( h2 )) ; /* 事前共通性ベクトルを横に */

/* 反復推定 */
	F = 0 ;
	do iter = 1 to maxiter while( max(abs(F - F2)) > conv ) ;
		F = F2 ;
		free A ;
		do i = 1 to p ;
			run lineqn( F, G, q, p, i ) ;
			A = A // _a_ ;
		end ;
		F2 = A      ; /* 新しい因子負荷の算出 */
		G  = F2 * F2` ;

		communal = communal // t( vecdiag( F2 * F2` ) ) ;
		change   = change   // max( abs( F - F2 ) )     ;
		iters    = iters    // iter ;
	end ;

	F  = F2       ;           /* 因子負荷の最終推定値 */
	Rf = F2 * F2` ;           /* 再生した相関係数行列 */
	Rr = R - Rf   ;           /* 最終的な残差相関行列  */
	h  = t( vecdiag( Rf ) ) ; /* 最終的な共通性推定値 */
finish ;


start lineqn( A, R, q, p, j ) global ( _a_ ) ;
/*-------------------------------------*/
/*  A   前の因子負荷行列               */
/*  R   （再生）相関行列               */
/*  q   因子の個数                     */
/*  p   観測変数の数                   */
/*  j   観測変数の番号                 */
/*  a_  次の因子負荷                   */
/*-------------------------------------*/
	_F = A ;
	_F[ j, 1:q ] = 0 ; /* j番変数を0に */
	rj = R[ 1:p, j ] ; /* j列目をとる  */
	rj[ j, 1 ] = 0   ; /* 共通性を0に  */
	_a_ = t( inv( _F` * _F ) * _F` * rj ) ;
finish ;


/* correlation matrix from harman, modern factor analysis,
   2nd edition, page 124, "eight physical variables"    */

r = { 1.000    .846    .805    .859    .473    .398    .301    .382 ,
       .846   1.000    .881    .826    .376    .326    .277    .415 ,
       .805    .881   1.000    .801    .380    .319    .237    .345 ,
       .859    .826    .801   1.000    .436    .329    .327    .365 ,
       .473    .376    .380    .436   1.000    .762    .730    .629 ,
       .398    .326    .319    .329    .762   1.000    .583    .577 ,
       .301    .277    .237    .327    .730    .583   1.000    .539 ,
       .382    .415    .345    .365    .629    .577    .539   1.000 
	} ;

ev = { "1" "2" "3" "4" "5" "6" "7" "8" } ;
nm = { V1 V2 V3 V4 V5 V6 V7 V8 } ;

run uls ;

title "SAS/IML による最小２乗法（ミンレス法）" ;
print / "相関行列の固有値"   , eigenval[ colname = ev format = 9.5 ] ,
        "共通性の事前推定値" , priors  [ colname = nm format = 9.6 ] ,
        "事前共通性の固有値" , preeigen[ colname = ev format = 9.4 ] ;
print / "反復過程" , iters change[ format = 9.6 ] communal[ format = 9.5 ] ;
print / "再生相関行列" , Rf[ format = 9.4               ] ,
        "因子負荷行列" , F [ format = 9.5 rowname = nm  ] ,
        "残差相関行列" , Rr[ format = 9.5 rowname = nm  ] ,
        "共通性推定値" , h [ format = 9.6               ] ;
quit  ;
title ;



/* Harman (1976) Modern factor analysis. 3rd Ed. p.22" */
data x( type  = corr
        label = "Eight Physiacl Variables for 305 girls" ) ;
    input _type_$ _name_$ v1 - v8 ;
    label
        v1 = 'Height'
        v2 = 'Arm span'
        v3 = 'Length of forearm'
        v4 = 'Length of lower leg'
        v5 = 'Weight'
        v6 = 'Bitrochanteric diameter'
        v7 = 'Chest girth'
        v8 = 'Chest width'
    ;
cards;
N    .    305     305     305     305     305     305     305     305     
MEAN .      0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000 
STD  .      1.000   1.000   1.000   1.000   1.000   1.000   1.000   1.000 
CORR v1     1.000    .846    .805    .859    .473    .398    .301    .382 
CORR v2      .846   1.000    .881    .826    .376    .326    .277    .415 
CORR v3      .805    .881   1.000    .801    .380    .319    .237    .345 
CORR v4      .859    .826    .801   1.000    .436    .329    .327    .365 
CORR v5      .473    .376    .380    .436   1.000    .762    .730    .629 
CORR v6      .398    .326    .319    .329    .762   1.000    .583    .577 
CORR v7      .301    .277    .237    .327    .730    .583   1.000    .539 
CORR v8      .382    .415    .345    .365    .629    .577    .539   1.000 
;
run ;

title "PROC FACTOR による ULS" ;
* proc factor data = x m = uls ev ;
run ;
title ;