/*                Alpha Factor Analysis                      */
/*  Ref: Kaiser et al., 1965 Psychometrika, pp. 12-13        */
/*  r correlation matrix (n.s.) already set up               */
/*  p number of variables                                    */
/*  q number of factors                                      */
/*  h communalities                                          */
/*  m eigenvalues                                            */
/*  e eigenvectors                                           */
/*  f factor pattern                                         */
/*  (IQ,H2,HI,G,MM) temporary use. freed up                  */
/*                                                           */

*ods trace output;
proc iml;
start alpha;
   p = ncol(r);
   q = 0;
   h = 0;                                      /* initialize */
   h2 = i(p)-diag(1/vecdiag(inv(r)));                /* smcs */
   do while(max(abs(h-h2))>.001); /* iterate until converges */
      h = h2;                      
      hi = diag(sqrt(1/vecdiag(h)));    
      g = hi*(r-i(p))*hi+i(p);     
      call eigen(m,e,g);         /* get eigenvalues and vecs */
      if q=0 then do;
         q = sum(m>1);                  /* number of factors */
         iq = 1:q;             
      end;                                   /* index vector */
      mm = diag(sqrt(m[iq,]));           /* collapse eigvals */
      e = e[,iq] ;                       /* collapse eigvecs */
      h2 = h*diag((e*mm) [,##]);        /* new communalities */
   end;
   hi = sqrt(h);
   h = vecdiag(h2);
   f = hi*e*mm;                         /* resulting pattern */
   free iq h2 hi g mm;                   /* free temporaries */
finish;

  /* Correlation Matrix from Harmon, Modern Factor Analysis, */
  /* Second 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};
nm = {Var1 Var2 Var3 Var4 Var5 Var6 Var7 Var8};
run alpha;
print ,"EIGENVALUES" , m;
print ,"COMMUNALITIES" , h[rowname=nm];
print ,"FACTOR PATTERN", f[rowname=nm];