PROGRAM Least3; { Gauss-Jordan least-squares fit } { From Borland Pascal Programs for Scientists and Engineers } { by Alan R. Miller, Copyright C 1993, SYBEX Inc } USES WinCrt; { Crt for non-windows version } CONST Maxr = 20; { data points } Maxc = 4; { polynomial terms } TYPE FltPt = Real; Ary = ARRAY[1..Maxr] OF FltPt; Arys = ARRAY[1..Maxc] OF FltPt; Ary2 = ARRAY[1..Maxr, 1..Maxc] OF FltPt; Ary2s = ARRAY[1..Maxc, 1..Maxc] OF FltPt; VAR X, Y, Y_Calc, Resid: Ary; Coef, Sig: Arys; Nrow, Ncol: Integer; Correl_Coef: Real; Done: Boolean; {\$I SQUARE} {Listing 3.2} {\$I GAUSSJ} {Listing 4.4} {\$I PLOT} {Listing 5.2} PROCEDURE Get_Data (VAR X : Ary; { independent variable } VAR Y : Ary; { dependent variable } VAR Nrow: Integer); { length of vectors } VAR I: Integer; BEGIN Nrow := 9; FOR I := 1 TO Nrow DO X[I] := I; Y[1]:= 2.07; Y[2]:= 8.6; Y[3]:= 14.42; Y[4]:= 15.80; Y[5]:= 18.92; Y[6]:= 17.96; Y[7]:= 12.98; Y[8]:= 6.45; Y[9]:= 0.27 END; { procedure Get_Data } PROCEDURE Linfit(X, { independent variable } Y: Ary; { dependent variable } VAR Y_Calc: Ary; { calculated dependent variable } VAR Resid : Ary; { array of residuals } VAR Coef : Arys; { coefficients } VAR Sig : Arys; { errors in coefficients } Nrow : Integer; { length of Ary } VAR Ncol : Integer); { number of terms } { least-squares fit to } { Nrow sets of X and Y pairs of points. } { Separate procedures needed: Square and Gaussj } VAR Xmatr: Ary2; { data matrix } A: Ary2s; { coefficient matrix } G: Arys; { constant vector } Error: Boolean; I, J, Nm: Integer; Xi, Yi, Yc, SRS, SEE, Sum_Y, Sum_Y2: Real; BEGIN { procedure Linfit } FOR I := 1 TO Nrow DO BEGIN { Setup X matrix } Xi := X[I]; Xmatr[I, 1] := 1.0; { first column } FOR J := 2 TO Ncol DO { other columns } Xmatr[I, J] := Xmatr[I,J-1] * Xi END; Square(Xmatr, Y, A, G, Nrow, Ncol); Gaussj(A, G, Coef, Ncol, Error); Sum_Y := 0.0; Sum_Y2 := 0.0; SRS := 0.0; FOR I := 1 TO Nrow DO BEGIN Yi := Y[I]; Yc := 0.0; FOR J := 1 TO Ncol DO Yc := Yc + Coef[J] * Xmatr[I, J]; Y_Calc[I] := Yc; Resid[I] := Yc - Yi; SRS := SRS + Sqr(Resid[I]); Sum_Y := Sum_Y + Yi; Sum_Y2 := Sum_Y2 + Yi * Yi END; Correl_Coef := Sqrt(1.0 - SRS /(Sum_Y2 - Sqr(Sum_Y) / Nrow)); IF Nrow = Ncol THEN Nm := 1 ELSE Nm := Nrow - Ncol; SEE := Sqrt(SRS / Nm); FOR I := 1 TO Ncol DO { errors in solution } Sig[I] := SEE * Sqrt(A[I, I]) END; { Linfit } PROCEDURE Write_Data; { print out the answers } VAR I: Integer; BEGIN WriteLn; WriteLn(' I X Y Y CALC RESID'); FOR I := 1 TO Nrow DO WriteLn(I:3, X[I]:8:1, Y[I]:9:2, Y_Calc[I]:9:2, Resid[I]:9:2); WriteLn; WriteLn('coefficients errors'); WriteLn(Coef[1]:12,' ', Sig[1]:12,' Constant term'); FOR I := 2 TO Ncol DO WriteLn (Coef[I]:12,' ', Sig[I]:12); { other terms } WriteLn; WriteLn (' Correlation coefficient is ', Correl_Coef:8:5) END; { Write_Data } BEGIN { main program } Done := False; WriteLn; Get_Data(X, Y, Nrow); REPEAT REPEAT Write(' Order of polynomial fit? '); ReadLn(Ncol) UNTIL Ncol < Maxc; IF Ncol < 1 THEN Done := True { quit if Ncol <1 } ELSE BEGIN Ncol := Ncol +1; { order is one less } Linfit(X, Y, Y_Calc, Resid, Coef, Sig, Nrow, Ncol); Write_Data; WriteLn(' Press Enter for plot'); REPEAT UNTIL KeyPressed; Plot(X, Y, Y_Calc, Nrow) END { ELSE } UNTIL Done; DoneWinCrt { for Windows version only } END.