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