/*----------------------------------------------*/
/*                  反復主因子法                */
/*    iterated  Principal Factor Analysis       */
/*                    [iPFA]                    */
/*----------------------------------------------*/
/*                                              */
/*     r  相関行列                              */
/*     p  変数の個数                            */
/*     q  因子の個数                            */
/*     h  共通性                                */
/*     m  固有値                                */
/*     e  固有ベクトル                          */
/*     f  因子負荷行列                          */
/*     m1 事前の固有値（１回目の固有値問題）    */
/*     e1 事前の固有ベクトル                    */
/*     conv 収束基準                            */
/*     iter 反復回数                            */
/*     (iq,h2,hi,g,mm) 一時的使用               */
/*                                              */
/*                                              */
/*  code by T.SUZUKI                            */
/*----------------------------------------------*/

title "SAS/IML による反復主因子法" ;

proc iml ;
start pfa ;

    /* 初期設定 */
    p       = ncol( r ) ;
    q       = 0         ;
    h       = 0         ;
    conv    = 0.001     ;

    smc = i(p) - diag(1 / vecdiag( inv(r) )); /*芝(1979)p62 */
    one = i(p)                              ;

    h2 = one ;                  /* 共通性の事前推定値の指定 */

    g  = ( r - diag(r) ) + h2 ; /* 共通性の事前推定値の代入 */

    if q = 0 then   /* 因子数が無指定の場合のデフォルト設定 */
        do ;
            call eigen( val, vec, r ) ;
            q  = sum( val > 1 ) ; /* 固有値＞１の個数を採用 */
            iq = 1:q            ;
        end ;

    eigenval = t( val           ) ; /* 相関行列の固有値を横ベクトルに */
    priors   = t( vecdiag( h2 ) ) ; /* 事前共通性推定値を横ベクトルに */
    call eigen(m1, e1, g )        ; /* 事前推定値の固有値問題を解く   */
    preeigen = t( m1            ) ; /* 事前の固有値を横ベクトルに     */

    /* 収束するまで共通性を反復推定する */
    iter = 0 ;
    do while( max( abs( h - h2 ) ) > conv ) ;
        h  = h2 ;
        g  = ( g - diag(g) ) + h ;     /* 共通性の代入       */
        call eigen(m, e, g )  ;        /* 固有値問題         */
        mm = diag( sqrt( m[iq, ] ) ) ; /* 固有値の対角化     */
        e  = e[ ,iq] ;                 /* ｑ個の固有ベクトル */
        f  = e * mm ;                  /* 因子負荷の推定値   */
        h2 = diag( f * f` ) ;          /* 新しい共通性の算出 */
        iter = iter + 1 ;

        vec_h  = vec_h  // t( vecdiag( h2 ) )   ;
        change = change // max( abs( h - h2 ) ) ;
        iters  = iters  // iter ;
    end;

    m  = t( m )             ; /* 最終的な縮退行列固有値 */
    h  = t( vecdiag( h2 ) ) ; /* 最終的な共通性の推定値 */ 
    f  = e * mm             ; /* 最終的な因子負荷推定値 */
    free iq h2 hi g mm ;      /* free temporaries       */
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 pfa ;


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 ] vec_h[ format = 9.5 ]  ;
print / "縮退した相関行列の固有値" , m[ format = 9.4 ] ,
        "因子負荷行列"             , f[ rowname = nm format = 9.5 ] ,
        "最終的な共通性推定値"     , 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 による反復主因子法" ;
proc factor data = x m = prinit ;
run ;
title ;