home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / TURBOPAS / PAS_ENG.ZIP / FITPOL.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-07-19  |  3.7 KB  |  147 lines
  1. program fitpol;  { -> 295 }
  2. { Pascal program to perform a linear least-squares fit }
  3. { to the ratio of 2 polynomials }
  4. { with Gauss-Jordan routine }
  5. { Sperate modules needed:
  6.    GAUSSJ}
  7.  
  8.  
  9. const maxr = 20;  { data prints }
  10.  maxc = 4;  { polynomial terms }
  11.  
  12. type
  13.  ary = array[1..maxr] of real;
  14.  arys = array[1..maxc] of real;
  15.  ary2 = array[1..maxr,1..maxc] of real;
  16.  ary2s = array[1..maxc,1..maxc] of real;
  17.  
  18. var
  19.  x,y,y_calc : ary;
  20.  resid  : ary;
  21.  coef,sig : arys;
  22.  nrow,ncol : integer;
  23.  correl_coef : real;
  24.  
  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.   { clausing factors }
  35.   nrow:=10;
  36.   x[1]:=0.1; y[1]:=0.9524;
  37.   x[2]:=0.2; y[2]:=0.9092;
  38.   x[3]:=0.5; y[3]:=0.8013;
  39.   x[4]:=1.0; y[4]:=0.6720;
  40.   x[5]:=1.2; y[5]:=0.6322;
  41.   x[6]:=1.5; y[6]:=0.5815;
  42.   x[7]:=2.0; y[7]:=0.5142;
  43.   x[8]:=3.0; y[8]:=0.4201;
  44.   x[9]:=4.0; y[9]:=0.3566;
  45.   x[10]:=6.0; y[10]:=0.2755;
  46. end;  { procedure get data }
  47.  
  48. procedure write_data;
  49. { print out the answers }
  50. var i : integer;
  51. begin
  52.   writeln;
  53.   writeln;
  54.   writeln('  I      X      Y      YCALC      RESID');
  55.   for i:=1 to nrow do
  56.     writeln(i:3,x[i]:8:1,y[i]:9:4,y_calc[i]:9:4,resid[i]:9:4);
  57.   writeln; writeln(' Coefficients errors ');
  58.   writeln(coef[1]:8:5,' ',sig[1],' constant term');
  59.   for i:=2 to ncol do
  60.     writeln(coef[i]:8:5,' ',sig[i]);  { other terms }
  61.   writeln;
  62.   writeln('Correlation coefficient is ',correl_coef:8:5)
  63. end;  { write_data }
  64.  
  65. {procedure square(x: ary2;
  66.    y: ary;
  67.       var a: ary2s;
  68.       var g: arys;
  69.   nrow,ncol: integer);}
  70. { matrix multiplication routine }
  71. { a= transpose x times x }
  72. { g= y times x }
  73. {$I SQUARE.LIB }
  74.  
  75. {external procedure gaussj(var b: ary2s;
  76.          y: arys;
  77.      var coef: arys;
  78.      ncol:  integer;
  79.      var error: boolean);
  80. }
  81. {$I GAUSSJ.LIB }
  82.  
  83. procedure linfit(x,  { independant variable }
  84.    y: ary; { dependent variable }
  85.    var y_calc: ary; { calculated dep. variable }
  86.    var resid:  ary; { array of residuals }
  87.    var coef:   arys; { coefficients }
  88.    var sig:    arys; { error on coefficients }
  89.    nrow:       integer; { length of array }
  90.    var ncol:   integer); { number of terms }
  91.  
  92. { least squares fit to nrow sets of x and y pairs of points }
  93. { Seperate procedures needed:
  94.  SQUARE -> form square coefficient matrix
  95.  GAUSSJ -> Gauss-Jordan elimination }
  96.  
  97. var xmatr  : ary2;  { data matrix }
  98.  a  : ary2s; { coefficient matrix }
  99.  g  : arys;  { constant vector }
  100.  error  : boolean;
  101.  i,j,nm  : integer;
  102.  xi,yi,yc,srs,see,
  103.  sum_y,sum_y2 : real;
  104.  
  105. begin  { procedure linfit }
  106.   ncol:=4; { number of terms }
  107.   for i:=1 to nrow do
  108.     begin  { setup matrix }
  109.       xi:=x[i];
  110.       yi:=y[i];
  111.       xmatr[i,1]:=1.0; { first column }
  112.       xmatr[i,2]:=-xi*yi; { second column }
  113.       xmatr[i,3]:=xi; { third column }
  114.       xmatr[i,4]:=-sqr(xi)*yi
  115.     end;
  116.   square(xmatr,y,a,g,nrow,ncol);
  117.   gaussj(a,g,coef,ncol,error);
  118.   sum_y:=0.0;
  119.   sum_y2:=0.0;
  120.   srs:=0.0;
  121.   for i:=1 to nrow do
  122.     begin
  123.       xi:=x[i];
  124.       yi:=y[i];
  125.       yc:=coef[1]+(-coef[2]*yi+coef[3]-coef[4]*xi*yi)*xi;
  126.       y_calc[i]:=yc;
  127.       resid[i]:=yc-yi;
  128.       srs:=srs+sqr(resid[i]);
  129.       sum_y:=sum_y+yi;
  130.       sum_y2:=sum_y2+yi*yi
  131.     end;
  132.   correl_coef:=sqrt(1.0-srs/(sum_y2-sqr(sum_y)/nrow));
  133.   if nrow=ncol then nm:=1
  134.   else nm:=nrow-ncol;
  135.   see:=sqrt(srs/nm);
  136.   for i:=1 to ncol do  { errors on solution }
  137.     sig[i]:=see*sqrt(a[i,i])
  138. end; { linfit }
  139.  
  140.  
  141. begin  { main program }
  142.   clrscr;
  143.   get_data(x,y,nrow);
  144.   linfit(x,y,y_calc,resid,coef,sig,nrow,ncol);
  145.   write_data
  146. end.
  147.