home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / TURBOPAS / PAS_ENG.ZIP / GD-LINF1.LIB < prev    next >
Encoding:
Text File  |  1985-07-18  |  2.0 KB  |  77 lines
  1.  
  2. { -> 216 }
  3. procedure get_data(var t : ary;  { independedt variable }
  4.      var cp : ary;  { dependent variable }
  5.      var nrow : integer); { length of vectors }
  6. var i : integer;
  7.  
  8. begin
  9.   nrow:=10;
  10.   for i:=1 to nrow do
  11.     t[i]:=(i+2)*100;
  12.   cp[1]:=7.02; cp[2]:=7.2;
  13.   cp[3]:=7.43; cp[4]:=7.67;
  14.   cp[5]:=7.88; cp[6]:=8.06;
  15.   cp[7]:=8.21; cp[8]:=8.34;
  16.   cp[9]:=8.44; cp[10]:=8.53
  17. end;  { procedure get_data }
  18.  
  19.  
  20. { -> 217 }
  21. procedure linfit(X, { independent variable }
  22.    y : ary; { dependent variable }
  23.  var y_calc : ary; { calculated dep. variable }
  24.  var resid : ary; { array of residuals }
  25.  var coef : arys; { coefficients }
  26.  var sig  : arys; { error on coefficients }
  27.   nrow : integer; { length of ary }
  28.  var ncol : integer); { number of terms }
  29.  
  30. { least-squares fit to nrow sets of x and y pairs of points }
  31. { Seperate  procedure needed:
  32.  SQUARE -> form square coefficient matrix
  33.  GAUSSJ -> Gauus-Jordan elimination }
  34.  
  35. var xmatr  : ary2;  { data matrix }
  36.  a  : ary2s; { coefficient matrix }
  37.  g  : arys;  { constant vector }
  38.  error  : boolean;
  39.  i,j,nm  : integer;
  40.  xi,yi,yc,srs,see,
  41.  sum_y,sum_y2 : real;
  42.  
  43. begin  { procedure linfit }
  44.   ncol:=3; { number of terms }
  45.   for i:=1 to nrow do
  46.     begin { setup x matrix }
  47.  xi:=x[i];
  48.  xmatr[i,1]:=1.0; { first column }
  49.  xmatr[i,2]:=xi;  { second column }
  50.  xmatr[i,3]:=1.0/sqr(xi) { third column }
  51.     end;
  52.   square(xmatr,y,a,g,nrow,ncol);
  53.   gaussj(a,g,coef,ncol,error);
  54.   sum_y:=0.0;
  55.   sum_y2:=0.0;
  56.   srs:=0.0;
  57.   for i:=1 to nrow do
  58.     begin
  59.       yi:=y[i];
  60.       yc:=0.0;
  61.       for j:=1 to ncol do
  62.  yc:=yc+coef[j]*xmatr[i,j];
  63.       y_calc[i]:=yc;
  64.       resid[i]:=yc-yi;
  65.       srs:=srs+sqr(resid[i]);
  66.       sum_y:=sum_y+yi;
  67.       sum_y2:=sum_y2+yi*yi
  68.     end;
  69.   correl_coef:=sqrt(1.0-srs/(sum_y2-sqr(sum_y)/nrow));
  70.   if nrow=ncol then nm:=1
  71.   else nm:=nrow-ncol;
  72.   see:=sqrt(srs/nm);
  73.   for i:=1 to ncol do  { errors on solution }
  74.     sig[i]:=see*sqrt(a[i,i])
  75. end;  { LINFIT }
  76.  
  77.