home *** CD-ROM | disk | FTP | other *** search
/ Oakland CPM Archive / oakcpm.iso / sigm / vol243 / lstsqr.pqs / LSTSQR.PAS
Encoding:
Pascal/Delphi Source File  |  1986-02-10  |  4.0 KB  |  112 lines
  1. { PROGRAM AUTHOR: Mark Aldon Weiss  PROGRAM DONATED TO PUBLIC DOMAIN }
  2.  
  3. CONST
  4.  
  5. MaxNumPts = 55;
  6.  
  7.  
  8.  
  9. VAR
  10.  
  11. x,y: Array [1..MaxNumPts] of Real;   i,numpts: 1..MaxNumPts;
  12.  
  13. understood, another, StillErrors: Boolean;    ch: Char;
  14.  
  15. xavg,yavg,varx,vary,covar,sumxy,sumxsqr,slope,int,sigma,devslope,devint: Real;
  16.  
  17.  
  18. BEGIN  { M A I N    P R O G R A M }
  19. Writeln;
  20. Writeln(' This program performs a linear least squares fit.  All input and');
  21. Writeln(' output is to the terminal.  You need not turn on the printer.  You');
  22. Writeln(' should keep paper and pencil handy to jot down the results.');
  23. Writeln(' You are allowed a maximum of ',MaxNumPts:3,' points per data set.  You');
  24. Writeln(' need change only one line in the source code to accomodate more.');
  25. Writeln;
  26. understood := TRUE;
  27. REPEAT
  28.   Writeln(' YOU SHOULD ENTER YOUR DATA IN THE FOLLOWING WAY:');
  29.   Writeln(' 1.  Type your first x value; type one or more spaces; type your y');
  30.   Writeln('          value that goes with this x.  Hit return.');
  31.   Writeln(' 2.  Repeat this procedure for all your (x,y) pairs EXCEPT FOR THE');
  32.   Writeln('          LAST (x,y) pair.  FOR THE LAST PAIR, see 3. below.');
  33.   Writeln(' 3.  For your last (x,y) pair, type the x; type one or more spaces;');
  34.   Writeln('          type the y; type a * with ONE space between the y value');
  35.   Writeln('          and the *.  Hit return.');
  36.   Writeln(#7);
  37.   Write(' Did you read these instructions carefully?     ');  Readln(ch);
  38.   Writeln;
  39.   IF ch IN ['y','Y'] THEN understood := TRUE ELSE understood := FALSE
  40. UNTIL understood;
  41. Writeln;
  42. Writeln(' Okay, ENTER YOUR DATA AS INSTRUCTED ABOVE [you will be given a');
  43. Writeln(' chance to correct errors after complete entry of all your data]:');
  44. REPEAT
  45.   Writeln;  Writeln(' ENTER DATA NOW . . .');
  46.   Writeln;
  47.   i := 0;
  48.   REPEAT
  49.     i := i + 1;
  50.     Readln( x[i], y[i], ch);
  51.   UNTIL ch = '*';
  52.   numpts := i;
  53.   Writeln;
  54.   Writeln(' These are your data as received:');
  55.   Writeln;
  56.   FOR i := 1 to numpts DO Writeln(i:3,'.)    x = ',x[i],'   y = ',y[i]);
  57.   Writeln;
  58.   Write(' Are there any errors?   ');  Readln(ch);   Writeln;
  59.   IF ch IN ['y','Y'] THEN
  60.      Begin
  61.        Writeln(' Begin by correcting your first error.');  Writeln;
  62.        StillErrors := TRUE;
  63.        WHILE StillErrors DO
  64.           Begin
  65.           Writeln(' Type the following (where the <> mean to strike a key:');
  66.           Write(' data point number  <space>  x  <space>  y  <return>  --->  ');
  67.           Readln(i,x[i],y[i]);
  68.           Writeln;
  69.           Write(' Any more errors?   ');  Readln(ch);
  70.           IF ch IN ['y','Y'] THEN StillErrors := TRUE ELSE StillErrors := FALSE;
  71.           Writeln
  72.           End
  73.      End;
  74.   xavg := 0;  yavg := 0;  sumxy := 0;  varx := 0;  vary := 0;  covar := 0;
  75.   sumxsqr := 0;
  76.   FOR  i := 1 to numpts DO
  77.      Begin
  78.      xavg := xavg + x[i];
  79.      yavg := yavg + y[i];
  80.      sumxy := sumxy +  x[i] * y[i];
  81.      sumxsqr := sumxsqr + SQR( x[i] )
  82.      End;
  83.   xavg := xavg / numpts;    yavg := yavg / numpts;
  84.   FOR  i := 1 to numpts DO
  85.      Begin
  86.      varx := varx + SQR( x[i] - xavg );
  87.      vary := vary + SQR( y[i] - yavg )
  88.      End;
  89.   varx := varx / numpts;    vary := vary / numpts;
  90.   covar := sumxy / numpts - ( xavg * yavg );
  91.   slope := covar / varx;
  92.   int := yavg - slope * xavg;
  93.   sigma := SQRT( numpts/(numpts-2) * (varx*vary - SQR(covar)) / varx );
  94.   devslope := sigma / SQRT( numpts * varx );
  95.   devint := sigma * SQRT( sumxsqr / ( SQR(numpts) * varx ) );
  96.   Writeln;
  97.   Writeln('          slope = ',   slope,'            intercept = ',   int);
  98.   Writeln(' st. dev. slope = ',devslope,'   st. dev. intercept = ',devint);
  99.   Writeln;
  100.   Write(' the correlation coefficient is');
  101.   Writeln( covar / SQRT(varx * vary) );    Writeln;
  102.   Write(' Do you have another data set for analysis?   ');  Readln(ch);
  103.   IF ch IN ['y','Y'] THEN another := TRUE ELSE another := FALSE
  104. UNTIL NOT ANOTHER
  105. END.   { M A I N    P R O G R A M }
  106.  
  107.  
  108.  
  109.  
  110.  
  111.  
  112.