home *** CD-ROM | disk | FTP | other *** search
/ ftp.barnyard.co.uk / 2015.02.ftp.barnyard.co.uk.tar / ftp.barnyard.co.uk / cpm / walnut-creek-CDROM / CPM / TURBOPAS / MAPSTATF.LBR / FACTOR.PZS / FACTOR.ÐAS
Text File  |  2000-06-30  |  7KB  |  219 lines
  1. (* Multivariate Analysis Package - Principle Component Factoring Module
  2.    Copyright 1985 Douglas L. Anderton.  This program may be freely
  3.    circulated so long as it is not sold for profit and any charge does
  4.    not exceed costs of reproduction *)
  5.  
  6. Program Factor(Input,Output);
  7. Const
  8.   N=30;
  9. Type
  10.   SUBS=1..N;
  11.   RVEC = Array [SUBS] Of Real;
  12.   NBYN = Array [SUBS] Of RVEC;
  13.   IVEC = Array [SUBS] Of Integer;
  14.   S8 = Array [SUBS] Of String[8];
  15. Var
  16.   dfile, ofile : Text;
  17.   sel : IVEC;
  18.   e, u, v, w, x, y: RVEC;
  19.   varn : S8;
  20.   cor, a : NBYN;
  21.   nf, nv, dv, i, j, k, rn, ot: Integer;
  22.   tol, temp, tmp2, nc : Real;
  23.   error : boolean;
  24.  
  25. Procedure openfiles(Var dfile, ofile:Text; Var nv, dv, ot:Integer);
  26. Var
  27.   dfl, ofl:String[12];
  28. Begin
  29.   ClrScr; Writeln(' *** FACTOR: PRINCIPLE COMPONENT FACTORING ***');
  30.   Writeln; Write('Name of the CORREL file? ');
  31.   Readln(dfl); Assign(dfile,dfl); Reset(dfile);
  32.   Write('Name of the output file? ');
  33.   Readln(ofl); Assign(ofile,ofl); Rewrite(ofile);
  34.   ot:= 0;
  35.   If (ofl='CON:') Or (ofl='con:') Then ot:=1;
  36.   If (ofl='LST:') Or (ofl='lst:') Then ot:=2;
  37.   If (ot=2) Then
  38.     Begin
  39.     Writeln(ofile,'Multivariate Analysis Package (1.6) - ',
  40.       'Program: FACTOR, Datafile: ',dfl); Writeln(ofile);
  41.     End;
  42.   Write('How many variables in CORREL matrix? '); Readln(nv);
  43.   Write('Number of variables to use in FACTOR? '); Readln(dv);
  44.   End; (* Of openfiles *)
  45.  
  46. Procedure selectvar2(Var sel:IVEC; dv:Integer);
  47. Var
  48.   i:Integer;
  49. Begin
  50.   Writeln;
  51.   For i:=1 To dv Do Begin
  52.     Write('Column number for variable ',i,'? '); Readln(sel[i]); End;
  53.   End; (* Of selectvar *)
  54.  
  55. Procedure wait(ot:Integer);
  56. Begin
  57.   If ot=1 Then Begin
  58.     Write('- Press any key to continue -'); While Not KeyPressed Do; ClrScr;
  59.     End;
  60.   End; (* of wait *)
  61.  
  62. {$I GETCOR.LIB}
  63. {$I EIGEN.LIB}
  64. {$I MATINV.LIB}
  65.  
  66. Procedure rotate(Var a:NBYN; nf, nv, meth, ot:Integer);
  67. Var
  68.   c:RVEC;
  69.   i,j,k,f,nr:Integer;
  70.   ep,em,a1,b1,c1,d1,u,v,qn,qd,cs,sn,cs1,sn1,sp,cp:Real;
  71. Begin
  72.   Writeln(ofile);
  73.   If(meth=0) Then Writeln(ofile,'Varimax Rotation of ',nf:2,' Factors')
  74.   Else Writeln(ofile,'Quartimax Rotation of ',nf:2,' Factors');
  75.   ep:=0.00116;
  76.   For j:=1 To nv Do
  77.     Begin
  78.     c[j]:=0.0;
  79.     For k:=1 To nf Do c[j]:=c[j]+Sqr(a[j,k]);
  80.     c[j]:=Sqrt(c[j]);
  81.     For k:=1 To nf Do a[j,k]:=a[j,k]/c[j];
  82.   End;
  83.   While((nr-((nf*(nf-1))/2))<>0) Do
  84.     Begin
  85.     nr:=0;
  86.     For i:=1 To (nf-1) Do
  87.       Begin
  88.       For j:=(i+1) To nf Do
  89.         Begin
  90.         a1:=0.0; b1:=0.0; c1:=0.0; d1:=0.0;
  91.         For k:=1 to nv Do Begin
  92.           u:=Sqr(a[k,i])-Sqr(a[k,j]); v:=a[k,i]*a[k,j]*2.0;
  93.           a1:=a1+u; b1:=b1+v; c1:=c1+Sqr(u)-Sqr(v); d1:=d1+(u*v*2.0); End;
  94.         If(meth=0) Then
  95.           Begin qn:=d1-((2*a1*b1)/nv); qd:=c1-((Sqr(a1)-Sqr(b1))/nv); End;
  96.         If(meth=1) Then Begin qn:=d1; qd:=c1; End;
  97.         f:=0;
  98.         If((Abs(qn)+Abs(qd))=0) Then Begin nr:=nr+1; f:=1; End;
  99.         If((Abs(qn)-Abs(qd))=0) and (f=0) Then
  100.           Begin cs:=0.70710678; sn:=cs; f:=2; End;
  101.         If((Abs(qn)-Abs(qd))>0) and (f=0) Then
  102.           Begin
  103.           em:=Abs(qd/qn);
  104.           If(em<ep) Then Begin cs:=0; sn:=1; End
  105.           Else Begin sn:=1/Sqrt(1+Sqr(em)); cs:=sn*em; End;
  106.           f:=2;
  107.           End;
  108.         If(f=0) Then
  109.           Begin
  110.           em:=Abs(qn/qd);
  111.           If(ep<=em) Then Begin
  112.             cs:=Cos(ArcTan(em)); sn:=Sin(ArcTan(em)); f:=2; End
  113.           Else If(qd>=0) Then Begin nr:=nr+1; f:=1; End
  114.             Else Begin sp:=0.70710678; cp:=sp; f:=3; End
  115.           End;
  116.         If(f=2) Then
  117.           Begin
  118.           em:=Sqrt((1+cs)*0.5); cs1:=Sqrt((1+em)*0.5); sn1:=sn/(4*cs1*em);
  119.           If(qd<0) Then
  120.             Begin cp:=0.70710678*(cs1+sn1); sp:=0.70710678*(cs1-sn1); End
  121.           Else Begin cp:=cs1; sp:=sn1; End;
  122.           If(qn<0) Then sp:=-sp;
  123.           End;
  124.         If((f=2) or (f=3)) Then
  125.           For k:=1 To nv Do
  126.             Begin em:=(a[k,i]*cp)+(a[k,j]*sp);
  127.             a[k,j]:=(a[k,j]*cp)-(a[k,i]*sp); a[k,i]:=em; End;
  128.         End;
  129.       End;
  130.     End;
  131.   For i:=1 To nv Do For j:=1 To nf Do a[i,j]:=a[i,j]*c[i];
  132.   For j:=1 To nf Do
  133.     Begin c[j]:=0.0; For k:=1 To nv Do c[j]:=c[j]+Sqr(a[k,j]); End;
  134.   ClrScr;
  135.   For i:=1 To nf Do Begin
  136.     Writeln(ofile); Writeln(ofile,'Rotated Factor Pattern For Factor: ',i:2);
  137.     For j:=1 To dv Do
  138.       Begin Writeln(ofile,'   ',varn[j]:8,': ',a[j,i]:10:8); End;
  139.     Writeln(ofile,'Factor Accounts For ',(c[i]*100/nv):8:5,'% of Variance');
  140.     wait(ot);
  141.     End;
  142.   End; (* of rotate *)
  143.  
  144. Begin (* main *)
  145.   openfiles(dfile, ofile, nv, dv, ot);
  146.   Write('Number of Factors to Extract (<=number of variables)? ');
  147.   Readln(nf); If (nf=dv) Then rn:=0 Else rn:=1;
  148.   tol:=1E-8;
  149.   Write('Tolerance for eigenvalue convergence <return for 1E-8>? ');
  150.   Readln(tol);
  151.   selectvar2(sel, dv);
  152.   (* read correlation matrix matrix *)
  153.   getcor(cor, x, y, varn, nc, sel, nv, dv, dfile);
  154.   (* full rank analysis *)
  155.   If(rn=0) Then
  156.     Begin
  157.     a:=cor;
  158.     (* invert *)
  159.     Writeln('full rank analysis inverting matrix...');
  160.     matinv(a,dv,temp,error);
  161.     ClrScr; Writeln(ofile,'Inverse:  (determinant=',temp,')');
  162.     If(error) Then Writeln(ofile,'Error reported in solving for Inverse');
  163.     temp:=-((nc-1.0)-((1.0/6.0)*(2.0*nc+5.0)))*Ln(temp);
  164.     Write(ofile,'Chi-Square for Sphericity = ',temp:10:3);
  165.     Writeln(ofile,'  degrees of freedom = ',((Sqr(dv)-dv)/2):6:0);
  166.     wait(ot);
  167.     For j:=1 To dv Do w[j]:=1.0-(1.0/a[j,j]);
  168.     End;
  169.   a:=cor;
  170.   eigen(cor, a, e, tol, dv, nf, ot);
  171.   wait(ot);
  172.   temp:=0.0; tol:=dv;
  173.   For j:=1 To nf Do Begin
  174.     u[j]:=(e[j]/tol)*100.0; temp:=temp+u[j]; v[j]:=temp;
  175.     tmp2:=Sqrt(e[j]);
  176.     For k:=1 To dv Do
  177.       Begin cor[k,j]:=a[k,j]*tmp2; a[k,j]:=a[k,j]/tmp2; End;
  178.     End;
  179.   ClrScr; Writeln(ofile);
  180.   Writeln(ofile,'Factor        Eigen          %Variance     %Cumulative');
  181.   Writeln(ofile);
  182.   For i:=1 To nf Do Writeln(ofile,'   ',i:2,'         ',e[i]:8:5,
  183.    '       ',u[i]:7:3,'       ',v[i]:7:3);
  184.   Writeln(ofile);
  185.   wait(ot);
  186.   For i:=1 To nf Do Begin
  187.     Writeln(ofile); Writeln(ofile,'Factor Pattern For Factor: ',i:2);
  188.     For j:=1 To dv Do
  189.       Begin Writeln(ofile,'   ',varn[j]:8,': ',cor[j,i]:10:8); End;
  190.     wait(ot);
  191.     End;
  192.   Writeln(ofile); Write(ofile,'Variable      Communality');
  193.   If(rn=0) Then Writeln(ofile,'     R-Square') Else Writeln(ofile);
  194.    For j:=1 To dv Do
  195.     Begin temp:=0.0;
  196.     For k:=1 To nf Do temp:=temp+Sqr(cor[j,k]);
  197. {DOUG    If temp>1.0 then temp:=1.0; }
  198.     Write(ofile,'   ',varn[j]:8,'    ',temp:10:7);
  199.     If(rn=0) Then Writeln(ofile,'     ',w[j]:10:7) Else Writeln(ofile);
  200.     End;
  201.   wait(ot);
  202.   For i:=1 To nf Do Begin
  203.     Writeln(ofile); Writeln(ofile,'Score Coefficients For Factor: ',i:2);
  204.      For j:=1 To dv Do
  205.       Begin Writeln(ofile,'   ',varn[j]:8,': ',a[j,i]:10:8); End;
  206.     wait(ot);
  207.     End;
  208.   If(nf>1) Then
  209.     Begin
  210.     Write('What type of Factor Rotation (0=none,1=varimax,2=quartimax)?');
  211.     Readln(i); i:=i-1; If(i>=0) Then rotate(cor, nf, dv, i, ot);
  212.     End;
  213.    Close(ofile); Close(dfile);
  214.    Assign(dfile,'MAPSTAT.COM'); Execute(dfile);
  215.   End.
  216.  
  217. :2);
  218.     For j:=1 To dv Do
  219.       Begin Writeln(ofile,'   ',va