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