home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / TURBOPAS / PAS_ENG.ZIP / LEAST1.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-07-19  |  3.6 KB  |  180 lines
  1. program least1;  { --> 191 }
  2. { Pascal Program to perform a liner least-squares fit using a parabolic }
  3. { curve. Aeperate procedure PLOT needed }
  4.  
  5. const maxr = 20;
  6.  maxc = 3;
  7.  
  8. type ary = array[1..maxr] of real;
  9.  arys = array[1..maxc] of real;
  10.  ary2s = array[1..maxc,1..maxc] of real;
  11.  
  12. var x,y,y_calc : ary;
  13.  coef  : arys;
  14.  nrow,ncol : integer;
  15.  correl_coef : real;
  16.  
  17.  
  18.  
  19. procedure get_data(var x,y: ary;
  20.      var nrow: integer);
  21. { get values for n and arrays x,y }
  22.  
  23. var i : integer;
  24.  
  25. begin
  26.   nrow:=9;
  27.   writeln;
  28.   for i:=1 to nrow do x[i]:=i;
  29.   y[1]:=2.07; y[2]:=8.6;
  30.   y[3]:=14.42; y[4]:=15.8;
  31.   y[5]:=18.92; y[6]:=17.96;
  32.   y[7]:=12.98; y[8]:=6.45;
  33.   y[9]:=0.27;
  34. end;  { procedure get_data }
  35.  
  36. procedure write_data;
  37. { print out the answers }
  38. var i : integer;
  39. begin
  40.   writeln;
  41.   writeln('  I      X       Y       YCALC');
  42.   for i:=1 to nrow do
  43.     writeln(i:3,x[i]:8:1,y[i]:9:2,y_calc[i]:9:2);
  44.   writeln; writeln(' Coefficients ');
  45.   for i:=1 to ncol do
  46.     writeln(coef[i]:8:4);
  47.   writeln;
  48.   writeln('Correlation coefficient is ',correl_coef:8:5)
  49. end;  { write_data }
  50.  
  51. procedure solve(a: ary2s;
  52.   y: arys;
  53.   var coef: arys;
  54.       nrow: integer;
  55.  var error: boolean);
  56.  
  57. var b : ary2s;
  58.  i,j : integer;
  59.  det : real;
  60.  
  61.  
  62. function deter(a: ary2s): real;
  63.  { calculate the determinant of a 3-by-3matrix }
  64. begin
  65.   deter:=a[1,1]*(a[2,2]*a[3,3]-a[3,2]*a[2,3])
  66.  -a[1,2]*(a[2,1]*a[3,3]-a[3,1]*a[2,3])
  67.  +a[1,3]*(a[2,1]*a[3,2]-a[3,1]*a[2,2])
  68. end;
  69.  
  70.  
  71.  
  72. procedure setup(var b : ary2s;
  73.   var coef: arys;
  74.   j : integer);
  75.  
  76. var i : integer;
  77.  
  78. begin { setup }
  79.   for i:=1 to nrow do
  80.     begin
  81.       b[i,j]:=y[i];
  82.       if j>1 then b[i,j-1]:=a[i,j-1]
  83.     end;
  84.   coef[j]:=deter(b)/det
  85. end;  { setup }
  86.  
  87. begin { procedure solve }
  88.   error:=false;
  89.   for i:=1 to nrow do
  90.     for j:=1 to nrow do
  91.       b[i,j]:=a[i,j];
  92.   det:=deter(b);
  93.   if det=0.0 then
  94.     begin
  95.       error:=true;
  96.       writeln(chr(7),'ERROR: matrix is singular')
  97.     end
  98.   else
  99.     begin
  100.       setup(b,coef,1);
  101.       setup(b,coef,2);
  102.       setup(b,coef,3)
  103.     end  { esle }
  104. end;  { procedure solve }
  105.  
  106.  
  107. procedure linfit(x,y: ary;
  108.    var y_calc: ary;
  109.    var coef:   arys;
  110.        nrow:   integer;
  111.    var ncol:   integer);
  112.  
  113. { least squares fit to a parabola }
  114. { nrow sets of x and y pair points }
  115.  
  116. var a  : ary2s;
  117.  g  : arys;
  118.  i  : integer;
  119.  error  : boolean;
  120.  
  121.  sum_x,sum_y,sum_xy,sum_x2,
  122.  sum_y2,xi,yi,sxy,syy,
  123.  sxx,sum_x3,sum_x4,sum_2y,
  124.  denom,srs,x2 : real;
  125.  
  126. begin   { linfit }
  127.   ncol:=3; { polynomial terms }
  128.   sum_x:=0.0;
  129.   sum_y:=0.0;
  130.   sum_xy:=0.0;
  131.   sum_x2:=0.0;
  132.   sum_y2:=0.0;
  133.   sum_x3:=0.0;
  134.   sum_x4:=0.0;
  135.   sum_2y:=0.0;
  136.   for i:=1 to nrow do
  137.     begin
  138.       xi:=x[i];
  139.       yi:=y[i];
  140.       x2:=xi*xi;
  141.       sum_x:=sum_x+xi;
  142.       sum_y:=sum_y+yi;
  143.       sum_xy:=sum_xy+xi*yi;
  144.       sum_x2:=sum_x2+x2;
  145.       sum_y2:=sum_y2+yi*yi;
  146.       sum_x3:=sum_x3+xi*x2;
  147.       sum_x4:=sum_x4+x2*x2;
  148.       sum_2y:=sum_2y+x2*yi
  149.     end;
  150.   a[1,1]:=nrow;
  151.   a[2,1]:=sum_x; a[1,2]:=sum_x;
  152.   a[3,1]:=sum_x2; a[1,3]:=sum_x2;
  153.   a[2,2]:=sum_x2; a[3,2]:=sum_x3;
  154.   a[2,3]:=sum_x3; a[3,3]:=sum_x4;
  155.   g[1]:=sum_y;
  156.   g[2]:=sum_xy;
  157.   g[3]:=sum_2y;
  158.   solve(a,g,coef,ncol,error);
  159.   srs:=0.0;
  160.   for i:=1 to nrow do
  161.     begin
  162.       y_calc[i]:=coef[1]+coef[2]*x[i]+coef[3]*sqr(x[i]);
  163.       srs:=srs+sqr(y[i]-y_calc[i])
  164.     end;
  165.   correl_coef:=sqrt(1.0-srs/(sum_y2-sqr(sum_y)/nrow))
  166. end;  { linfit }
  167.  
  168. { external procedure plot(x,y,y_calc: ary; nrow: integer);
  169. }
  170.  
  171. {$I PLOT.LIB }  { get ptocedure PLOT }
  172.  
  173. begin  { MAIN program }
  174.   clrscr;
  175.   get_data(x,y,nrow);
  176.   linfit(x,y,y_calc,coef,nrow,ncol);
  177.   write_data;
  178.   plot(x,y,y_calc,nrow)
  179. end.  { MAIN }
  180.