Non-Linear Regression


Language: Turbo Pascal.

Objective: This program performs a non-linear regression between a set of (X,Y) points, according to a function previously defined by the user.

The program has graphical capability, so you can plot dispersion graphics to verify grafically the fit between data and the proposed function.

Notes: The compilation of this program must be done under the Turbo Pascal compiler, version 5.0 or superior. It requires the Turbo Pascal Graphics Toolbox. After compiling, the following files must be included in the same directory of the executable file:


***** Begin of Program Listing *****

Program Merlin;

{--------------------------------------------------------------------------}
{-                                                                        -}
{-          M         E         R         L         I         N           -}
{-                                                                        -}
{-                                                                        -}
{-        Non-Linear Regression Program with Graphical Capability         -}
{-                                                                        -}
{-                                                                        -}
{-    Reference:                                                          -}
{-                                                                        -}
{-      - CACECI, M.S. & CACHERIS, W.P.: Byte, May 1984, 340-62.          -}
{-                                                                        -}
{-                                                                        -}
{-    Unidades Utilizadas:                                                -}
{-    LeastSqr.TPU, GDriver.TPU, GKernel.TPU, GWindow.TPU, GShell.TPU     -}
{-                                                                        -}
{-    Arquivos Utilizados:                                                -}
{-    4x6.fon, 8x8.fon (for IBM color graphics),                          -}
{-    4x9.fon (Hercules mono graphics) and Error.msg.                     -}
{-                                                                        -}
{-                                                                        -}
{-            Antonio Augusto Gorni  ---  March 6, 1990                   -}
{-                                                                        -}
{-                   Lattest Revision: 11/06/1992                         -}
{-                                                                        -}
{--------------------------------------------------------------------------}
{$R+}


uses
  Dos, Crt, Printer, GDriver, GKernel, GWindow, GShell;


type
ColorArray = array[1..4] of byte;
GorniArray = array[1..600] of Real;


const
  IOerr : boolean = false;       { Detecta Erro de E/S }
  LinSup : integer = 1;          { EGA =  1; CGA =  2 }
  LinInf : integer = 25;         { EGA = 25; CGA = 24 }
  PrinterMode = 6;               { EGA =  6; CGA =  1 }
  WindowColor   : ColorArray = (Green, Magenta, LightCyan, White);
  BackColor     : ColorArray = (Black, Black, Black, Black);
  NVPP = 2;     {Numero Total de Variaveis por Ponto}
  MNP = 350;    {Numero Maximo Possivel de Dados}
  ALFA = 1.0;   {Coeficiente de Reflexao, >0}
  BETA = 0.5;   {Coeficiente de Contracao, 0~1}
  GAMMA = 2.0;  {Coeficiente de Expansao, >1}
  ROOT2 = 1.414214;


TYPE
VECTOR = ARRAY[1..8] OF REAL;
DATAROW = ARRAY[1..2] OF REAL;
INDEX = 0..255;
String40 = String[40];


VAR
DONE        :   BOOLEAN;
I,J,K       :   INTEGER;
H,L         :   ARRAY[1..8] OF INDEX;
NP, M, N,
MAXITER,
NITER,NDAT  :   INTEGER;
NEXT,
CENTER,
MEAN, ERROR,
MAXERR,
P, Q,
STEP        :   VECTOR;
SIMP        :   ARRAY[1..8] OF VECTOR;
DATA        :   ARRAY[1..600] OF DATAROW;


var
  Func : integer;
  Narq : string;
  Bandeira, SalvaDados : boolean;
  XData, YData : GorniArray;

GraphArray : PlotArray;


function Potencia(x, y  : real) : real;

begin
  Potencia:=exp(Ln(x)*y);
end; {Potencia}


Function F(X: Vector; D: Datarow) : REAL;

begin
  begin
    F := X[1]/(1+(X[2]-X[3]*D[1]))/X[4];
    Bandeira:=True;
  end
end;


procedure IOCheck;

{---------------------------------------------------}
{- Check for I/O error; print message if needed.   -}
{---------------------------------------------------}

type
  String80 = string[80];

var
  IOcode : integer;

procedure Error(Msg : String80);
begin
  Writeln;
  Write(^G); { Beep! }
  Writeln(Msg);
  Writeln;
end; { procedure Error }

begin { procedure IOCheck }
  IOcode := IOresult;
  IOerr := IOcode <> 0;
  if IOerr then
    case IOcode of
      2   : Error('File not found.');
      3   : Error('Path not found.');
      4   : Error('Too many open files.');
      5   : Error('File access denied.');
      6   : Error('Invalid file handle.');
      12  : Error('Invalid file access code.');
      15  : Error('Invalid drive number.');
      16  : Error('Cannot remove current directory.');
      17  : Error('Cannot rename across drives.');
      100 : Error('Disk read error.');
      101 : Error('Disk write error.');
      102 : Error('File not assigned.');
      103 : Error('File not open.');
      104 : Error('File not open for input.');
      105 : Error('File not open for output.');
      106 : Error('Invalid numeric format.');
    else
      begin
	Writeln;
	Writeln(^G);
	Writeln('Unidentified error message = ', IOcode, '. See manual.');
	Writeln;
      end;
    end; { case }
end; { procedure IOCheck }


function FileExists(Fname : String40) : boolean;
var
  CheckFile : file;
begin
  Assign(CheckFile, Fname);
  {$I-} Reset(CheckFile); {$I+}
   if IOresult = 0 then
   begin
     FileExists := true;
     Close(CheckFile)
   end
 else
   FileExists := false;
end; { function FileExists }


procedure GetGraphData(var X, Y       : GorniArray;
		       var GraphArray : PlotArray;
		       NP : Integer);
var
  Index : Integer;

begin
  for Index := 1 to NP do
  begin
    GraphArray[Index, 1] := X[Index];
    GraphArray[Index, 2] := Y[Index];
  end;
end; { procedure GetGraphData }


procedure SetColors(Fore, Back : Integer);

begin
  if MaxForeGround = 15 then
  begin
    SetForeGroundColor(Fore);
    SetBackGroundColor(Back);
  end;
end;    { Procedure SetColors }


procedure Tela1;

begin
  TextMode(0);
  GoToXY(1,6);
  write('****************************************');
  write('*                                      *');
  write('*      M    E    R    L    I    N      *');
  write('*                                      *');
  write('*                                      *');
  write('*            Developed por:            *');
  write('*                                      *');
  write('*    Antonio Augusto Gorni - EGT/T     *');
  write('*                                      *');
  write('*     Lattest Version: 11/06/1992      *');
  write('*                                      *');
  write('****************************************');
  delay(5000);
  TextMode(2)
end; { procedure Tela1 }


procedure Flash;

begin
  TextBackground(0);
  TextColor(31);
end;


procedure Normal;

begin
  TextBackground(0);
  TextColor(7);
end;


procedure Inverse;

begin
   TextBackground(7);
   TextColor(0);
end;


procedure Centraliza(var Buffer : string);

begin
   GotoXY(round((80-Length(Buffer))/2+1),1);
   Inverse;
   Writeln(Buffer);
   Normal;
end;  {procedure Centraliza}


procedure Tabula(var Buffer, BufImpr : string;
		 Posicao : integer);

var
i, Compr : integer;

begin
  for i := 1 to Posicao - Length(BufImpr) do
    Buffer := ' ' + Buffer;
  BufImpr:=BufImpr+Buffer;
end;   {Tabula}


procedure CabImpr(var BufImpr, Buffer : string;
		  var NroPag : integer);

var
Ano, Mes, Dia, DiaSemana, Hora, Minuto, Segundo, CentSegundo : word;

begin
  Writeln(Lst,BufImpr);
  Writeln(Lst,Buffer);
  GetDate(Ano, Mes, Dia, DiaSemana);
  GetTime(Hora, Minuto, Segundo, CentSegundo);
  Write(Lst,'Page no. ',NroPag,' ***** ');
  Write(Lst,Dia,'/',Mes,'/',Ano,', ',Hora,':');
  If Minuto<10 then Write('0');
  Writeln(Lst,Minuto);
  Writeln(Lst); Writeln(Lst); Writeln(Lst);
end; {CabImpr}


Procedure Eject;

begin
Writeln(Lst,#$C);
end;   {Eject}


procedure Dir(var Caminho : string;
		  Espec : string);

var
Contador, i, j, NroLin : integer;
Flag : boolean;
Buffer : string;
Tipo : word;
ArqDados : SearchRec;
ArqRes : Array [1..255] of string[12];

begin
  Contador := 0;
  Tipo := 0;
  FindFirst(Caminho + Espec, Tipo, ArqDados);
  If DosError = 0 then
  begin
    ArqRes[Contador+1] := ArqDados.Name;
    Contador := Contador + 1;
    Flag := true;
    While Flag do
    begin
      FindNext(ArqDados);
      If DosError = 0 then
      begin
	ArqRes[Contador+1] := ArqDados.Name;
	Contador := Contador + 1;
      end
      else Flag := false;
    end;
  end;
  ClrScr;
  Write('DIRECTORY: '); Inverse; Writeln(Caminho); Normal; Writeln;
  if Contador = 0 then
  begin
    GotoXY(1,7); Flash;
    Writeln('No File Was Found in this Directory!');
    Normal;
    GotoXY(1,14); Write('Press  to Continue... '); Readln(Buffer);
  end
  else
  begin
    For i:=1 to Trunc(Contador/5)+1 do
    begin
      if i mod 20 = 0 then
      begin
	Writeln; Writeln;
	Write('Press  to Continue... '); Readln(Buffer); ClrScr;
	GotoXY(1,4);
	Write('DIRECTORY: '); Inverse; Writeln(Caminho); Normal; Writeln;
      end;
      for j := 1 to 5 do
      begin
	if 5*(i-1)+j <= Contador then
	begin
	  GotoXY((j-1)*15+1,WhereY);
	  write(ArqRes[5*(i-1)+j]);
	end;
      end;
      Writeln;
      end;
    Writeln; Writeln;
    Write('Press  to Continue... '); Readln(Buffer);
  end;
end;     {Dir}


PROCEDURE SUM_OF_RESIDUALS (VAR X : VECTOR);

VAR I : INTEGER;

    BEGIN
      X[N] := 0.0;
      FOR I := 1 TO NP DO
	X[N] := X[N] + SQR(F(X,DATA[I]) - DATA [I,2]);
    END;


Procedure New_Vertex;

var Media : Real;

begin
  Media:=0.0;
  For i:=1 to N do
  begin
    Media:=Media+Error[i];
    Simp[H[N],i]:=Next[i];
  end;
  Media:=Media/N;
  GotoXY(1,15);
  Writeln('# ',Niter:4,' -> ',Media);
end;



PROCEDURE ORDER;

   VAR I,J      : INDEX;

   BEGIN
	FOR J:=1 TO N DO
	   BEGIN
		FOR I:=1 TO N DO
		   BEGIN
			IF SIMP[I,J]SIMP[H[J],J] THEN H[J]:=I
		   END
	   END
   END;


procedure MenuPrincipal(var Funcao : integer);

var
   NroFuncoes : integer;
   Buffer : string;

begin
   Window(1,1,80,25);
   ClrScr;
   Buffer := 'NON-LINEAR REGRESSION';
   Centraliza(Buffer);
   GotoXY(29,6);  Writeln('<1> Data Input');
   GotoXY(29,8);  Writeln('<2> Data Correction');
   GotoXY(29,10); Writeln('<3> Data Suppression');
   GotoXY(29,12); Writeln('<4> Data Output');
   GotoXY(29,14); Writeln('<5> Non-Linear Regression');
   GotoXY(29,16); Writeln('<6> Graphics Plotting');
   GotoXY(29,18); Writeln('<7> End');
   Funcao:=100;
   While (Funcao < 1) or (Funcao > 7) do
   Begin
     GotoXY(29,22); Write('Your Choice? '); Readln(Funcao);
   End;
end;  {procedure MenuPrincipal}


procedure EntradaDados;

var
   Buffer, Xtr, Ytr, Espec, Caminho : string;
   Funcao, Cod, Nr, i : integer;
   InFile : text;

begin
  ClrScr;
  Buffer := 'DATA INPUT';
  Centraliza(Buffer);
  Window(1,4,80,25);
  if SalvaDados then
  begin
    ClrScr;
    write(^G);
    GotoXY(1,9); Flash; Writeln('Data Were Not Saved Yet!'); Normal;
    Buffer := '';
    While ((Buffer <> 'Y') and (Buffer <> 'N')) do
    begin
      GotoXY(1,13); Write('Do You Want to Continue (Y/N)? '); Readln(Buffer);
    end;
    ClrScr;
  end;
  if (Buffer = 'Y') or (not SalvaDados) then
  begin
    GotoXY(32,6);  Writeln('<1> New Data');
    GotoXY(32,10); Writeln('<2> More Data');
    Funcao := 100;
    While ((Funcao < 1) or  (Funcao > 2)) do
    begin
      GotoXY(32,18); Write('Your Choice? '); Readln(Funcao);
    end;
    If Funcao = 1 then
    begin
      ClrScr;
      Buffer:='';
      While (Buffer<>'Y') and (Buffer<>'N') do
      begin
	GotoXY(1,10);
	Write('Proposed Function is Already Programmed (Y/N)? ');
	Readln(Buffer);
      end;
      If Buffer='N' then
      begin
	ClrScr; GotoXY(1,8);
	Writeln('Define the Proposed Function as "Function F" in this Listing;');
	Writeln;
	Writeln('X[i] is a Array with the Constants to be Fitted;');
	Writeln('D[i] is a Array with the Experimental Data;');
	Writeln('ESTADO is a Flag Variable that Signs Problems in the Calculation.');
	Writeln;
	Writeln('Press  to Continue...'); Readln(Buffer); ClrScr;
	Writeln; Writeln;
	Halt;
      end;
      SalvaDados := False;
      NP := 0;
      Narq:='';
    end;
    Buffer := 'INPUT OPTIONS';
    Window(1,1,80,25);
    ClrScr;
    Centraliza(Buffer);
    GotoXY(35,9);  Writeln('<1> Keyboard');
    GotoXY(35,13); Writeln('<2> Disk');
    Funcao := 100;
    While ((Funcao < 1) or (Funcao > 2)) do
    begin
      GotoXY(35,22); Write('Your Choice? '); Readln(Funcao);
    end;
    If Funcao = 1 then
      begin
	Buffer := 'INPUT VIA KEYBOARD';
	Window(1,1,80,25);
	ClrScr;
	Centraliza(Buffer);
	Window(1,4,80,25);
	Writeln('Now Enter the Experimental Data.');
	Writeln('Digit  to Finish.');
	Writeln; Xtr := ''; Ytr := '';
	While ((Xtr <> 'END') and (Ytr <> 'END')) do
	begin
	  NP := NP + 1;
	  Writeln;
	  Writeln('Point #', NP);
	  Write('X? ');
	  Readln(Xtr);
	  If Xtr <> 'END' then
	  begin
	    Write('Y? ');
	    Readln(Ytr);
	    Val(Xtr,XData[NP],Cod);
	    Val(Ytr,YData[NP],Cod);
	  end;
	end;
	NP := NP - 1;
	if NP <> 0 then SalvaDados := true;
      end
      else
      begin
	Buffer := 'INPUT VIA DISK';
	Window(1,1,80,25);
	ClrScr;
	Centraliza(Buffer);
	Window(1,4,80,25);
	Narq := '?';
	While Narq = '?' do
	begin
	  GotoXY(1,6);
	  Write('File Name (? to List Directory)? ');
	  Readln(Narq);
	  If Narq = '?' then
	  begin
	    Writeln; Write('Drive or Path? '); Readln(Caminho);
	    Dir(Caminho,'*.ZDZ'); ClrScr;
	  end;
	end;
	Narq := Narq + '.ZDZ';
	If not FileExists(Narq) then
	begin
	  ClrScr;
	  GotoXY(1,9);
	  Write(^G);
	  Flash; Writeln('File ',Narq,' Does Not Exist!'); Normal;
	  Writeln; GotoXY(1,13);
	  Write('Press  to Continue... ');
	  Readln(Buffer);
	end
	else
	begin
	  If NP <> 0 then SalvaDados := true;
	  GotoXY(1,16); Write('Reading '); Flash; Writeln(Narq);
	  Assign(InFile,Narq);
	  Reset(InFile);
	  {I-} Readln(InFile, Nr); {I+}
	  IOCheck;
	  for i := 1 to Nr do
	  begin
	    {I-} Readln(InFile, XData[NP+i],
			YData[NP+i]); {i+}
	    IOCheck;
	  end;
	  Close(InFile);
	  NP := NP + Nr
	end;
      end;
  end;
end;   {EntradaDados}


Procedure CorrigeDados;

var
Buffer, Buffer1 : string;
Cod, Nponto : integer;

begin
  ClrScr;
  Buffer := 'DATA CORRECTION';
  Centraliza(Buffer);
  Window(1,4,80,25);
  While Buffer <> 'END' do
  begin
    ClrScr;
    GotoXY(1,7); Writeln('Point Index ( to Finish)? ');
    Readln(Buffer);
    If Buffer <> 'END' then
    begin
      Val(Buffer, Nponto, Cod);
      If ((Nponto >= 1) and (Nponto <= NP)) then
      begin
	writeln;
	write('X = ',XData[Nponto],'.   New Value? ');
	Readln(Buffer1);
	If Buffer1 <> '' then Val(Buffer1, XData[Nponto], Cod);
	write('Y = ',YData[Nponto],'.   New Value? ');
	Readln(Buffer1);
	If Buffer1 <> '' then Val(Buffer1, YData[Nponto], Cod);
	SalvaDados := True;
       end;
    end;
  end;
end; {Corrige Dados}


Procedure SuprimeDados;

var
Cod, Nponto, NroSupr, i, j, Troca : integer;
Vsupr : array [1..30] of integer;
Buffer, Buffer1 : string;

begin
  NroSupr := 0;
  ClrScr;
  Buffer := 'SUPRESSAO DE DADOS';
  Centraliza(Buffer);
  Window(1,4,80,25);
  While Buffer <> 'END' do
  begin
    ClrScr;
    GotoXY(1,7); Writeln('Point Index ( to Finish)? ');
    Readln(Buffer);
    If Buffer <> 'END' then
    begin
      Val(Buffer, Nponto, Cod);
      If ((Nponto >= 1) and (Nponto <= NP)) then
      begin
	writeln;
	write('X = ',XData[Nponto],' ***** ',YData[Nponto]);
	Buffer1 := '';
	While ((Buffer1 <> 'Y') and (Buffer1 <> 'N')) do
	begin
	  GotoXY(1,12); Write(^G); Write('Confirm (Y/N)! ');
	  Readln(Buffer1);
	end;
	If Buffer1 = 'Y' then
	begin
	  NroSupr := NroSupr + 1; Vsupr[NroSupr] := Nponto;
	end
	else ClrScr;
      end;
    end;
  end;
  If NroSupr > 0 then
  begin
    GotoXY(1,14); Write(^G); Flash;
    Writeln('Making Suppressions...');
    Normal;
    for i := 1 to NroSupr - 1 do for j := i + 1 to NroSupr do
    If Vsupr[j] > Vsupr[i] then
    begin
      Troca := Vsupr[j]; Vsupr[j] := Vsupr[i]; Vsupr[i] := Troca;
    end;
    for j := 1 to NroSupr do
    begin
      For i := Vsupr[j] to NP - 1 do
      begin
	XData[i] := XData[i+1]; YData[i] := YData[i+1];
      end;
    NP := NP - 1;
    end;
  end;
  SalvaDados := True;
end; {SuprimeDados}


Procedure SaidaDados;

var
Opcao, Contador, i, NroPag : integer;
Buffer, BufImpr, Mensagem, Caminho : string;
OutFile : text;

begin
  ClrScr;
  Buffer := 'DATA OUTPUT';
  Centraliza(Buffer);
  Window(1,4,80,25);
  GotoXY(35,4);  Writeln('<1> Screen');
  GotoXY(35,8);  Writeln('<2> Printer');
  GotoXY(35,12); Writeln('<3> Disk');
  Opcao := 100;
  While ((Opcao < 1) or (Opcao > 3)) do
    begin
    GotoXY(35,18); Write('Your Choice? '); Readln(Opcao);
    end;
  Case Opcao of
    1 : begin
	  Window(1,1,80,25); ClrScr;
	  Buffer := 'DATA OUTPUT';
	  Centraliza(Buffer);
	  Window(1,4,80,25);
	  Contador := 0;
	  For i := 1 to NP do
	  begin
	    Writeln('Point # ',i);
	    Write('X = ',XData[i]);
	    GotoXY(30,WhereY);
	    Writeln('Y = ',YData[i]);
	    Writeln;
	    Contador := Contador + 3;
	    If ((Contador >= 17) and (i < NP)) then
	    begin
	      GotoXY(1,21);
	      Flash; Write('Press  to Continue...'); Normal;
	      Contador := 0;
	      Readln(Buffer);
	      ClrScr; GotoXY(1,1);
	    end;
	  end;
	  GotoXY(1,21);
	  Flash; Write('Press  to Continue.'); Normal;
	  Readln(Buffer);
	  ClrScr;
	end;
    2 : begin
	  Window(1,1,80,25); ClrScr;
	  Buffer := 'PRINTING DATA';
	  Centraliza(Buffer);
	  Window(1,4,80,25); GotoXY(1,8);
	  Writeln('Enter Data Identification Message.');
	  Readln(Mensagem);
	  Writeln; Writeln; Writeln;
	  Write('Prepare Printer; Mark Start of Report.');
	  Writeln; Writeln;
	  Write('Press  to Continue...');
	  NroPag := 1; Readln(Buffer);
	  BufImpr := 'Experimental Data for Non-Linear Regression';
	  CabImpr(BufImpr,Mensagem,NroPag); Contador := 0;
	  For i:=1 to NP do
	  begin
	    Str(i,Buffer);
	    BufImpr := '#' + Buffer;
	    Str(XData[i],Buffer);
	    Tabula(Buffer,BufImpr,10);
	    Str(YData[i],Buffer);
	    Tabula(Buffer,BufImpr,35);
	    Writeln(Lst,BufImpr);
	    Contador := Contador + 1;
	    If ((Contador >= 55) and (i < NP)) then
	    begin
	      NroPag := NroPag + 1; Eject;
	      BufImpr := 'Dados para Ajuste de Curvas';
	      CabImpr(BufImpr,Mensagem,NroPag); Contador := 0;
	    end;
	  end;
	  Eject;
	  ClrScr;
	end;
    3 : begin
	  Window(1,1,80,25); ClrScr;
	  Buffer := 'DATA SAVING IN DISK';
	  Centraliza(Buffer);
	  Window(1,4,80,25);
	  GotoXY(1,6);
	  If Narq<>'' then
	  begin
	    Buffer:='';
	    While ((Buffer<>'Y') and (Buffer<>'N')) do
	    begin
	      GotoXY(1,8); Write('Do You Want to Save in '); Flash; Write(Narq);
	      Normal; Write(' (Y/N)? '); Readln(Buffer); ClrScr;
	    end;
	  end;
	  If ((Narq='') or (Buffer='N')) then
	  begin
	    Narq := '?';
	    While Narq = '?' do
	    begin
	      GotoXY(1,6);
	      Write('File Name (? to List Directory)? ');
	      Readln(Narq);
	      If Narq = '?' then
	      begin
		Writeln; Write('Drive or Path? '); Readln(Caminho);
		Dir(Caminho,'*.ZDZ'); ClrScr;
	      end
			    else Narq:=Narq+'.ZDZ'
	    end;
	  end;
	  GotoXY(1,16); Write('Saving '); Flash; Writeln(Narq);
	  Assign(OutFile,Narq);
	  Rewrite(OutFile);
	  {I-} Writeln(OutFile, NP); {I+}
	  IOCheck;
	  for i := 1 to NP do
	  begin
	    {I-} Writeln(OutFile, XData[i], YData[i]); {i+}
	    IOCheck;
	  end;
	  Close(OutFile); SalvaDados := false;
	end;
      end; {Case}
end; {SaidaDados}


Procedure AjustaCurvas;

var
Contador : index;
Console, Impressora, Buffer, BufImpr, BufAux : string;
a1, b1, s1, s2, y, st, d, p1, d2, be, xe, z, ff, d1, ae, fe, ye,
z1, st1, st2, MediaY : real;
Ano, Mes, Dia, DiaSemana, Hora, Minuto, Segundo, CentSegundo : word;
ep, pn, c : real;
Impressao : array[1..7,1..2] of string[80];

begin
  Console:='Con';
  Impressora:='Lst';
  ClrScr;
  Buffer := 'NON-LINEAR REGRESSION';
  Centraliza(Buffer);
  Window(1,4,80,25);
  GotoXY(1,5);
  Write('Maximum Number of Iteractions? '); Readln(Maxiter); Writeln;
  Write('Number of Constants to Fit? '); Readln(m); n:=m+1;
  Writeln;
  For i:=1 to M do
  begin
    Write('Guess of the Constant ',i,'? '); Readln(Simp[1,i]);
  end;
  Writeln;
  For i:=1 to M do
  begin
    Write('Increment of the Constant ',i,'? '); Readln(Step[i]);
  end;
  Writeln;
  For i:=1 to N do
  If i = n then
  begin
    Write('Precision of the Residuals? '); Readln(Maxerr[i]);
  end
	   else
  begin
    Write('Parameter Precision ',i,'? '); Readln(Maxerr[i]);
  end;
  ClrScr;
  GotoXY(1,8); Writeln('Cogito, Ergo Sum!');
  GotoXY(1,10); Writeln('Maximum Number of Interactions: ',MAXITER);
  for j:=1 to Np do
  begin
    Data[j,1]:=XData[j];
    Data[j,2]:=YData[j]
  end;
  SUM_OF_RESIDUALS(SIMP[1]);
  FOR I:=1 TO M DO
  BEGIN
    P[I]:=STEP[I]*(SQRT(N)+M-1)/(M*ROOT2);
    Q[I]:=STEP[I]*(SQRT(N)-1)/(M*ROOT2)
  END;
  FOR I:=2 TO N DO
  BEGIN
    FOR J:=1 TO M DO SIMP[I,J]:= SIMP[1,J]+Q[J];
    SIMP[I,I-1]:=SIMP[1,I-1]+P[I-1];
    SUM_OF_RESIDUALS(SIMP[I])
  END;
  FOR I:=1 TO N DO
  BEGIN
    L[I]:=1; H[I]:=1
  END;
  ORDER;
  NITER:=0;
  REPEAT
    DONE:=TRUE;
    NITER:=SUCC(NITER);
    FOR I:=1 TO N DO CENTER[I]:=0.0;
    FOR I:=1 TO N DO
      IF I<>H[N] THEN
	FOR J:=1 TO M DO
	  CENTER[J]:=CENTER[J]+SIMP[I,J];
    FOR I:=1 TO N DO
    BEGIN
      CENTER[I]:=CENTER[I]/M;
      NEXT[I]:=(1.0+ALFA)*CENTER[I]-ALFA*SIMP[H[N],I]
    END;
    SUM_OF_RESIDUALS(NEXT);
    IF NEXT[N]<=SIMP[L[N],N] THEN
    BEGIN
      NEW_VERTEX;
      FOR I:=1 TO M DO
	NEXT[I]:=GAMMA*SIMP[H[N],I]+(1.0-GAMMA)*CENTER[I];
      SUM_OF_RESIDUALS(NEXT);
      IF NEXT[N]<=SIMP[L[N],N] THEN NEW_VERTEX
    END
			     ELSE
    IF NEXT[N]<=SIMP[H[N],N] THEN NEW_VERTEX
			     ELSE
    BEGIN
      FOR I:=1 TO M DO
      NEXT[I]:=BETA*SIMP[H[N],I]+(1.0-BETA)*CENTER[I];
      SUM_OF_RESIDUALS(NEXT);
      IF NEXT[N]<=SIMP[H[N],N] THEN NEW_VERTEX
			       ELSE
      FOR I:=1 TO N DO
	BEGIN
	  FOR J:=1 TO M DO SIMP[I,J]:=(SIMP[I,J]+SIMP[L[N],J])*BETA;
	  SUM_OF_RESIDUALS(SIMP[I])
	END
    END;
    ORDER;
    FOR J:=1 TO N DO
    BEGIN
      ERROR[J]:=(SIMP[H[J],J]-SIMP[L[J],J])/SIMP[H[J],J];
      IF DONE THEN
	IF ABS(ERROR[J])>ABS(MAXERR[J]) THEN DONE:=FALSE
    END
  UNTIL (DONE OR (NITER=MAXITER));
  FOR I:=1 TO N DO
  BEGIN
    MEAN[I]:=0.0;
    FOR J:=1 TO N DO MEAN[I]:=MEAN[I]+SIMP[J,I];
    MEAN[I]:=MEAN[I]/N
  END;
  st:=0.; st1:=0.;
  for j:=1 to NP do st1:=st1+YData[j];
  MediaY:=st1/NP;
  st1:=0.;
  for j:= 1 to NP do
  begin
    y:=F(Mean,Data[j]);
    d:=YData[j]-y; st:=st+d*d; d:=YData[j]-MediaY; st1:=st1+d*d;
  end;
  c:=abs(1-(st/(NP-2))/(st1/(NP-1)));
  if c>1 then c:=0.001;
  ep:=sqrt(st/NP); ff:=c*(NP-2)/(1-c);
  p1:=1; d1:=1; d2:=NP-2;
  if ((d1=0) or (d2=0) or (ff=0)) then pn:=p1
				  else
  begin
    if ff < 1 then
    begin
      ae:=d2; be:=d1; fe:=1/ff
    end
	      else
    begin
      ae:=d1; be:=d2; fe:=ff
    end;
    a1:=2/(9*ae); b1:=2/(9*be); xe:=(1-b1)*Potencia(fe,0.3333333)-1+a1;
    ye:=sqrt(b1*Potencia(fe,0.66666667)+a1); z:=abs(xe/ye);
    if be < 4 then z:=z*(1+0.08*Potencia(z,4)/Potencia(be,3));
    z1:=0.115194+z*(0.000344+z*0.019527);
    p1:=0.5/Potencia(1+z*(0.196854+z*z1),4);
    if ff<0 then pn:=p1*100
	    else pn:=(1-p1)*100;
  end;
  if pn=100 then pn:=99.999999; Bandeira:=True;
  ClrScr; GotoXY(1,6); Inverse; Writeln('End of Calculations!'); Normal;
  GotoXY(1,9); Writeln('Maximum Number of Iteractions: ',MaxIter);
  GotoXY(1,11); Writeln('Effective Number of Interactions: ', NIter);
  GotoXY(1,18); Flash; Write('Press  to Continue...'); Normal;
  Readln(Buffer);
  ClrScr; GotoXY(1,6); Writeln('Solutions:'); GotoXY(1,9);
  For i:=1 to M do
  begin
    Writeln('X(',i,') = ',Mean[i]);
    Writeln;
  end;
  GotoXY(1,WhereY+4);
  Flash; Write('Press  to Continue...'); Normal; Readln(Buffer);
  ClrScr; GotoXY(1,6); Writeln('Estimated Fractionary Errors:');
  GotoXY(1,9);
  For i:=1 to M do
  begin
    Writeln('Err(',i,') -> ',Error[i]);
    Writeln;
  end;
  GotoXY(1,WhereY+4);
  Flash; Write('Press  to Continue...'); Normal; Readln(Buffer);
  ClrScr; GotoXY(1,6); Writeln('R^2: ',c); GotoXY(1,9);
  Write('E.P.E.: ',ep); GotoXY(1,12);
  Writeln('% Confidence: ',pn); GotoXY(1,WhereY+4);
  Flash; Write('Press  to Continue...'); Normal; Readln(Buffer);
  ClrScr;
  repeat
    GotoXY(1,6); Write('Do Your Want to Print Results (Y/N)? ');
    Readln(Buffer);
  until (Buffer='Y') or (Buffer='N');
  if Buffer='Y' then
  begin
    ClrScr; GotoXY(1,6); Writeln('Enter Results Identification Message.');
    Readln(Buffer); ClrScr;
    GotoXY(1,6); Write('Prepare Printer and Press !');
    Readln(BufAux);
    Writeln(Lst,'NON-LINEAR REGRESSION');
    Writeln(Lst,'Versao 1.1990, IBM-XT'); Writeln(Lst);
    If Narq<>'' then
    begin
      Write(Lst,'Data from');
      If SalvaDados then Write(Lst,' (partially)');
      Write(Lst,' the File ',Narq,' *** ')
    end
		else
    Write(Lst,'Data Entered via Keyboard *** ');
    GetDate(Ano, Mes, Dia, DiaSemana);
    GetTime(Hora, Minuto, Segundo, CentSegundo);
    Write(Lst,Dia,'/',Mes,'/',Ano,', ',Hora,':');
    If Minuto<10 then Write(Lst,'0');
    Writeln(Lst,Minuto);
    Writeln(Lst,''); Writeln(Lst,Buffer); Writeln(Lst,''); Writeln(Lst,'');
    Writeln(Lst,'Maximum Number of Iteractions: ',MaxIter);
    Writeln(Lst,'Effective Number of Iteractions: ', NIter);
    Writeln(Lst); Writeln(Lst);
    Writeln(Lst,'Solutions:'); Writeln(Lst);
    For i:=1 to M do Writeln(Lst,'X(',i,') = ',Mean[i]);
    Writeln(Lst); Writeln(Lst);
    Writeln(Lst,'Estimated Fractionary Errors:'); Writeln(Lst);
    For i:=1 to M do Writeln(Lst,'Err(',i,') -> ',Error[i]);
    Writeln(Lst); Writeln(Lst);
    Writeln(Lst,'R^2: ',c);
    Writeln(Lst,'E.P.E.: ',ep);
    Writeln(Lst,'% Confidence: ',pn); Writeln(Lst); Writeln(Lst);
  end;
end; {AjustaCurvas}


Procedure PlotaCurvas;

var
k, Opcao, Funcao, Cod, NX, NY, CharCode : integer;
X0, X1, Y0, Y1, XCalc, Incr : real;
Buffer, Console, LabelX, LabelY : string;
Loop : boolean;

begin
   Loop:=true;
   If Bandeira then
   begin
     Opcao:=0;
     X0:=XData[1]; X1:=X0; Y0:=YData[1]; Y1:=Y0;
     for i:=1 to NP-1 do
       for j:=i+1 to NP do
       begin
	if XData[j]X1 then X1:=XData[j];
	if YData[j]Y1 then Y1:=YData[j]
       end;
     NX:=5; NY:=5; LabelX:=''; LabelY:=''; CharCode:=8;
     While Loop do
     begin
       Window(1,1,80,25); ClrScr;
       Buffer := 'GRAPHICS PLOTTING';
       Centraliza(Buffer);
       GotoXY(29,6);  Writeln('<1> Axis Definition');
       GotoXY(29,9); Writeln('<2> Plot Graphic');
       GotoXY(29,12); Writeln('<3> Show Graphic');
       GotoXY(29,15); Writeln('<4> Print Graphic');
       GotoXY(29,18); Writeln('<5> Main Menu');
       Opcao:=0;
       While (Opcao < 1) or (Opcao > 5) do
       begin
	 GotoXY(29,23); Write('Your Choice? '); Readln(Opcao);
       end;
       case Opcao of
	 1 : begin
	       ClrScr; Buffer:='AXIS DEFINITION'; Centraliza(Buffer);
	       Window(1,4,80,25);
	       GotoXY(1,3);
	       Write('Current Minimum X: ',X0,' - New: '); Readln(Buffer);
	       If Buffer<>'' then Val(Buffer,X0,Cod);
	       Write('Current Maximum X: ',X1,' - New: '); Readln(Buffer);
	       If Buffer<>'' then Val(Buffer,X1,Cod);
	       Writeln; Writeln;
	       Write('Current Minimum Y: ',Y0,' - New: '); Readln(Buffer);
	       If Buffer<>'' then Val(Buffer,Y0,Cod);
	       Write('Current Maximum Y: ',Y1,' - New: '); Readln(Buffer);
	       If Buffer<>'' then Val(Buffer,Y1,Cod);
	       Writeln; Writeln;
	       Write('Current Tick Density in the X-Axis: ',NX,
		     ' - New (-9 a +9): ');
	       Readln(Buffer);
	       If Buffer<>'' then Val(Buffer,NX,Cod);
	       Write('Current Tick Density in the Y-Axis: ',NY,
		     ' - New (-9 a +9): ');
	       Readln(Buffer);
	       If Buffer<>'' then Val(Buffer,NY,Cod);
	       If (Abs(NX)>9) or (Abs(NY)>9) then
	       begin
		 Writeln(^G); Inverse;
		 Writeln('Tick Density Illegal!'); Normal
	       end;
	       Writeln; Writeln;
	       Write('Current X-Axis Label: ',LabelX,' - New: ');
	       Readln(Buffer);
	       If Buffer<>'' then LabelX:=Buffer;
	       Write('Current Y-Axis Label: ',LabelY,' - New: ');
	       Readln(Buffer);
	       If Buffer<>'' then LabelY:=Buffer;
	       Writeln; Writeln;
	       Write('Current Point Type: ',CharCode,' - New: ');
	       Readln(Buffer);
	       If Buffer<>'' then Val(Buffer,Charcode,Cod);
	     end;
	 2 : begin
	       RAMScreenGlb := true;
	       Incr:=(X1-X0)/MaxPlotGlb; XCalc:=X0;
	       for k:=1 to MaxPlotGlb do
	       begin
		 GraphArray[k,1]:=XCalc;
		 Data[Np+1,1]:=XCalc;
		 GraphArray[k,2]:=F(Mean, Data[Np+1]);
		 XCalc:=XCalc+Incr;
		 If (GraphArray[k,1]X1) or
		    (GraphArray[k,2]Y1) or
		    not Bandeira then
		    begin
		    GraphArray[k,1]:=X0;
		    GraphArray[k,2]:=Y0;
		    end;
	       end;
	       EnterGraphic;
	       DefineWindow(1, 0, 0, XMaxGlb, YMaxGlb);
	       DefineWorld(1,X0,Y0,X1,Y1);
	       SelectWorld(1); SelectWindow(1); SetBackground(0);
	       ClearScreen;
	       GotoXY(1,LinSup); Writeln(LabelY);
	       GotoXY(80-Length(LabelX),LinInf); Writeln(LabelX);
	       DrawBorder;
	       DrawAxis(NX, NY, 0, 15, 0, 15, -1, -1, true);
	       DrawPolygon(GraphArray, 1, MaxPlotGlb, -9, 1, 0);
	       ResetAxis;
	       GetGraphData(XData, YData, GraphArray, NP);
	       DrawPolygon(GraphArray, 1, NP, -CharCode, 2, 0);
	       Readln(Buffer);
	       SwapScreen;
	       LeaveGraphic;
	     end;
	 3 : begin
	       EnterGraphic;
	       SwapScreen;
	       Readln(Buffer);
	       SwapScreen;
	       LeaveGraphic;
	     end;
	 4 : begin
	       Window(1,1,80,25); ClrScr;
	       Buffer:='GRAPHIC PRINTING';
	       Centraliza(Buffer);
	       GotoXY(1,12);
	       Writeln('Prepare Printer and Press  to Continue!');
	       Readln(Buffer);
	       EnterGraphic;
	       SwapScreen;
	       HardCopy(false,PrinterMode);
	       SwapScreen;
	       LeaveGraphic;
	       Writeln(Lst,Chr(27),'@');
	     end;
	 5 : Loop:=false;
       end {case}
    end
  end;
end;  {PlotaCurvas}


procedure Fim;

var
Buffer : string;

begin
  If SalvaDados then
  begin
    ClrScr;
    Buffer := 'END OF PROGRAM RUN';
    Centraliza(Buffer);
    Window(1,4,80,25);
    GotoXY(1,9); Flash;
    Writeln(^G'Data Not Saved Yet!');
    Normal;
    Buffer := '';
    While ((Buffer<>'Y') and (Buffer<>'N')) do
    begin
      GotoXY(1,12); Write('Quit Program (Y/N)? '); Readln(Buffer);
    end;
    if Buffer = 'Y' then
    begin
      ClrScr;
      Halt;
    end;
  end
		     else
  begin
    ClrScr;
    Halt
  end;
end;   {Fim}


begin { program Merlin }
  InitGraphic;
  SetColors(WindowColor[1], BackColor[1]);
  LeaveGraphic;
  Tela1;
  SalvaDados := false;
  Bandeira := false;
  While true do
  begin
    MenuPrincipal(Func);
    case Func of
    1 : EntradaDados;
    2 : CorrigeDados;
    3 : SuprimeDados;
    4 : SaidaDados;
    5 : AjustaCurvas;
    6 : PlotaCurvas;
    7 : Fim;
    end; {case}
  end;
end. { program Merlin }

***** End of Program Listing ******

Return to the Software Menu.

Last Update: 28 June 1996
© Antonio Augusto Gorni