home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / TURBOPAS / PAS_ENG.ZIP / CFIT4.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-07-18  |  1.9 KB  |  105 lines
  1. program cfit4;  {164}
  2.  
  3. { plot service included }
  4. { Pascal program to perform a linear least-squares fit }
  5.  
  6. const max = 20;
  7.  
  8. type index = 1..max;
  9.  ary = array[index] of real;
  10.  
  11. var x,y,y_calc : ary;
  12.  n  : integer;
  13.  first,done : boolean;
  14.  a,b,correl_coef,
  15.  sigma_a,sigma_b,
  16.  see,seed : real;
  17.  
  18.  
  19. function random(dummy: integer): real;
  20. { random number 0-1 }
  21. { define seed=4.0 as global }
  22.  
  23. const pi = 3.14159;
  24.  
  25. var x : real;
  26.  i : integer;
  27.  
  28. begin { RANDOM }
  29.   x:=seed+pi;
  30.   x:=exp(5.0*ln(x));
  31.   seed:=x-trunc(x);
  32.   random:=seed
  33. end; { RANDOM }
  34.  
  35.  
  36.  
  37. procedure get_data(var x,y: ary;
  38.      var n: integer);
  39. { get values for n and arrays x,y }
  40. { y is randomly scattered about a straight line }
  41.  
  42. const a = 2.0;
  43.  b = 5.0;
  44.  
  45. var i,j : integer;
  46.  fudge : real;
  47.  
  48. begin
  49.   write('Fudge? ');
  50.   readln(fudge);
  51.   if fudge<0.0 then done:=true
  52.   else
  53.     begin
  54.       repeat
  55.  write('How many points? ');
  56.  readln(n)
  57.       until (n>2) and (n<=max);
  58.       if first then first:=false else clrscr;
  59.       for i:=1 to n do
  60.  begin
  61.    j:=n+1-i;
  62.    x[i]:=j;
  63.    y[i]:=(a+b*j)*(1.0+(2.0*random(0)-1.0)*fudge)
  64.       end { for-loop }
  65.     end  { if }
  66. end;  { procedure get_data }
  67.  
  68.  
  69. procedure write_data;
  70. { print out the answers }
  71. var i,j : integer;
  72.  
  73. begin
  74.   writeln;
  75.   writeln('  I      X       Y      YCALC');
  76.   for i:=1 to n do
  77.     writeln(i:3,x[i]:8:1,y[i]:9:2,y_calc[i]:9:2);
  78.   writeln;
  79.   writeln('Intercept is ',a:8:3,', sigma is ',sigma_a:8:3);
  80.   writeln('    Slope is ',b:8:2,', sigma is ',sigma_b:8:3);
  81.   writeln;
  82.   writeln('Correlation coeffizient is ',correl_coef:7:4);
  83.   for j:=1 to 40 do for i:=1 to 10000 do; clrscr
  84. end;  { write_data }
  85.  
  86. {$I LINFIT2.LIB }
  87.  
  88. {$I PLOT.LIB }
  89.  
  90. begin { MAIN program }
  91.   seed:=4.0;
  92.   clrscr;
  93.   first:=true;
  94.   done:=false;
  95.   repeat
  96.     get_data(x,y,n);
  97.     if not done then
  98.       begin
  99.  linfit(x,y,y_calc,a,b,n);
  100.  write_data;
  101.  plot(x,y,y_calc,n);
  102.     end
  103.   until done
  104. end.
  105.