home *** CD-ROM | disk | FTP | other *** search
/ CP/M / CPM_CDROM.iso / lambda / soundpot / f / map12.lbr / FACTOR.PZS / FACTOR.PAS
Encoding:
Pascal/Delphi Source File  |  1993-10-26  |  10.5 KB  |  379 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. Type
  8.   R35 = Array [1..35] Of Real;
  9.   M35 = Array [1..35,1..35] Of Real;
  10.   I35 = Array [1..35] Of Integer;
  11.   S8 = Array [1..35] Of String[8];
  12. Var
  13.   dfile, ofile : Text;
  14.   sel : I35;
  15.   e, u, v, w, x, y: R35;
  16.   varn : S8;
  17.   cor, a, vc : M35;
  18.   i, j, k, nf, rn, nv, dv, ot, nc: Integer;
  19.   temp : Real;
  20.  
  21. Procedure openfiles(Var dfile, ofile:Text; Var nv, dv, ot:Integer);
  22. Var
  23.   dfl, ofl:String[12];
  24. Begin
  25.   ClrScr; Writeln(' *** FACTOR: PRINCIPLE COMPONENT FACTORING ***');
  26.   Writeln; Write('Name of the CORREL file? ');
  27.   Readln(dfl); Assign(dfile,dfl); Reset(dfile);
  28.   Write('Name of the output file? ');
  29.   Readln(ofl); Assign(ofile,ofl); Rewrite(ofile);
  30.   ot:= 0;
  31.   If (ofl='CON:') Or (ofl='con:') Then ot:=1;
  32.   If (ofl='LST:') Or (ofl='lst:') Then ot:=2;
  33.   If (ot=2) Then
  34.     Begin
  35.     Writeln(ofile,'Multivariate Analysis Package (1.2) - ',
  36.       'Program: FACTOR, Datafile: ',dfl); Writeln(ofile);
  37.     End;
  38.   Write('How many variables in CORREL matrix? ');
  39.   Readln(nv);
  40.   Write('Number of variables to use in FACTOR? ');
  41.   Readln(dv);
  42.   End; (* Of openfiles *)
  43.  
  44. Procedure selectvar(Var sel:I35; Var varn:S8; dv:Integer);
  45. Var
  46.   i:Integer;
  47. Begin
  48.   Writeln;
  49.   For i:=1 To dv Do
  50.     Begin
  51.     Write('Column number for variable ',i,'? ');
  52.     Read(sel[i]); Write('  Name? '); Readln(varn[i]);
  53.     End;
  54.   End; (* Of selectvar *)
  55.  
  56. Procedure wait(ot:Integer);
  57. Begin
  58.   If ot=1 Then
  59.     Begin
  60.     Write('- Press any key to continue -'); While Not KeyPressed Do;
  61.     ClrScr;
  62.     End;
  63.   End; (* of wait *)
  64.  
  65. Procedure getcor(Var cor:M35; Var mean, stdev:R35; Var nc:Integer;
  66.                sel:I35; nv, dv:Integer; Var dfile:Text);
  67. Var
  68.   i,j,k,l,z:Integer;
  69.   x,y:Real;
  70.   dum15:String[15];
  71.   dum21:String[21];
  72.   dum8:String[8];
  73.   dum36:String[36];
  74. Begin
  75.   Readln(dfile); Readln(dfile);
  76.   For i := 1 To nv Do
  77.     Begin
  78.     Readln(dfile,dum15,x,dum21,y,dum8,z);
  79.     If (i=1) Then nc:=z;
  80.     If (nc>z) Then nc:=z;
  81.     For j := 1 To dv Do
  82.       Begin
  83.       If (sel[j]=i) Then
  84.         Begin
  85.         mean[j] := x; stdev[j] := y;
  86.         End;
  87.       End;
  88.     End;
  89.     For i :=1 to 3 Do Readln(dfile);
  90.     For i := 1 To nv-1 Do
  91.       For k := i+1 To nv Do
  92.         Begin
  93.         Readln(dfile,dum36,x,dum15,y);
  94.         For j := 1 To dv Do
  95.           Begin
  96.           cor[j,j] := 1.0;
  97.           If (sel[j]=i) Then
  98.             For l := 1 To dv Do
  99.               If (sel[l]=k) Then
  100.               Begin
  101.               cor[j,l] := y; cor[l,j] := y;
  102.               End;
  103.           End;
  104.         End;
  105.   End; (* Of getcor *)
  106.  
  107. Procedure eigen(Var a, v:M35; Var e:R35; Var dv, nf, ot:Integer);
  108. (* eigenvalues and vectors by iteration and exhaustion - a is input
  109.   matrix of dimension dv, nf the number of vectors to extract, e is
  110.   returned eigenvalues and v is returned eigen vectors *)
  111. Var
  112.   x,y:R35;
  113.   j,k,l,it:Integer;
  114.   t,s1,s2:Real;
  115. Begin
  116.   l:=1; t:=0.00001;
  117.   While (((dv-l)>0) and (l<=nf)) Do
  118.     Begin
  119.     it:=0; s1:=1.0;
  120.     For j:=1 To dv Do y[j]:=1.0;
  121.     While ((s1-t)>0.0) Do
  122.       Begin
  123.       it:=it+1;
  124.       For j:=1 To dv Do
  125.         Begin
  126.         x[j]:=0.0;
  127.         For k:=1 To dv Do x[j]:=x[j]+(a[j,k]*y[k]);
  128.         End;
  129.       e[l]:=x[1]; s1:=0.0;
  130.       For j:=1 To dv Do
  131.         Begin
  132.         v[j,l]:=x[j]/x[1]; s1:=s1+(Abs(y[j]-v[j,l])); y[j]:=v[j,l];
  133.         End;
  134.       If ((it=20) and ((s2-s1)<=0)) Then Begin
  135.         Writeln;
  136.         Writeln('*** Roots Not Converging in Eigen - Error Abort ***');
  137.         Bdos(0); End;
  138.       s2:=s1;
  139.       End;
  140.     s1:=0.0;
  141.     For j:=1 To dv Do s1:=s1+(v[j,l]*v[j,l]); s1:=Sqrt(s1);
  142.     For j:=1 To dv Do v[j,l]:=v[j,l]/s1;
  143.     For j:=1 To dv Do
  144.       For k:=1 To dv Do a[j,k]:=a[j,k]-(v[j,l]*v[k,l]*e[l]);
  145.     If(e[l]<1.0) Then
  146.       Begin
  147.       Writeln;
  148.       nf:=l;
  149.       Writeln('** Number of Eigenvalues Extracted Limited to ',nf:2,' **');
  150.       Writeln('** Variance Explained by Factor too small For Accuracy **');
  151.       wait(ot);
  152.       End;
  153.     l:=l+1;
  154.     End;
  155. End; (* of eigen *)
  156.  
  157. Procedure minv(Var a:M35; Var d:Real; m:Integer);
  158. (* matrix inversion in pascal: gauss - reduction method
  159.   inverse returned in a and determinant in d, m is dimension *)
  160. Var
  161.   j,k,l:Integer;
  162.   pvt,t:Real;
  163. Begin
  164.   d := 1.0;
  165.   For j:=1 To m Do
  166.     Begin
  167.     pvt := a[j,j]; d := d*pvt; a[j,j] := 1.0;
  168.     If (pvt=0) Then 
  169.       Begin
  170.       Writeln; Writeln('  *** Multicollinear Matrix ***'); Writeln;
  171.       End;
  172.     For k:=1 To m Do a[j,k] := a[j,k]/pvt;
  173.     For k:=1 To m Do
  174.       If ((k-j) <> 0) Then
  175.         Begin
  176.         t := a[k,j]; a[k,j] := 0.0;
  177.         For l:=1 To m Do a[k,l] := a[k,l] - (a[j,l]*t);
  178.         End;
  179.     End;
  180.   End; (* of minv *)
  181.  
  182. Procedure rotate(Var a:M35; nf, nv, meth, ot:Integer);
  183. Var
  184.   c:R35;
  185.   i,j,k,f,nr:Integer;
  186.   ep,em,a1,b1,c1,d1,u,v,qn,qd,cs,sn,cs1,sn1,sp,cp:Real;
  187. Begin
  188.   Writeln(ofile);
  189.   If(meth=0) Then Writeln(ofile,'Varimax Rotation of ',nf:2,' Factors')
  190.   Else Writeln(ofile,'Quartimax Rotation of ',nf:2,' Factors');
  191.   ep:=0.00116;
  192.   For j:=1 To nv Do
  193.     Begin
  194.     c[j]:=0.0;
  195.     For k:=1 To nf Do c[j]:=c[j]+Sqr(a[j,k]);
  196.     c[j]:=Sqrt(c[j]);
  197.     For k:=1 To nf Do a[j,k]:=a[j,k]/c[j];
  198.   End;
  199.   While((nr-((nf*(nf-1))/2))<>0) Do
  200.     Begin
  201.     nr:=0;
  202.     For i:=1 To (nf-1) Do
  203.       Begin
  204.       For j:=(I+1) To nf Do
  205.         Begin
  206.         a1:=0; b1:=0; c1:=0; d1:=0;
  207.         For k:=1 to nv Do
  208.           Begin
  209.           u:=Sqr(a[k,i])-Sqr(a[k,j]); v:=a[k,i]*a[k,j]*2;
  210.           a1:=a1+u; b1:=b1+v; c1:=c1+Sqr(u)-Sqr(v); d1:=d1+(u*v*2);
  211.           End;
  212.         If(meth=1) Then
  213.           Begin 
  214.           qn:=d1-((2*a1*b1)/nv); qd:=c1-((Sqr(a1)-Sqr(b1))/nv);
  215.           End;
  216.         If(meth=0) Then
  217.           Begin
  218.           qn:=d1; qd:=c1;
  219.           End;
  220.         f:=0;
  221.         If((Abs(qn)+Abs(qd))<1) Then
  222.           Begin
  223.           nr:=nr+1; f:=1;
  224.           End;
  225.         If((Abs(qn)-Abs(qd))=0) Then
  226.           Begin
  227.           cs:=0.70710678; sn:=cs; f:=2;
  228.           End;
  229.         If((Abs(qn)-Abs(qd))>1) Then
  230.           Begin
  231.           em:=Abs(qd/qn);
  232.           If((em-ep)>=0) Then
  233.             Begin
  234.             cs:=0; sn:=1;
  235.             End
  236.           Else
  237.             Begin
  238.             sn:=1/Sqrt(1+Sqr(em)); cs:=sn*em;
  239.             End;
  240.           f:=2;
  241.           End;
  242.         If(f=0) Then
  243.           Begin
  244.           em:=Abs(qn/qd);
  245.           If((em-ep)>=0) Then 
  246.             Begin
  247.             cs:=Cos(ArcTan(em)); sn:=Sin(ArcTan(em)); f:=2;
  248.             End
  249.           Else
  250.             If(qd>=0) Then
  251.               Begin
  252.               nr:=nr+1; f:=1;
  253.               End
  254.             Else
  255.               Begin
  256.               sp:=0.70710678; cp:=sp; f:=3;
  257.               End
  258.           End;
  259.         If(f=2) Then
  260.           Begin
  261.           em:=Sqrt((1+cs)*0.5); cs1:=Sqrt((1+em)*0.5); sn1:=sn/(4*cs1*em);
  262.           If(qd<0) Then
  263.             Begin
  264.             cp:=0.70710678*(cs1+sn1); sp:=0.70710678*(cs1-sn1);
  265.             End
  266.           Else
  267.             Begin
  268.             cp:=cs1; sp:=sn1;
  269.             End;
  270.           If(qn<0) Then sp:=-sp;
  271.           End;
  272.         If((f=2) or (f=3)) Then
  273.           For k:=1 To nv Do
  274.             Begin
  275.             em:=(a[k,i]*cp)-(a[k,j]*sp);
  276.             a[k,j]:=(a[k,j]*cp)-(a[k,i]*sp);
  277.             a[k,i]:=em;
  278.             End;
  279.         End;
  280.       End;
  281.     End;
  282.   For i:=1 To nv Do
  283.     For j:=1 To nf Do a[i,j]:=a[i,j]*c[i];
  284.   For j:=1 To nf Do
  285.     Begin
  286.     c[j]:=0;
  287.     For k:=1 To nv Do c[j]:=c[j]+Sqr(a[k,j]);;
  288.     End;
  289.   em:=nv;
  290.   For j:=1 To nv Do c[j]:=c[j]/nv;
  291.   ClrScr;
  292.   For i:=1 To nf Do
  293.     Begin
  294.     Writeln(ofile); Writeln(ofile,'Rotated Factor Pattern For Factor: ',i:2);
  295.     For j:=1 To dv Do
  296.       Begin
  297.       Writeln(ofile,'   ',varn[j]:8,': ',a[i,j]:12);
  298.       End;
  299.     Writeln(ofile,'Factor Accounts For ',(c[i]*100):8:5,'% of Variance');
  300.     wait(ot);
  301.     End;
  302.   End; (* of rotate *)
  303.  
  304. Begin (* main *)
  305.   openfiles(dfile, ofile, nv, dv, ot);
  306.   Write('Number of Factors to Extract (<=number of variables)? ');
  307.   Readln(nf); If (nf=dv) Then rn:=0 Else rn:=1;
  308.   selectvar(sel, varn, dv);
  309.   (* read correlation matrix matrix *)
  310.   getcor(cor, x, y, nc, sel, nv, dv, dfile);
  311.   (* full rank analysis *)
  312.   If(rn=0) Then
  313.     Begin
  314.     For i:=1 To dv Do For j:=1 To dv Do a[i,j]:=cor[i,j];
  315.     (* invert *)
  316.     minv(a,temp,dv);
  317.     temp:=-((nc-1.0)-((1.0/6.0)*(2.0*nc+5.0)))*Ln(temp);
  318.     ClrScr;
  319.     Writeln(ofile);
  320.     Write(ofile,'Chi-Square for Sphericity = ',temp:10:3);
  321.     Writeln(ofile,'  degrees of freedom = ',((Sqr(dv)-dv)/2):6:0);
  322.     wait(ot);
  323.     For j:=1 To dv Do w[j]:=1.0-(1.0/a[j,j]);
  324.     End;
  325.   eigen(cor, vc, e, dv, nf, ot);
  326.   temp:=0.0;
  327.   For j:=1 To nf Do
  328.     Begin
  329.     u[j]:=(e[j]/dv)*100; temp:=temp+u[j]; v[j]:=temp;
  330.     For k:=1 To dv Do
  331.       Begin
  332.       a[k,j]:=cor[k,j]*(1/Sqr(e[j]));; cor[k,j]:=cor[k,j]*Sqr(e[j]);
  333.       End;
  334.     End;
  335.   ClrScr; Writeln(ofile);
  336.   Writeln(ofile,'Factor        Eigen          %Variance     %Cumulative');
  337.   Writeln(ofile);
  338.   For i:=1 To nf Do Writeln(ofile,'   ',i:2,'         ',e[i]:8:5,
  339.    '       ',u[i]:7:3,'       ',v[i]:7:3);
  340.   Writeln(ofile);
  341.   wait(ot);
  342.   For i:=1 To nf Do
  343.     Begin
  344.     Writeln(ofile); Writeln(ofile,'Factor Pattern For Factor: ',i:2);
  345.     For j:=1 To dv Do
  346.       Begin
  347.       Writeln(ofile,'   ',varn[j]:8,': ',cor[i,j]:12);
  348.       End;
  349.     wait(ot);
  350.     End;
  351.   Writeln(ofile); Write(ofile,'Variable      Communality');
  352.   If(rn=0) Then Writeln(ofile,'     R-Square') Else Writeln(ofile);
  353.    For j:=1 To dv Do
  354.     Begin
  355.     temp:=0.0;
  356.     For k:=1 To nf Do temp:=temp+Sqr(cor[j,k]);
  357.     Write(ofile,'   ',varn[j]:8,'    ',temp:10:7);
  358.     If(rn=0) Then Writeln(ofile,'     ',w[j]:10:7) Else Writeln(ofile);
  359.     End;
  360.   wait(ot);
  361.   For i:=1 To nf Do
  362.     Begin
  363.     Writeln(ofile); Writeln(ofile,'Score Coefficients For Factor: ',i:2);
  364.      For j:=1 To dv Do
  365.       Begin
  366.       Writeln(ofile,'   ',varn[j]:8,': ',a[i,j]:12);
  367.       End;
  368.     wait(ot);
  369.     End;
  370.   If(nf>1) Then
  371.     Begin
  372.     Write('What type of Factor Rotation (0=none,1=varimax,2=quartimax)?');
  373.     Readln(i); i:=i-1; If(i>=0) Then rotate(cor, nf, dv, i, ot);
  374.     End;
  375.    Close(ofile); Close(dfile);
  376.   End.
  377. ▌~═W╓08)■
  378. 8
  379.  !╓■
  380. 8■0T])╪)╪
  381.  T]╪)╪_╪╦y╖╚|ç╔▌ßßσ}&#9Nü8wδ!DφB9∙δσφ░δß+Oφ╕δ#∙▌Θ>├$▌ß═yWß═j    _ßσ}ô8(║8