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 / CORREL.PZS / CORREL.ÐAS
Text File  |  2000-06-30  |  6KB  |  198 lines
  1. (* Multivariate Analysis Package - Correlation and Covariance 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 Correl(Input,Output);
  7. Const
  8.    N=50;
  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.   cvj, muj, nc, miss, vars, mean, stdev, sm : RVEC;
  19.   varn : S8;
  20.   mu, cv : NBYN;
  21.   nv, i, j, k, wgt : Integer;
  22.   ot, dv, wt : SUBS;
  23.   vj, vw, vr, wg2, wht, corr : Real;
  24.   wflg : Boolean;
  25.  
  26. Procedure openfiles(Var dfile, ofile:Text; Var ot:SUBS);
  27. Var
  28.   dfl, ofl:String[12];
  29. Begin
  30.   ClrScr; Writeln(' *** CORREL: CORRELATION AND COVARIANCE MATRIX ***');
  31.   Writeln; Write('Name of the data file? ');
  32.   Readln(dfl); Assign(dfile,dfl); Reset(dfile);
  33.   Write('Name of the output file? ');
  34.   Readln(ofl); Assign(ofile,ofl); Rewrite(ofile);
  35.   ot:= 0;
  36.   If (ofl='CON:') Or (ofl='con:') Then ot:=1;
  37.   If (ofl='LST:') Or (ofl='lst:') Then ot:=2;
  38.   If (ot=2) Then
  39.     Begin
  40.     Writeln(ofile,'Multivariate Analysis Package (1.6) - ',
  41.       'Program: CORREL, Datafile: ',dfl); Writeln(ofile);
  42.     End;
  43.   End; (* Of openfiles *)
  44.  
  45. Procedure wait(ot:Integer);
  46. Begin
  47.   If ot=1 Then Begin
  48.     Write('- Press any key to continue -'); While Not KeyPressed Do; ClrScr;
  49.     End;
  50.   End; (* of wait *)
  51.  
  52. Procedure selcvar(Var sel:IVEC;Var varn:S8;Var miss:RVEC;
  53.                   Var nv:Integer; Var dv, wt:SUBS);
  54. Var
  55.   cfile:Text;
  56.   cfl:String[12];
  57.   i,j,f:Integer;
  58.   mis:Real;
  59.   van:String[8];
  60. Begin
  61.   Write('Name of the codebook file (or NONE)? '); Readln(cfl);
  62.   If (cfl<>'NONE') And (cfl<>'none') Then f:=1 Else f:=0;
  63.   If f=1 Then Begin Assign(cfile,cfl); Reset(cfile); End;
  64.   Writeln;
  65.   Write('How many variables in data file? '); Readln(nv);
  66.   Write('Number of variables to use in CORREL? '); Readln(dv);
  67.   For i := 1 To dv Do
  68.     Begin
  69.     Write('Column number for variable ',i,'? '); Readln(sel[i]);
  70.     Str(sel[i]:3,varn[i]); miss[i]:=-1E37;
  71.     End;
  72.   Write('Of these Column numbers which is weight (0=none)? '); Readln(wt);
  73.   If f=1 Then Begin
  74.     For j:=1 to nv Do Begin
  75.       mis:=-1E37;
  76.       Readln(cfile,f,van,mis);
  77.       For i:=1 to dv Do
  78.         If f=sel[i] Then Begin
  79.           varn[i]:=van; miss[i]:=mis;
  80.           Writeln('Col: ',sel[i],'  Name: ',varn[i],' Missing: ',miss[i]:6);
  81.           End;
  82.       End;
  83.     Close(cfile);
  84.     wait(1);
  85.     End;
  86.   End; (* Of selcvar *)
  87.  
  88. Procedure getcase(Var vars:RVEC;Var sel:IVEC; nv:Integer; dv:SUBS;
  89.                   Var dfile:Text);
  90. Var
  91.   i, j:SUBS;
  92.   x:Real;
  93. Begin
  94.   For i := 1 To nv Do
  95.     Begin
  96.     Read(dfile,x); For j:=1 to dv Do If (sel[j]=i) Then vars[j]:=x;
  97.     End;
  98.   End; (* Of getcase *)
  99.  
  100. Begin (* main *)
  101.   openfiles(dfile, ofile, ot);
  102.   selcvar(sel, varn, miss, nv, dv, wt);
  103.   (* initialize *)
  104.   wgt:=0; wflg:=false;
  105.   FillChar(nc,6*N,0); FillChar(stdev,6*N,0); FillChar(sm,6*N,0);
  106.   FillChar(cv,6*N*N,0); FillChar(mu,6*N*N,0);
  107.   For i := 1 To dv Do
  108.     If wt=sel[i] Then Begin wgt:=i; wflg:=true; End;
  109.   (* accumulate sscp *)
  110.   Writeln;
  111.   k:=0;
  112.  
  113.   While Not EOF(dfile) Do
  114.     Begin
  115.     k:=k+1;
  116.     getcase(vars, sel, nv, dv, dfile);
  117.     If Frac(k/10)=0.0 Then Write('+');
  118.     { core computational loop }
  119.     vw:=vars[wgt]; { min subscript refs }
  120.     If Not EOF(dfile) Then
  121.       If Not (wflg And (vw=miss[wgt])) Then
  122.       For j:=1 To dv Do
  123.         If vars[j]<>miss[j] Then
  124.           Begin
  125.           vj:=vars[j]; { min subscript refs }
  126.           If wflg And (j<>wgt) Then wht:=vw Else wht:=1.0;
  127.           { spicer algorithm for descriptive stats }
  128.           nc[j]:=nc[j]+wht;
  129.           If wflg Then Begin { avoid mult and div by 1 if no weight }
  130.             vr:=(vj-mean[j])*wht/nc[j];
  131.             stdev[j]:=stdev[j]+(nc[j]*(nc[j]-wht)/wht)*Sqr(vr);
  132.           End Else Begin
  133.             vr:=(vj-mean[j])/nc[j];
  134.             stdev[j]:=stdev[j]+(nc[j]*(nc[j]-1.0))*Sqr(vr);
  135.           End;
  136.           mean[j]:=mean[j]+vr;
  137.           { loop for covariance and pairwise deletion means stored
  138.           in same matrix as triangulars plus diagonal sm }
  139.           cvj:=cv[j]; muj:=mu[j]; {slicing to avoid subs nc*dv*dv*6 times}
  140.           For i:=j To dv Do
  141.             If vars[i]<>miss[i] Then
  142.               Begin
  143.               If wflg And (i<>wgt) Then wg2:=wht*vw Else wg2:=wht;
  144.               If i=j Then Begin sm[i]:=sm[i]+wg2; vr:=wg2/sm[i]; End
  145.               Else {i>j}
  146.                 Begin
  147.                 cv[i,j]:=cv[i,j]+wg2; vr:=wg2/cv[i,j];
  148.                 mu[i,j]:=mu[i,j]+(vars[i]-mu[i,j])*vr;
  149.                 End;
  150.               muj[i]:=muj[i]+(vj-muj[i])*vr;
  151.               cvj[i]:=cvj[i]+(vj*vars[i]-cvj[i])*vr;
  152.               End;
  153.           cv[j]:=cvj; mu[j]:=muj;
  154.           End;
  155.     End;
  156.  
  157.   (* compute And Output *)
  158.   ClrScr;
  159.   Writeln(ofile,'Means and Standard Deviations:');
  160.   Writeln(ofile);
  161.   For j:=1 To dv Do
  162.     Begin
  163.     If nc[j] > 1.0 Then
  164.       Begin
  165.       stdev[j]:=Sqrt(stdev[j]/(nc[j]-1.0));
  166.       Writeln(ofile,varn[j]:8,'  Mean:',mean[j]:13:5,
  167.        '  Standard Deviation:',stdev[j]:13:5,'  Cases:',nc[j]:9:0);
  168.       End;
  169.     End;
  170.   Writeln(ofile);
  171.   wait(ot);
  172.   Writeln(ofile,'Covariance and Correlations (pairwise deletion):');
  173.   Writeln(ofile);
  174.   For j:= 1 To dv Do
  175.     If sm[j] > 1 Then
  176.       cv[j,j]:=cv[j,j]-Sqr(mu[j,j]);
  177.   For j:= 1 To dv Do
  178.     Begin
  179.     cvj:=cv[j]; muj:=mu[j]; {slicing}
  180.     For i:=j+1 To dv Do
  181.       If cv[i,j]>0 Then cvj[i]:=cvj[i]-(muj[i]*mu[i,j]);
  182.     For i:=j+1 To dv Do
  183.       Begin
  184.       If cv[i,j]>0.0 Then
  185.         Begin
  186.         corr:=cvj[i]/Sqrt(cv[i,i]*cvj[j]);
  187.         Writeln(ofile,varn[j]:8,' with ',varn[i]:8,'   Covariance: ',
  188.           cvj[i]:13,'   Correlation:',corr:10:6);
  189.         End;
  190.       End;
  191.     End;
  192.   Close(dfile); Close(ofile);
  193.   Assign(dfile,'MAPSTAT.COM'); Execute(dfile);
  194. End.
  195.  Deviations:');
  196.   Writeln(ofile);
  197.   For j:=1 To dv Do
  198.