home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 September / Simtel20_Sept92.cdr / msdos / turbopas / pas_sci.arc / LEAST5.PAS < prev    next >
Pascal/Delphi Source File  |  1985-08-20  |  4KB  |  161 lines
  1. program least5;        { --> 220 }
  2. { Pascal program to perform a linear least-squares fit }
  3. { with Gauss-Jordan routine to Vapor Pressure Equation }
  4. { Separate modules needed:
  5.             GAUSSJ,
  6.             PLOT }
  7.  
  8. const    maxr    = 20;        { data prints }
  9.     maxc    = 5;        { polynomial terms }
  10.  
  11. type    ary    = array[1..maxr] of real;
  12.     arys    = array[1..maxc] of real;
  13.     ary2    = array[1..maxr,1..maxc] of real;
  14.     ary2s    = array[1..maxc,1..maxc] of real;
  15.  
  16. var    x,y,y_calc    : ary;
  17.     resid        : ary;
  18.     coef,sig    : arys;
  19.     nrow,ncol    : integer;
  20.     correl_coef    : real;
  21.     first,done    : boolean;
  22.  
  23. { -> 220 }
  24. procedure get_data(var t    : ary;        { independedt variable }
  25.            var p    : ary;        { dependent variable }
  26.            var nrow    : integer);    { length of vectors }
  27. var    i    : integer;
  28.  
  29. begin
  30.   nrow:=10;
  31.   for i:=1 to nrow do
  32.     t[i]:=(i+6.0)*100.0;
  33.   p[1]:=1.0E-9;        p[2]:=5.598E-8;
  34.   p[3]:=1.234E-6;    p[4]:=1.507E-5;
  35.   p[5]:=1.138E-4;    p[6]:=6.067E-4;
  36.   p[7]:=2.512E-3;    p[8]:=8.337E-3;
  37.   p[9]:=2.371E-2;    p[10]:=5.875E-2;
  38.   for i:=1 to nrow do
  39.     p[i]:=ln(p[i])    { take log data }
  40. end;        { procedure get_data }
  41.  
  42.  
  43. procedure write_data;
  44. { print out the answers }
  45. var    i    : integer;
  46. begin
  47.   if first then first:=false else ClrScr;
  48.   writeln;
  49.   writeln;
  50.   writeln('  I      X       Y      YCALC    RESID');
  51.   for i:=1 to nrow do
  52.     writeln(i:3,x[i]:8:1,y[i]:9:2,y_calc[i]:9:2,resid[i]:9:2);
  53.   writeln; writeln(' Coefficients errors ');
  54.   writeln(coef[1],' ',sig[1],' Constant term');
  55.   for i:=2 to ncol do
  56.     writeln(coef[i],' ',sig[i]);        { other terms }
  57.   writeln;
  58.   writeln('Correlation coefficient is ',correl_coef:8:5)
  59. end;        { write_data }
  60.  
  61. {procedure square(x: ary2;
  62.          y: ary;
  63.          var a: ary2s;
  64.          var g: arys;
  65.      nrow,ncol: integer);}
  66. { matrix multiplication routine }
  67. { a= transpose x times x }
  68. { g= y times x }
  69. {$I SQUARE.LIB }
  70.  
  71. {external procedure gaussj(var b:    ary2s;
  72.                   y:    arys;
  73.               var coef:    arys;
  74.               ncol:        integer;
  75.               var error:    boolean);
  76. }
  77. {$I GAUSSJ.LIB }
  78.  
  79.  
  80. procedure linfit(X,    { independent variable }
  81.          y    : ary;    { dependent variable }
  82.     var y_calc    : ary;    { calculated dep. variable }
  83.     var resid    : ary;    { array of residuals }
  84.     var coef    : arys;    { coefficients }
  85.     var sig        : arys;    { error on coefficients }
  86.         nrow    : integer;    { length of ary }
  87.     var ncol    : integer);    { number of terms }
  88.  
  89. { least-squares fit to nrow sets of x and y pairs of points }
  90. { Seperate  procedure needed:
  91.     SQUARE -> form square coefficient matrix
  92.     GAUSSJ -> Gauus-Jordan elimination }
  93.  
  94. var    xmatr        : ary2;        { data matrix }
  95.     a        : ary2s;    { coefficient matrix }
  96.     g        : arys;        { constant vector }
  97.     error        : boolean;
  98.     i,j,nm        : integer;
  99.     xi,yi,yc,srs,see,
  100.     sum_y,sum_y2    : real;
  101.  
  102. begin        { procedure linfit }
  103.   ncol:=3;    { number of terms }
  104.   for i:=1 to nrow do
  105.     begin    { setup x matrix }
  106.     xi:=x[i];
  107.     xmatr[i,1]:=1.0;    { first column }
  108.     xmatr[i,2]:=1.0/xi;        { second column }
  109.     xmatr[i,3]:=ln(xi)    { third column }
  110.     end;
  111.   square(xmatr,y,a,g,nrow,ncol);
  112.   gaussj(a,g,coef,ncol,error);
  113.   sum_y:=0.0;
  114.   sum_y2:=0.0;
  115.   srs:=0.0;
  116.   for i:=1 to nrow do
  117.     begin
  118.       yi:=y[i];
  119.       yc:=0.0;
  120.       for j:=1 to ncol do
  121.     yc:=yc+coef[j]*xmatr[i,j];
  122.       y_calc[i]:=yc;
  123.       resid[i]:=yc-yi;
  124.       srs:=srs+sqr(resid[i]);
  125.       sum_y:=sum_y+yi;
  126.       sum_y2:=sum_y2+yi*yi
  127.     end;
  128.   correl_coef:=sqrt(1.0-srs/(sum_y2-sqr(sum_y)/nrow));
  129.   if nrow=ncol then nm:=1
  130.   else nm:=nrow-ncol;
  131.   see:=sqrt(srs/nm);
  132.   for i:=1 to ncol do        { errors on solution }
  133.     sig[i]:=see*sqrt(a[i,i])
  134. end;        { LINFIT }
  135.  
  136. {external procedure plot(x,y,z: ary; nrow: integer);
  137. }
  138. {$I C:PLOT.LIB }
  139.  
  140. begin        { main program }
  141.   ClrScr;
  142.   first:=true;
  143.   done:=false;
  144.   writeln;
  145.   get_data(x,y,nrow);
  146.   repeat
  147.     repeat
  148.       write('Order of polynomial fit? ');
  149.       readln(ncol)
  150.     until ncol<5;
  151.     if ncol<1 then done:=true    { quit if ncol<1 }
  152.     else
  153.       begin
  154.     ncol:=ncol+1;        { order is one less }
  155.     linfit(x,y,y_calc,resid,coef,sig,nrow,ncol);
  156.     write_data;
  157.     plot(x,y,y_calc,nrow)
  158.     end    { else }
  159.   until done
  160. end.
  161.