home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / TURBOPAS / PAS_ENG.ZIP / LEAST2.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-07-19  |  3.9 KB  |  163 lines
  1. program least2;  { --> 203 }
  2. { Pascal program to perform a linear least-squares fit }
  3. { with Gauss-Jordan routine }
  4. { Sperate modules needed:
  5.    GAUSSJ,
  6.    PLOT }
  7.  
  8.  
  9. const maxr = 20;  { data prints }
  10.  maxc = 4;  { polynomial terms }
  11.  
  12. type ary = array[1..maxr] of real;
  13.  arys = array[1..maxc] of real;
  14.  ary2 = array[1..maxr,1..maxc] of real;
  15.  ary2s = array[1..maxc,1..maxc] of real;
  16.  
  17. var x,y,y_calc : ary;
  18.  resid  : ary;
  19.  coef,sig : arys;
  20.  nrow,ncol : integer;
  21.  correl_coef : real;
  22.  
  23.  
  24.  
  25. procedure get_data(var x: ary;  { independant variable }
  26.      var y: ary;  { dependant variable }
  27.      var nrow: integer); { length of vectors }
  28. { get values for n and arrays x,y }
  29.  
  30. var i : integer;
  31.  
  32. begin
  33.   nrow:=9;
  34.   for i:=1 to nrow do x[i]:=i;
  35.   y[1]:=2.07; y[2]:=8.6;
  36.   y[3]:=14.42; y[4]:=15.8;
  37.   y[5]:=18.92; y[6]:=17.96;
  38.   y[7]:=12.98; y[8]:=6.45;
  39.   y[9]:=0.27;
  40. end;  { proceddure get data }
  41.  
  42. procedure write_data;
  43. { print out the answers }
  44. var i : integer;
  45. begin
  46.   writeln;
  47.   writeln;
  48.   writeln('  I      X      Y        YCALC   RESID');
  49.   for i:=1 to nrow do
  50.     writeln(i:3,x[i]:8:1,y[i]:9:2,y_calc[i]:9:2,resid[i]:9:2);
  51.   writeln; writeln(' Coefficients errors ');
  52.   writeln(coef[1],' ',sig[1],' constant term');
  53.   for i:=2 to ncol do
  54.     writeln(coef[i],' ',sig[i]);  { other terms }
  55.   writeln;
  56.   writeln('Correlation coefficient is ',correl_coef:8:5)
  57. end;  { write_data }
  58.  
  59. procedure square(x: ary2;
  60.    y: ary;
  61.       var a: ary2s;
  62.       var g: arys;
  63.   nrow,ncol: integer);
  64. { matrix multiplication routine }
  65. { a= transpose x times x }
  66. { g= y times x }
  67.  
  68. var i,k,l : integer;
  69.  
  70. begin  { square }
  71.   for k:=1 to ncol do
  72.     begin
  73.       for l:=1 to k do
  74.  begin
  75.    a[k,l]:=0.0;
  76.    for i:=1 to nrow do
  77.      begin
  78.        a[k,l]:=a[k,l]+x[i,l]*x[i,k];
  79.        if k<>l then a[l,k]:=a[k,l]
  80.      end
  81.  end;  { l-loop }
  82.       g[k]:=0.0;
  83.       for i:=1 to nrow do
  84.  g[k]:=g[k]+y[i]*x[i,k]
  85.       end { k-loop }
  86. end;  { SQUARE }
  87.  
  88. {external procedure gaussj(var b: ary2s;
  89.          y: arys;
  90.      var coef: arys;
  91.      ncol:  integer;
  92.      var error: boolean);
  93. }
  94. {$I GAUSSJ.LIB }
  95.  
  96. procedure linfit(x,  { independant variable }
  97.    y: ary; { dependent variable }
  98.    var y_calc: ary; { calculated dep. variable }
  99.    var resid:  ary; { array of residuals }
  100.    var coef:   arys; { coefficients }
  101.    var sig:    arys; { error on coefficients }
  102.    nrow:       integer; { length of array }
  103.    var ncol:   integer); { number of terms }
  104.  
  105. { least squares fit to nrow sets of x and y pairs of points }
  106. { Seperate procedures needed:
  107.  SQUARE -> form square coefficient matrix
  108.  GAUSSJ -> Gauss-Jordan elimination }
  109.  
  110. var xmatr  : ary2;  { data matrix }
  111.  a  : ary2s; { coefficient matrix }
  112.  g  : arys;  { constant vector }
  113.  error  : boolean;
  114.  i,j,nm  : integer;
  115.  xi,yi,yc,srs,see,
  116.  sum_y,sum_y2 : real;
  117.  
  118. begin  { procedure linfit }
  119.   ncol:=3; { number of terms }
  120.   for i:=1 to nrow do
  121.     begin  { setup matrix }
  122.       xi:=x[i];
  123.       xmatr[i,1]:=1.0; { first column }
  124.       xmatr[i,2]:=xi; { second column }
  125.       xmatr[i,3]:=xi*xi { third column }
  126.     end;
  127.   square(xmatr,y,a,g,nrow,ncol);
  128.   gaussj(a,g,coef,ncol,error);
  129.   sum_y:=0.0;
  130.   sum_y2:=0.0;
  131.   srs:=0.0;
  132.   for i:=1 to nrow do
  133.     begin
  134.       yi:=y[i];
  135.       yc:=0.0;
  136.       for j:=1 to ncol do
  137.  yc:=yc+coef[j]*xmatr[i,j];
  138.       y_calc[i]:=yc;
  139.       resid[i]:=yc-yi;
  140.       srs:=srs+sqr(resid[i]);
  141.       sum_y:=sum_y+yi;
  142.       sum_y2:=sum_y2+yi*yi
  143.     end;
  144.   correl_coef:=sqrt(1.0-srs/(sum_y2-sqr(sum_y)/nrow));
  145.   if nrow=ncol then nm:=1
  146.   else nm:=nrow-ncol;
  147.   see:=sqrt(srs/nm);
  148.   for i:=1 to ncol do  { errors on solution }
  149.     sig[i]:=see*sqrt(a[i,i])
  150. end; { linfit }
  151.  
  152. {external procedure plot(x,y,z: ary; nrow: integer);
  153. }
  154. {$I PLOT.LIB }
  155.  
  156. begin  { main program }
  157.   clrscr;
  158.   get_data(x,y,nrow);
  159.   linfit(x,y,y_calc,resid,coef,sig,nrow,ncol);
  160.   write_data;
  161.   plot(x,y,y_calc,nrow)
  162. end.
  163.