home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / pascal / library / dos / math / matrix / fileregr.pas next >
Encoding:
Pascal/Delphi Source File  |  1988-06-13  |  10.7 KB  |  259 lines
  1. program REGRESS_FILE;
  2.         uses MATRIXU,CRT;
  3.  
  4.      {by       Josh Mitteldorf
  5.                6800 Scotforth Rd
  6.                Philadelphia, PA 19119
  7.                70441,1305 }
  8.  
  9.  
  10.    {Performs a linear regression (best line fit) on an ASCII-formatted file.
  11.     File should be organized in columns, with only whitespace between the
  12.     columns.  You may specify at run time which column contains the dependent
  13.     variable and which contains the independent variables.  Columns which
  14.     are not specified as either are read in and skipped.  They may contain
  15.     non-numeric data without causing an error.  However, non-numeric data in
  16.     any variable cell will cause an entire row of numbers to be skipped over,
  17.     with a message to the user.
  18.  
  19.     If a row is found with more or less than the specified number of columns,
  20.     the program attempts to skip that row and move on.}
  21.  
  22. const Nmax=20;  NSQ=400;    cmax=50;  maxstr=80;
  23.  
  24. type float=double;
  25.      vector=array[1..Nmax] of float;
  26.      matrix=array[1..NSQ] of float;
  27.  
  28. const errval:float=9.87654321987654321E37;  ef=#26;   space=' ';   tab=#9;
  29.       cr=#13;  lf=#10;   crlf:array[1..2] of char=(cr,lf);
  30.  
  31. var  f                                   :text;
  32.      xxsum,xxavg,xxvar,xinv              :matrix;
  33.      xavg,xyavg,xsum,xysum,x,
  34.      xavgyavg,aw,num                     :vector;
  35.      y,ysum,yysum,yavg,yyavg,norm,
  36.      yvar,dy,dyf                         :float;
  37.      n,d,count,t1,t2,t3,cols,ycol,nmiss  :integer;
  38.      g                                   :file of float;
  39.      xcol                                :array[1..Nmax] of integer;
  40.      filename                            :string[maxstr];
  41.      vertical                            :boolean;
  42.  
  43. procedure CumRegistersToZero;
  44.           var i,j  :integer;
  45.           begin
  46.           ysum:=0; yysum:=0; n:=0;  nmiss:=0;
  47.           for i:=1 to d do begin
  48.               xsum[i]:=0;
  49.               xysum[i]:=0;
  50.               for j:=1 to d do xxsum[d*pred(i)+j]:=0;
  51.               end;
  52.           end;
  53.  
  54. procedure Accumulate;
  55.           var k,i,j   :integer;
  56.               xi      :float;
  57.           begin
  58.           inc(n);
  59.          {$ifdef DIAGNOSTIC} write(n:4,'  '); {$endif}
  60.           ysum:=ysum+y;    yysum:=yysum+sqr(y);
  61.           for i:=1 to d do begin
  62.                   xi:=x[i];
  63.                   xsum[i]:=xsum[i]+xi;
  64.                   xysum[i]:=xysum[i]+y*xi;
  65.                   for j:=i to d do xxsum[d*pred(i)+j]:=xxsum[d*pred(i)+j]+xi*x[j];
  66.                   end;                        {Upper half only, for speed}
  67.               end;
  68.  
  69. procedure SumsToAverages;
  70.           var i,j  :integer;
  71.  
  72.           begin
  73.           if (n-d-1<=0) then begin
  74.              writeln(crlf,'Not enough points for regression.');
  75.              writeln('You have ',n,' lines of valid data and ',d,' regression variables.');
  76.              halt; end;
  77.           norm:=1/n;
  78.           yavg:=ysum*norm;    yyavg:=yysum*norm;
  79.             {$ifdef DIAGNOSTIC} write('XSUM: '); vecprint(xsum,d); {$endif}
  80.           ScalarMultV(addr(xsum),norm,addr(xavg),d);
  81.             {$ifdef DIAGNOSTIC} write('XAVG: '); vecprint(xavg,d); {$endif}
  82.           ScalarMultV(addr(xavg),yavg,addr(xavgyavg),d);
  83.           ScalarMultV(addr(xysum),norm,addr(xyavg),d);
  84.           ScalarMult(addr(xxsum),norm,addr(xxavg),d);
  85.    {$ifdef DIAGNOSTIC} write('XXAVG: '); for i:=1 to d do write('[',i,']=',xxavg[{pred(1)*d+} i]:11:6,'  '); writeln; {$endif}
  86.           for i:=1 to d do for j:=succ(i) to d do xxavg[d*pred(j)+i]:=xxavg[d*pred(i)+j];
  87.                                 {Because we saved time and computed upper half before}
  88.           yvar:=yyavg-sqr(yavg);
  89.           end;
  90.  
  91. procedure ReportMissing;
  92.           begin
  93.           inc(nmiss);
  94.           write(crlf,'Missing data in line ',n + nmiss,'.');
  95.           end;
  96.  
  97. procedure Regression;
  98.           var xavg2          :matrix;
  99.               wv             :vector;
  100.               i,j            :integer;
  101.               wr             :float;
  102.               c              :array[1..cmax] of float;
  103.               ok,warned      :boolean;
  104.               ch             :char;
  105.  
  106.    procedure DistributeAndCheckData;
  107.              var i :integer;
  108.              begin
  109.              y:=c[ycol];
  110.              ok:=(y<>errval);
  111.              for i:=1 to d do begin
  112.                  x[i]:=c[xcol[i]];
  113.                  ok:=(ok and (x[i]<>errval));
  114.                  end;
  115.              repeat read(f,ch) until ((eof(f) or not (ch in [space,tab])));
  116.              if ((ch<>cr) and (ch<>ef)) then begin
  117.                 write(crlf,'Warning!  Unexpected number of columns in line ',1+n+nmiss,'.');
  118.                 ok:=false;  warned:=true; inc(nmiss);
  119.                 end
  120.              else warned:=false;
  121.              repeat read(f,ch) until ((ch=lf) or (eof(f)));
  122.              if (eof(f) and (x[1]=0) and (y=0)) then ok:=false;
  123.                 {A blank line at the end of the file reads in as a row of zeros.}
  124.              end;
  125.  
  126.  
  127.           begin {Regression}
  128.           CumRegistersToZero;
  129.           repeat
  130.              i:=1;
  131.              for i:=1 to cols do begin
  132. {$I-}            read(f,c[i]);
  133.                  if (ioresult>0) then c[i]:=errval;
  134.                  { write(c[i]:16); }
  135. {$I+}            end;
  136.              DistributeAndCheckData;
  137. {$ifdef DIAGNOSTIC}   write(y:8:5); for i:=1 to d do write(x[i]:8:5); writeln; {$endif}
  138.              if (ok) then Accumulate else if (warned) or (eof(f) and (x[1]=0) and (y=0)) then else ReportMissing;
  139.           until (ioresult>0) or (eof(f));
  140.           SumsToAverages;
  141.           Diadic(addr(xavg),addr(xavg),addr(xavg2),d);
  142.           MatDif(addr(xxavg),addr(xavg2),addr(xxvar),d);
  143.  
  144.           VecDif(addr(xyavg),addr(xavgyavg),addr(num),d);
  145.             {$ifdef DIAGNOSTIC}
  146.              writeln('Y AVG: ',yavg:11);
  147.              write(crlf,'XY AVG:'); for i:=1 to d do write(xyavg[i]:11,' ');
  148.              write(crlf,'XAVYAV:'); for i:=1 to d do write(xavgyavg[i]:11,' ');
  149.              write(crlf,'NUM: '); for i:=1 to d do write(num[i]:11,' ');
  150.              write(crlf,'DEN: '); for i:=1 to d do begin
  151.                  for j:=1 to d do write(xxvar[i,j]:11,' '); write(crlf,'     ');
  152.                  end; {$endif}
  153.           Solve(addr(xxvar),addr(num),addr(aw),d);
  154.           Prod(addr(xxvar),addr(aw),addr(wv),d);  {wv:=σ²(x) · a}
  155.           dy:=sqrt(yvar - dot(addr(wv),addr(aw),d));   {δy²= σ²(y) - a²σ²(x)}
  156.           dyf:=dy/sqrt(n-d-1);
  157. {          for i:=1 to d do write(dyf/sqrt(xxavg[d*pred(i)+i]-xavg2[d*pred(i)+i]):9:5); writeln;  (This is my own measure.)  }
  158.           Inverse(addr(xxvar),addr(xinv),d);
  159.           end;
  160.  
  161. procedure Report;
  162.           var i    :integer;
  163.           begin
  164.           clrscr;
  165.           write('INDEPENDENT: col ',ycol,'  DEPENDENT:'); for i:=1 to d do write(' ',xcol[i]);
  166.           writeln(crlf,n,' point regression with ',d,' variables.  ',n-d-1,' degrees of freedom.');
  167.           writeln('R² = ',Dot(addr(aw),addr(num),d)/yvar:9:5);
  168.           writeln('δY (RMS deviance from line) = ',dy:9:5,',');
  169.           writeln('Compare σ(Y) = ',sqrt(yvar):9:5);
  170.           if (vertical) then begin
  171.             writeln(crlf,'   COL    COEFFICIENT      STD ERROR     STD DEVIATION    PROP. OF VARIANCE');
  172.             writeln(     '--------  -----------      ---------     -------------    -----------------');
  173.             writeln(     ' const  ',yavg-dot(addr(xavg),addr(aw),d):13:7,dy:15:7);
  174.             for i:=1 to d do
  175.               writeln(xcol[i]:5,aw[i]:16:7,dyf*sqrt(xinv[d*pred(i)+i]):15:7,sqrt(xxvar[pred(i)*d+i]):18:5,
  176.                   xxvar[pred(i)*d+i]*sqr(aw[i])/yvar:21:5);
  177.             end
  178.           else begin
  179.             writeln('B = ',yavg-dot(addr(xavg),addr(aw),d):9:5,'  (Std error=δY)');
  180.             writeln('Coefficient vector A:');
  181.             for i:=1 to d do write(aw[i]:9:5); writeln;
  182.             writeln('Std errors of coefficients:');
  183.             for i:=1 to d do write(dyf*sqrt(xinv[d*pred(i)+i]):9:5); writeln;
  184.             writeln('Std Deviations σ(X):');
  185.             for i:=1 to d do write(sqrt(xxvar[pred(i)*d+i]):9:5); writeln;
  186.             writeln('Proportion of Y''s variance attributable to:');
  187.             for i:=1 to d do write(xxvar[pred(i)*d+i]*sqr(aw[i])/yvar:9:5); writeln;
  188.             end;
  189.           end;
  190.  
  191.  
  192. procedure Parse(xstr :string);
  193.           var i,j,k,ta  :integer;
  194.               substr    :string[80];
  195.               v         :float;
  196.           begin
  197.           for i:=succ(length(xstr)) to maxstr do xstr[i]:=' ';
  198.           i:=0; k:=1;
  199.           while (i<length(xstr)) do begin
  200.                 substr:=''; j:=1;
  201.                 repeat inc(i) until (xstr[i] in ['0'..'9']);
  202.                 repeat substr[j]:=xstr[i]; inc(i); inc(j); until not (xstr[i] in ['0'..'9']);
  203.                 substr[0]:=char(pred(j));
  204.                 val(substr,v,ta);
  205.                 xcol[k]:=round(v); inc(k);
  206.                 end;
  207.           d:=pred(k);
  208.           write('There are ',d,' variables in cols ');
  209.           for k:=1 to d do write(xcol[k],' ');
  210.           end;
  211.  
  212. procedure Input;
  213.           var ki  :char;
  214.              xstr :string[80];
  215.              i    :integer;
  216.              ok   :boolean;
  217.           begin
  218.           clrscr;
  219.           writeln('ASCII FILE -> Regression',crlf);
  220. {$I-}     repeat
  221.              write('What is the name of the file containing the data? ');
  222.                 readln(filename); assign(f,filename); reset(f);
  223.           until (ioresult=0);
  224. {$I+}     write('Total number of ');
  225.            textattr:=112; write('columns'); textattr:=7;
  226.            write(' in file: '); readln(cols);
  227.           repeat
  228.             write('Which column contains the ');
  229.              textattr:=112; write('dependent'); textattr:=7;
  230.              write(' variable? '); readln(ycol);
  231.              ok:=((ycol<=cols) and (ycol>0));
  232.           until ok;
  233.           repeat
  234.              write('Which columns contain the ');
  235.               textattr:=112; write('independent'); textattr:=7;
  236.               write(' variables? '); readln(xstr);
  237.              Parse(xstr);
  238.              ok:=true;
  239.              for i:=1 to d do ok:=(ok and ((xcol[i]<=cols) and (xcol[i]>0)));
  240.              if ok then begin
  241.                 write(crlf,'  Correct? '); readln(ki);
  242.                 end
  243.              else begin ki:='N'; writeln; end;
  244.           until upcase(ki) in ['Y','y'];
  245.           textattr:=112; write('H'); textattr:=7;
  246.           write('orizontal or ');
  247.           textattr:=112; write('V'); textattr:=7;
  248.           write('ertical output form? ');
  249.           repeat ki:=upcase(readkey) until ki in ['H','V']; writeln(ki);
  250.           if (ki='V') then vertical:=true else vertical:=false;
  251.           {window(1,8,80,25); clrscr; window(1,1,80,25); gotoxy(1,8);}
  252.           write('Computing...');
  253.           end;
  254.  
  255. begin
  256. Input;
  257. Regression;
  258. Report;
  259. end.