Web pdp-10.trailing-edge.com

Trailing-Edge - PDP-10 Archives - decuslib10-08 - 43,50504/mulreg.sim
There are 2 other files named mulreg.sim in the archive. Click here to see a list.
```BEGIN	COMMENT     Title : Multiple Linear Regression Analysis,
Author : Marten van Gelderen,
Version : 5H(246), Simula 67, DEC-10, Date 800216,
Copyright : Foundation Mathematical Centre, Amsterdam;

COMMENT ************************************* external procedures ********************************;

EXTERNAL REAL PROCEDURE Maxreal;
EXTERNAL BOOLEAN PROCEDURE Dotypeout;
EXTERNAL REF (Infile) PROCEDURE Findinfile;
EXTERNAL REF (Printfile) PROCEDURE Findprintfile;
EXTERNAL REF (Directfile) PROCEDURE Finddirectfile;
EXTERNAL LONG REAL PROCEDURE Scanreal, Lmin, Lmax;
EXTERNAL INTEGER PROCEDURE Maxint, Scanint, Imin, Imax;
EXTERNAL TEXT PROCEDURE Upcase, Compress, Rest, Scanto, Filename;

INTEGER scalelimit, optionlimit; scalelimit:= 3; optionlimit:= 10; Lowten ('#'); Simset

BEGIN	INTEGER ARRAY hashprime [1 : scalelimit]; BOOLEAN ARRAY optionwanted [1 : optionlimit];

COMMENT ************************************* error messages *************************************;

BOOLEAN erroneous;

PROCEDURE overflow (string); NAME string; TEXT string;
INSPECT printstream DO
BEGIN	Eject (Line + 5); Outtext (string); Outtext (" overflow");
IF tablescale LT scalelimit THEN
BEGIN	Outtext ("    (message)"); Outimage;
Eject (Line + 5); Outtext ("Rerun of job :"); Outint (jobnumber, 3); Outimage;
Eject (1); tablescale:= tablescale + 1; stacklimit:= stacklimit * 2; GOTO rerun;
END; Outimage; erroneous:= TRUE; GOTO endjob;
END *** overflow ***;

PROCEDURE errortext (errnum); INTEGER errnum;
INSPECT printstream DO
BEGIN	Eject (Line + 1); Outtext ("Error : "); INSPECT errorstream DO
BEGIN	Locate (errnum - errnum // 10 + 11); Inimage;
printstream . Outtext (Image . Strip); printstream . Outimage;
END OTHERWISE Outint (errnum, 2); erroneous:= TRUE;
END *** errortext ***;

PROCEDURE error (errnum); INTEGER errnum;
INSPECT printstream DO
BEGIN	REF (atomcell) atom; errortext (errnum);
FOR atom:- lastname, lastnumber, lastsymbol DO
BEGIN Setpos (Pos + 13); printitem (atom); END; Outimage;
END *** error ***;

PROCEDURE error1 (errnum, term); INTEGER errnum, term;
INSPECT printstream DO
BEGIN	errortext (errnum); Outint (term, 14); Outimage; END;

PROCEDURE error2 (errnum, wrong, case); INTEGER errnum; REAL wrong, case;
INSPECT printstream DO
BEGIN	errortext (errnum); Outfix (wrong, 5, 20); Outfix (case, 5, 20); Outimage; GOTO endjob; END;

COMMENT ************************************* file handling **************************************;

INTEGER imagelength, errorlength;

REF (Infile) inputstream;	REF (Printfile) printstream;
REF (Infile) datastream;	REF (Printfile) outputstream;	   REF (Directfile) errorstream;

TEXT PROCEDURE enterspec (stream); NAME stream; TEXT stream;
BEGIN Outtext (stream); Breakoutimage; Inimage; enterspec:- Sysin . Image . Strip; END;

REF (Infile) PROCEDURE getinpfile (filespec, stream); NAME stream; TEXT filespec, stream;
IF filespec == NOTEXT THEN getinpfile:- NONE ELSE INSPECT Findinfile (filespec) DO
getinpfile:- THIS Infile OTHERWISE getinpfile:- getinpfile (enterspec (stream), stream);

REF (Printfile) PROCEDURE getprifile (filespec, stream); NAME stream; TEXT filespec, stream;
IF filespec == NOTEXT THEN getprifile:- NONE ELSE INSPECT Findprintfile (filespec) DO
getprifile:- THIS Printfile OTHERWISE getprifile:- getprifile (enterspec (stream), stream);

PROCEDURE enterfilespecs;
BEGIN	Eject (Line + 1); Outtext ("Enter file specifications"); Outimage;
inputstream  :- getinpfile (enterspec ("Inputstream  : "), "Inputstream  ? ");
printstream  :- getprifile (enterspec ("Printstream  : "), "Printstream  ? ");
datastream   :- getinpfile (enterspec ("Datastream   : "), "Datastream   ? ");
outputstream :- getprifile (enterspec ("Outputstream : "), "Outputstream ? ");
errorstream  :- Finddirectfile ("Hlp: Mulerr.hlp", FALSE);
Eject (Line + 1); Outtext ("Executing ..."); Outimage; Outimage;
END *** enterfilespecs ***;

IF filespec . More AND filespec . Strip NE Filename (inputstream) . Strip THEN
BEGIN	IF inputstream =/= Sysin THEN inputstream . Close;
inputstream:- getinpfile (filespec . Strip, "Inputstream  ? "); runable:= FALSE;
INSPECT inputstream DO Open (Blanks (imagelength)) OTHERWISE inputstream:- Sysin;

PROCEDURE readprispec (filespec); TEXT filespec; IF filespec . More THEN
BEGIN	IF printstream =/= Sysout THEN printstream . Close;
printstream:- getprifile (filespec . Strip, "Printstream  ? ");
INSPECT printstream DO Open (Blanks (imagelength)) OTHERWISE printstream:- Sysout;

PROCEDURE readdatspec (filespec); TEXT filespec; IF filespec . More THEN
BEGIN	INSPECT datastream DO Close; datastream:- getinpfile (filespec . Strip, "Datastream   ? ");
INSPECT datastream DO Open (Blanks (imagelength)); datarecord:= 0; datafileread:= FALSE;

PROCEDURE readoutspec (filespec); TEXT filespec; IF filespec . More THEN
BEGIN	INSPECT outputstream DO Close; outputstream:- getprifile (filespec . Strip, "Outputstream ? ");
INSPECT outputstream DO Open (Blanks (imagelength));

BEGIN	runable:= TRUE; IF lastchar EQ '/' AND lastline . More THEN
BEGIN	TEXT scanline;
lastline:- Copy (Rest (lastline)); scanline:- Scanto (lastline, '=');
IF NOT runable THEN readline (inputstream);
END;

PROCEDURE openfiles;
BEGIN	INSPECT inputstream  DO Open (Blanks (imagelength)) OTHERWISE inputstream:- Sysin;
INSPECT printstream  DO Open (Blanks (imagelength)) OTHERWISE printstream:- Sysout;
INSPECT datastream   DO Open (Blanks (imagelength)); Sysout . Linesperpage (-1);
INSPECT outputstream DO Open (Blanks (imagelength)); Sysout . Image:- Blanks (imagelength);
INSPECT errorstream  DO Open (Blanks (errorlength));
END *** openfiles ***;

PROCEDURE closefiles;
BEGIN	INSPECT inputstream DO Close;	INSPECT printstream  DO Close;
INSPECT datastream  DO Close;	INSPECT outputstream DO Close;	INSPECT errorstream DO Close;
END *** closefiles ***;

COMMENT ************************************* text & data storage ********************************;

INTEGER bufferlength, resetindex;

Link CLASS textbuffer; BEGIN TEXT line; END;

CLASS textstorage;
BEGIN	INTEGER size;
REF (textbuffer) pointer; REF (Head) base;

PROCEDURE reset;
BEGIN pointer:- base . First; endtext:= FALSE; endline:= TRUE; END;

PROCEDURE store (line); TEXT line;
BEGIN	pointer:- NEW textbuffer;
pointer . Into (base); pointer . line:- Copy (line);
pointer . line . Setpos (line . Pos - (IF endline THEN 0 ELSE 1));
size:= size + (IF endline THEN 0 ELSE 1);
END *** store ***;

TEXT PROCEDURE nextline;
BEGIN	nextline:- lastline:- pointer . line;
pointer:- pointer . Suc; endtext:= pointer == NONE;
END *** nextline ***;

size:= 0; pointer:- NONE; base:- NEW Head;
END *** textstorage ***;

BEGIN LONG REAL ARRAY buffer [1 : bufferlength]; END;

CLASS datastorage;
BEGIN	INTEGER index, size;
REF (databuffer) pointer; REF (Head) base;

PROCEDURE store (number); LONG REAL number;
BEGIN	IF index GT bufferlength THEN
BEGIN	pointer:- NEW databuffer;
pointer . Into (base); index:= 1;
END; pointer . buffer [index]:= number;
index:= index + 1; size:= size + 1;
END *** store ***;

PROCEDURE bufferror; INSPECT printstream DO
BEGIN errortext (3); Outimage; GOTO endjob; END;

size:= 0; index:= bufferlength + 1;
END *** datastorage ***;

datastorage CLASS inputstorage;
BEGIN	PROCEDURE reset;
BEGIN pointer:- base . First; index:= 0; END;

LONG REAL PROCEDURE nextnumber;
BEGIN	IF index LT bufferlength THEN index:= index + 1
ELSE BEGIN pointer:- pointer . Suc; index:= 1; END;
IF pointer == NONE THEN bufferror;
nextnumber:= pointer . buffer [index];
END *** nextnumber ***;

PROCEDURE skipnumbers (count); INTEGER count;
BEGIN	index:= index + count; WHILE index GE bufferlength DO
BEGIN	index:= index - bufferlength; pointer:- pointer . Suc;
IF pointer == NONE THEN bufferror;
END;
END *** skipnumbers ***;
END *** inputstorage ***;

datastorage CLASS rightstorage;

LONG REAL PROCEDURE lastnumber;
BEGIN	IF index GT 1 THEN index:= index - 1 ELSE
BEGIN pointer:- pointer . Pred; index:= bufferlength; END;
IF pointer == NONE THEN bufferror;
lastnumber:= pointer . buffer [index];
END *** lastnumber ***;
END *** rightstorage ***;

rightstorage CLASS leftstorage;
BEGIN	INTEGER repeats;

PROCEDURE fixrepeats;
BEGIN	repeats:= size - repeats;
store (repeats); repeats:= size;
END *** fixrepeats ***;

repeats:= 0;
END *** leftstorage ***;

rightstorage CLASS designstorage;
BEGIN	LONG REAL sum, squares, min, max;

PROCEDURE reset;
BEGIN pointer:- base . Last; index:= resetindex; END;

sum:= squares:= 0; min:= Maxreal; max:= - Maxreal;
END *** designstorage ***;

COMMENT ************************************* atom cells *****************************************;

CLASS atomcell;
BEGIN TEXT item, entry; INTEGER index, type; END;

atomcell CLASS identifier;;	atomcell CLASS number;;
atomcell CLASS functionname;;	atomcell CLASS optionname;;

REF (atomcell)	lastatom, lastname, lastnumber, lastsymbol,
addsign, subsign, mulsign, divsign, idisign, ttpsign, withsign,
equal, comma, lparent, rparent, lbracket, rbracket, smaller, greater, semicolon;

COMMENT ************************************* item & atom hashing ********************************;

INTEGER maxrank; BOOLEAN declared; REF (hashing) atomhash;

CLASS hashing (hashprime); INTEGER hashprime;
BEGIN	INTEGER twinprime, key, maxkey, vacant, hashlimit, index, shift; TEXT entry;

REF (atomcell) ARRAY hashtable [1 : hashprime];

PROCEDURE hasher;
BEGIN	key:= maxkey - maxrank; WHILE entry . More DO
key:= Mod (entry . Pos + key * Rank (entry . Getchar), maxkey);
END *** hasher ***;

BOOLEAN PROCEDURE present;
BEGIN	entry:- Upcase (Compress (Copy (lastitem), ' ')); hasher;
index:= Mod (key, hashprime) + 1; shift:= Mod (key, twinprime) + 1;
declared:= FALSE; WHILE NOT declared AND hashtable [index] =/= NONE DO
IF hashtable [index] . entry EQ entry THEN
declared:= TRUE ELSE index:= Mod (index + shift, hashprime) + 1;
present:= declared;
END *** present ***;

REF (atomcell) PROCEDURE insert (newcell); REF (atomcell) newcell;
BEGIN	IF vacant LT hashlimit THEN overflow ("Hash table"); vacant:= vacant - 1;
newcell . item:- lastitem; newcell . entry:- entry; newcell . index:= index;
insert:- hashtable [index]:- newcell; newcell . entry . Setpos (1);
END *** insert ***;

REF (atomcell) PROCEDURE retrieve; retrieve:- hashtable [index];

REF (atomcell) PROCEDURE reserve (newcell, item, type);
VALUE item; REF (atomcell) newcell; TEXT item; INTEGER type;
BEGIN	REF (atomcell) atomtype; lastitem:- item;
atomtype:- reserve:- IF present THEN retrieve ELSE insert (newcell);
atomtype . type:= type;
END *** reserve ***;

twinprime:= hashprime - 2; maxkey:= Maxint // maxrank;
vacant:= hashprime; hashlimit:= Entier (0.2 * hashprime);
FOR index:= 1 STEP 1 UNTIL hashprime DO hashtable [index]:- NONE;

addsign   :- reserve (NEW atomcell, "+", 0);
subsign   :- reserve (NEW atomcell, "-", 0);
mulsign   :- reserve (NEW atomcell, "*", 0);
divsign   :- reserve (NEW atomcell, "/", 0);
idisign   :- reserve (NEW atomcell, ":", 0);
ttpsign   :- reserve (NEW atomcell, "^", 0);
withsign  :- reserve (NEW atomcell, "&", 0);
equal     :- reserve (NEW atomcell, "=", 0);
comma     :- reserve (NEW atomcell, ",", 0);
lparent   :- reserve (NEW atomcell, "(", 0);
rparent   :- reserve (NEW atomcell, ")", 0);
lbracket  :- reserve (NEW atomcell, "[", 0);
rbracket  :- reserve (NEW atomcell, "]", 0);
smaller   :- reserve (NEW atomcell, "<", 0);
greater   :- reserve (NEW atomcell, ">", 0);
semicolon :- reserve (NEW atomcell, ";", 0);
reserve (NEW functionname, "Abs", 1) . index:= 1;
reserve (NEW functionname, "Sign", 2) . index:= 1;
reserve (NEW functionname, "Sqrt", 3) . index:= 1;
reserve (NEW functionname, "Sin", 4) . index:= 1;
reserve (NEW functionname, "Cos", 5) . index:= 1;
reserve (NEW functionname, "Tan", 6) . index:= 1;
reserve (NEW functionname, "Ln", 7) . index:= 1;
reserve (NEW functionname, "Log", 8) . index:= 1;
reserve (NEW functionname, "Exp", 9) . index:= 1;
reserve (NEW functionname, "Entier", 10) . index:= 1;
reserve (NEW functionname, "Round", 11) . index:= 1;
reserve (NEW functionname, "Mod", 12) . index:= 2;
reserve (NEW functionname, "Min", 13) . index:= 2;
reserve (NEW functionname, "Max", 14) . index:= 2;
reserve (NEW functionname, "Arcsin", 15) . index:= 1;
reserve (NEW functionname, "Arccos", 16) . index:= 1;
reserve (NEW functionname, "Arctan", 17) . index:= 1;
reserve (NEW functionname, "Sinh", 18) . index:= 1;
reserve (NEW functionname, "Cosh", 19) . index:= 1;
reserve (NEW functionname, "Tanh", 20) . index:= 1;
reserve (NEW functionname, "Indicator", 21) . index:= 3;
reserve (NEW optionname, "Transformed data matrix", 1);
reserve (NEW optionname, "Correlation matrix", 2);
reserve (NEW optionname, "Residual analysis", 3);
reserve (NEW optionname, "No regression analysis", 4);
reserve (NEW optionname, "Process submodels", 5);
reserve (NEW optionname, "Print input data", 6);
reserve (NEW optionname, "No input data rewind", 7);
reserve (NEW optionname, "Save original model", 8);
reserve (NEW optionname, "Test reduced model", 9);
reserve (NEW optionname, "Missing values", 10);
END *** hashing ***;

COMMENT ************************************* text & data handling *******************************;

INTEGER datarecord; TEXT lastline, lastitem; CHARACTER lastchar;
REF (textstorage) currenttext, modeltext, inputtext, optiontext; REF (inputstorage) inputdata;

NAME sourcetext; REF (textstorage) sourcetext;
BEGIN	sourcetext:- NEW textstorage; INSPECT sourcetext DO
BEGIN	WHILE NOT endtext DO
BEGIN store (lastline); readline (inputstream); END;
IF size EQ 0 THEN sourcetext:- NONE;
END;

PROCEDURE readline (stream); REF (Infile) stream;
INSPECT stream DO
BEGIN	IF stream == Sysin THEN BEGIN Outchar ('*'); Breakoutimage; END; Inimage;
lastline:- Image . Strip; WHILE nextchar LT ' ' AND NOT endline DO;
endtext:= lastchar EQ '"' OR Endfile;

PROCEDURE readdata (stream); REF (Infile) stream;
BEGIN	INTEGER start; LONG REAL nextreal; erroneous:= FALSE;
inputdata:- NEW inputstorage; INSPECT inputdata DO
BEGIN	WHILE NOT endtext DO
BEGIN	WHILE NOT endline DO
BEGIN	IF IF integral (lastchar) THEN TRUE ELSE numerical (lastchar) THEN
BEGIN	start:= lastline . Pos - 1; lastline . Setpos (start);
nextreal:= Scanreal (lastline); IF start EQ lastline . Pos
THEN error1 (4, size + 1) ELSE store (nextreal);
lastline . Setpos (lastline . Pos + 1);
END; WHILE nextchar LT ' ' AND NOT endline DO;
IF dataread THEN BEGIN store (size); reset; END;
END; INSPECT datastream DO IF Endfile AND datafileread THEN
BEGIN Close; datastream:- NONE; END ELSE datafileread:= stream == datastream;

INSPECT datastream DO
WHILE (endtext OR endline) AND NOT Endfile DO
BEGIN	IF endtext THEN
BEGIN	datarecord:= datarecord + 1; endtext:= FALSE;
END; IF endline THEN readline (datastream);
END; readdata (datastream); datarecord:= datarecord + 1;

BEGIN	readkeyword:- Upcase (Copy (nextitem (2)));
WHILE nextchar NE '"' AND NOT endline DO;
WHILE nextchar LT ' ' AND NOT endline DO;

runable:= FALSE; WHILE NOT runable DO
BEGIN	echotext; keyword:- readkeyword; endtext:= FALSE;
IF keyword EQ "RU" THEN readfilespecs ELSE
IF keyword EQ "MO" THEN readtext (modeltext) ELSE
IF keyword EQ "IN" THEN readtext (inputtext) ELSE
IF keyword EQ "OP" THEN readtext (optiontext) ELSE
IF keyword EQ "DA" THEN readdata (inputstream) ELSE
IF keyword EQ "HE" THEN printhelp ELSE
IF keyword EQ "EX" THEN GOTO exit;
END; jobnumber:= jobnumber + 1;

PROCEDURE printtext;
INSPECT printstream DO INSPECT currenttext DO
BEGIN	reset; Eject (Line + 5);
WHILE NOT endtext DO BEGIN Outtext (nextline); Outimage; END; reset;
END *** printtext ***;

PROCEDURE printdata;
INSPECT printstream DO INSPECT inputdata DO
BEGIN	INTEGER i; Eject (Line + 5); Outtext ("""Data"""); Outimage; Eject (Line + 1);
reset; FOR i:= size - 1 STEP - 1 UNTIL 1 DO Outfix (nextnumber, 3, 12); Outimage;
END *** printdata ***;

PROCEDURE printitem (atom); REF (atomcell) atom;
INSPECT printstream DO INSPECT atom DO Outtext (item . Sub (1, Imin (item . Length, 8)));

PROCEDURE printhelp;
BEGIN	TEXT filespec; filespec:- Copy ("Hlp: Mulreg.hlp");
IF lastchar EQ '/' THEN filespec . Sub (9, 3):= nextitem (3);
INSPECT printstream DO INSPECT Findinfile (filespec) DO
BEGIN	Open (Blanks (imagelength)); Inimage; WHILE NOT Endfile DO
BEGIN Outtext (Image . Strip); Outimage; Inimage; END; Close;
IF printstream =/= Sysout THEN Eject (1); Dotypeout (Sysout);
INSPECT Sysout DO BEGIN Outtext ("End of helptext"); Outimage; END;
END OTHERWISE INSPECT Sysout DO
BEGIN Outtext ("No help information available"); Outimage; END;
END *** printhelp ***;

INTEGER PROCEDURE field (w); INTEGER w;
field:= IF w EQ 0 THEN 1 ELSE Entier (Ln (Abs (w)) / 2.3) + (IF w LT 0 THEN 2 ELSE 1);

PROCEDURE echotext;
INSPECT printstream DO
WHILE NOT endtext OR inputstream . Endfile DO
BEGIN	WHILE NOT endtext DO
BEGIN	IF inputstream == Sysin AND printstream == Sysout THEN
BEGIN	IF lastline =/= NOTEXT THEN Outchar ('?')
ELSE Outtext ("For help type: ""Help""");
END; Outtext (lastline); Outimage; readline (inputstream);
END; INSPECT inputstream DO IF Endfile THEN
BEGIN Close; inputstream:- Sysin; readline (inputstream); END;
END *** echotext ***;

CHARACTER PROCEDURE nextchar;
BEGIN	endline:= NOT lastline . More; lastchar:= ' ';
WHILE lastchar EQ ' ' AND lastline . More DO lastchar:= lastline . Getchar;
nextchar:= lastchar;
END *** nextchar ***;

BOOLEAN PROCEDURE integral (char); CHARACTER char;
integral:= Digit (char) OR char EQ '+' OR char EQ '-';

BOOLEAN PROCEDURE numerical (char); CHARACTER char;
numerical:= char EQ '.' OR char EQ '#';

TEXT PROCEDURE nextitem (length); INTEGER length;
BEGIN	WHILE nextchar LT ' ' AND NOT endline DO; nextitem:-
lastline . Sub (lastline . Pos - 1, Imin (length, lastline . Length - lastline . Pos + 2));
END *** nextitem ***;

REF (atomcell) PROCEDURE nextatom;
INSPECT currenttext DO INSPECT atomhash DO
BEGIN	WHILE endline AND NOT endtext DO BEGIN nextline; nextchar; END;
IF endline AND endtext THEN BEGIN nextatom:- lastatom; END ELSE
BEGIN	INTEGER start, width; start:= lastline . Pos - 1;
IF Letter (lastchar) THEN
BEGIN	WHILE Letter (nextchar) OR Digit (lastchar) DO;
width:= lastline . Pos - start - 1;
IF endline THEN width:= width + 1;
lastitem:- lastline . Sub (start, width) . Strip;
nextatom:- lastatom:- lastname:-
IF present THEN retrieve ELSE insert (NEW identifier);
END ELSE
IF IF Digit (lastchar) THEN TRUE ELSE numerical (lastchar) THEN
BEGIN	IF Digit (lastchar) THEN WHILE Digit (nextchar) DO;
IF lastchar EQ '.' THEN
BEGIN	IF NOT Digit (nextchar) THEN error (5)
ELSE WHILE Digit (nextchar) DO;
END;
IF lastchar EQ '#' THEN
BEGIN	IF NOT integral (nextchar) THEN error (6)
ELSE WHILE Digit (nextchar) DO;
END;
width:= lastline . Pos - start - 1;
IF endline THEN width:= width + 1;
lastitem:- lastline . Sub (start, width) . Strip;
nextatom:- lastatom:- lastnumber:-
IF present THEN retrieve ELSE insert (NEW number);
END ELSE
BEGIN	nextchar; lastitem:- lastline . Sub (start, 1);
nextatom:- lastatom:- lastsymbol:-
IF present THEN retrieve ELSE insert (NEW atomcell);
END; WHILE lastchar LT ' ' AND NOT endline DO nextchar;
END;
END *** nextatom ***;

COMMENT ************************************* system variables ***********************************;

INTEGER stk, add, sub, mul, div, idi, ttp, neg, eql, rep, dep, res, tst, ini,
frs, fls, smc, tcv, tav, rav, sav, lav, rtn, jmp, fix, skp, fun, savedf, jobnumber,
termcounter, leftcounter, rightcounter, designcounter, weightcounter, programcounter,
tablescale, programscale, stacklimit, tablelimit, programspace, stackspace, stacklength;
LONG REAL precision, Ln10, maxobs, maxexp, savess, savems; BOOLEAN runable, savemodel, savefit;

PROCEDURE initializesystem;
BEGIN	Outtext ("Multiple Linear Regression Analysis"); Outimage;
hashprime [1]:= 271; hashprime [2]:= 523; hashprime [3]:= 1033;
jobnumber:= datarecord:= 0; tablescale:= 1; programscale:= 1; stacklimit:= 32;
bufferlength:= 128; imagelength:= 132; errorlength:= 75; maxrank:= 131;
precision:= &&-18; Ln10:= 2.30258509299404572&&0; maxobs:= &&19; maxexp:= 88;
stk:= tcv:= 1; add:= tav:= 2; sub:= rav:= 3; mul:= sav:= 4; div:= lav:= 5;
idi:= rtn:= 6; ttp:= jmp:= 7; neg:= fix:= 8; eql:= skp:= 9; rep:= fun:= 10;
dep:= 11; res:= 12; tst:= 13; ini:= 14; frs:= 15; fls:= 16; smc:= 17;
modeltext:- inputtext:- optiontext:- NONE; inputdata:- NONE;
END *** initializesystem ***;

PROCEDURE initializerun;
BEGIN	lastatom:- lastname:- lastnumber:- lastsymbol:- NONE;
programcounter:= 1; leftcounter:= rightcounter:= 0; erroneous:= FALSE;
tablelimit:= hashprime [tablescale]; programspace:= tablelimit * programscale;
atomhash:- NEW hashing (tablelimit); stacklength:= stackspace:= tablelimit + stacklimit;
END *** initializerun ***;

initializesystem; enterfilespecs; openfiles;

job : readjob;  rerun : initializerun;

BEGIN	INTEGER ARRAY program [1 : programspace]; LONG REAL ARRAY stack [1 : stackspace];

IF programcounter LE programspace THEN
BEGIN	program [programcounter]:= address * 100 + instruction;
programcounter:= programcounter + 1;
END ELSE overflow ("Program space");

COMMENT ************************************* model statement ************************************;

INTEGER modelinstructions, repeatinstruction, idtype;

PROCEDURE modelstatement;
INSPECT modeltext DO
IF lastatom == equal THEN nextatom ELSE error (11);
END OTHERWISE error (10);

PROCEDURE lefthandpart;
repeatinstruction:= programcounter - 1; load (smc, 0);
idtype:= 1; IF lastatom == withsign THEN
BEGIN nextatom; expression; END ELSE load (tcv, 1);
END *** lefthandpart ***;

PROCEDURE righthandpart;
BEGIN	term; load (res, 0); termcounter:= 1;
BEGIN nextatom; term; load (res, 0); termcounter:= termcounter + 1; END;
END *** righthandpart ***;

PROCEDURE expression;
BEGIN	IF lastatom == subsign THEN
BEGIN nextatom; term; load (neg, 0); END ELSE
BEGIN IF lastatom == addsign THEN nextatom; term; END;
nextterm;
END *** expression ***;

PROCEDURE nextterm;
nextterm;
END ELSE
IF lastatom == subsign THEN
nextterm;
END *** nextterm ***;

PROCEDURE term; BEGIN factor; nextfactor; END;

PROCEDURE nextfactor;
IF lastatom == mulsign THEN
nextfactor;
END ELSE
IF lastatom == divsign THEN
nextfactor;
END ELSE
IF lastatom == idisign THEN
nextfactor;
END *** nextfactor ***;

PROCEDURE factor; BEGIN primary; nextprimary; END;

PROCEDURE nextprimary;
WHILE lastatom == ttpsign DO

PROCEDURE primary;
IF lastatom == lparent THEN
BEGIN	nextatom; expression;
IF lastatom == rparent THEN nextatom ELSE error (12);
END ELSE INSPECT lastatom
WHEN   identifier DO arithname
WHEN       number DO arithvalue
WHEN functionname DO functiondesignator (lastatom)
WHEN   optionname DO BEGIN error (13); nextatom; END
OTHERWISE BEGIN error (14); nextatom; END;

PROCEDURE arithname;
INSPECT lastatom DO IF idtype EQ 6 THEN
BEGIN	IF IF declared THEN type EQ 6 ELSE FALSE THEN
arithvalue ELSE BEGIN error (15); nextatom; END;
END ELSE BEGIN arithvalue; type:= idtype; END;

PROCEDURE arithvalue;
INSPECT lastatom DO BEGIN load (tav, index); nextatom; END;

PROCEDURE functiondesignator (atom); REF (atomcell) atom;
BEGIN	INTEGER params; params:= 0;
IF nextatom == lparent THEN
BEGIN	nextatom; parameterlist (params);
IF lastatom == rparent THEN nextatom ELSE error (16);
END; INSPECT atom DO
IF params EQ index THEN load (fun, type) ELSE error (17);
END *** functiondesignator ***;

PROCEDURE parameterlist (params); NAME params; INTEGER params;
BEGIN	expression; params:= 1;
WHILE lastatom == comma DO
BEGIN load (stk, 0); nextatom; expression; params:= params + 1; END;
END *** parameterlist ***;

COMMENT ************************************* input statement ************************************;

INTEGER inputinstructions; BOOLEAN three, four, five;

PROCEDURE inputstatement;
INSPECT inputtext DO
BEGIN	inputpart;
WHILE lastatom == comma DO BEGIN nextatom; inputpart; END;
END OTHERWISE error (20);

PROCEDURE inputpart;
IF lastatom == smaller THEN
BEGIN	idtype:= 6; nextatom; expression;
IF lastatom == greater THEN nextatom ELSE error (21); control;
END ELSE INSPECT lastatom
WHEN   identifier DO inputname
WHEN       number DO inputvalue
WHEN functionname DO BEGIN functiondesignator (lastatom); control; END
WHEN   optionname DO BEGIN error (22); nextatom; control; END
OTHERWISE description (FALSE);

PROCEDURE inputname;
INSPECT lastatom DO IF declared THEN inputvalue
ELSE BEGIN load (sav, index); nextatom; type:= 6; END;

PROCEDURE inputvalue;
INSPECT lastatom DO BEGIN load (tav, index); nextatom; control; END;

PROCEDURE control;
IF lastatom == mulsign THEN
program [address]:= programcounter * 100 + jmp;

PROCEDURE description (repeat); BOOLEAN repeat;
IF lastatom == lparent THEN
BEGIN	nextatom; inputstatement;
IF lastatom == rparent THEN nextatom ELSE error (23);
END ELSE
IF lastatom == lbracket THEN
BEGIN	three:= four:= five:= FALSE;
nextatom; variablelist;
IF lastatom == rbracket THEN nextatom ELSE error (24);
END ELSE
BEGIN error (25); nextatom; END;
END *** description ***;

PROCEDURE variablelist;
BEGIN	variable;
WHILE nextatom == comma DO BEGIN nextatom; variable; END;
END *** variablelist ***;

PROCEDURE variable;
INSPECT lastatom WHEN identifier DO
BEGIN	IF declared THEN
BEGIN	IF type EQ 1 OR type EQ 3 THEN
BEGIN	IF type EQ 1 THEN rightcounter:= rightcounter + 1;
load (rav, index); type:= 3; three:= TRUE;
END ELSE
IF type EQ 2 OR type EQ 5 THEN
BEGIN	IF type EQ 2 THEN leftcounter:= leftcounter + 1;
load (lav, index); type:= 5; five:= TRUE;
END ELSE
IF type EQ 4 THEN load (sav, index) ELSE error (26);
END ELSE
BEGIN load (sav, index); type:= 4; four:= TRUE; END;
END OTHERWISE error (27);

IF NOT (three OR five) THEN skip (address) ELSE

BEGIN	INTEGER count; count:= programcounter - address;
END *** skip ***;

COMMENT ************************************* option statement ***********************************;

INTEGER optioninstructions, optionnumber, submodel;

PROCEDURE optionstatement;
BEGIN	INTEGER i; FOR i:= 1 STEP 1 UNTIL optionlimit DO optionwanted [i]:= FALSE;
INSPECT optiontext DO IF lastatom =/= semicolon THEN
BEGIN	option; WHILE lastatom == comma DO BEGIN nextatom; option; END;
END ELSE currenttext:- optiontext:- NONE;
END *** optionstatement ***;

PROCEDURE option;
BEGIN	INSPECT lastatom WHEN optionname DO
BEGIN	optionnumber:= type; optionwanted [type]:= TRUE; END
WHEN number DO
BEGIN	optionnumber:= Scanint (entry); entry . Setpos (1);
IF 1 LE optionnumber AND optionnumber LE optionlimit
THEN optionwanted [optionnumber]:= TRUE ELSE error (30);
END OTHERWISE error (31); IF nextatom == lparent THEN
BEGIN	nextatom; specifier;
WHILE nextatom == comma DO BEGIN nextatom; specifier; END;
IF lastatom == rparent THEN nextatom ELSE error (32);
END;
END *** option ***;

PROCEDURE specifier;
INSPECT lastatom WHEN number DO IF optionnumber EQ 5 THEN
BEGIN	submodel:= Abs (Scanint (entry)); entry . Setpos (1);
IF submodel LE termcounter - 1 THEN load (submodel, 0) ELSE error (33);
END ELSE IF optionnumber EQ 10 THEN
BEGIN	stack [stacklength]:= Scanreal (entry); stacklength:= stacklength - 1;
IF entry . Pos EQ 1 THEN error (33); entry . Setpos (1);
END ELSE error (34) OTHERWISE error (35);

COMMENT ************************************* translation ****************************************;

PROCEDURE translate (sourcetext, statement, restofstatement, instructions); NAME instructions;
REF (textstorage) sourcetext; PROCEDURE statement, restofstatement; INTEGER instructions;
BEGIN	currenttext:- sourcetext; printtext; nextatom; statement; INSPECT currenttext DO
BEGIN	IF lastatom =/= semicolon AND NOT (endline AND endtext) THEN error (36);
WHILE lastatom =/= semicolon AND NOT (endline AND endtext) DO restofstatement;
IF lastatom =/= semicolon THEN error (37);
END; instructions:= programcounter - 1;
END *** translate ***;

PROCEDURE prepare;
BEGIN	IF rightcounter EQ 0 THEN BEGIN rightcounter:= 1; error1 (40, 1); END;
IF leftcounter EQ 0 THEN
designcounter:= termcounter + 1; weightcounter:= termcounter + 2;
END *** prepare ***;

translate (modeltext, modelstatement, expression, modelinstructions);
translate (inputtext, inputstatement, inputstatement, inputinstructions);
translate (optiontext, optionstatement, optionstatement, optioninstructions);
IF erroneous THEN GOTO endjob; prepare;

BEGIN	INTEGER stackpointer, termpointer, casepointer, repeats, repeatsize, designsize;
LONG REAL weight, weightroot, weightsize, repeatroot, repeatsum, repeatsquares;

REF (designstorage) ARRAY designdata [1 : weightcounter];
BEGIN	REF  (rightstorage) ARRAY  rightdata [1 :  rightcounter];
REF   (leftstorage) ARRAY   leftdata [1 :   leftcounter];

COMMENT ************************************* checking *******************************************;

PROCEDURE checkmodel;
BEGIN	INTEGER count, param, term, right, left,
instruction, address, last, next; BOOLEAN constant;
param:= term:= right:= left:= 0; constant:= FALSE;
FOR count:= 3 STEP 1 UNTIL inputinstructions DO
BEGIN	instruction:= program [count]; address:= instruction // 100;
IF instruction GT 100 AND Mod (instruction, 100) EQ tav THEN
BEGIN	INSPECT atomhash . hashtable [address] WHEN identifier DO
IF type EQ 1 THEN
BEGIN	IF term GE 1 THEN
BEGIN	param:= param + 1; stack [address]:= 1;
next:= program [count + 1];
IF last EQ res AND next EQ res THEN constant:= TRUE;
IF last EQ res EQV next EQ mul THEN error1 (41, term);
END ELSE error1 (42, term);
END ELSE
IF type EQ 3 AND stack [address] EQ 0 THEN
BEGIN	right:= right + 1; stack [address]:= right;
END ELSE
IF type EQ 5 AND stack [address] EQ 0 THEN
BEGIN	left:= left + 1; stack [address]:= left;
END ELSE
IF type EQ 2 THEN error1 (43, term)
WHEN number DO IF entry . Pos EQ 1 THEN
IF entry . Pos EQ 1 THEN error1 (44, term);
END;
END ELSE
IF instruction EQ res THEN
BEGIN	IF last NE mul THEN
BEGIN	IF constant THEN constant:= FALSE
ELSE error1 (45, term)
END; IF param GT 1 THEN error1 (46, term)
ELSE IF param EQ 0 THEN error1 (47, term);
param:= 0; term:= term + 1;
END ELSE
IF instruction EQ dep THEN BEGIN instruction:= res; term:= 1; END;
last:= instruction;
END;
IF erroneous THEN GOTO endjob;
IF NOT dataread THEN INSPECT printstream DO
BEGIN	errortext (IF dataskipped THEN 2 ELSE 1); INSPECT datastream DO
BEGIN	printstream . Setpos (14); Outtext ("In record ");
Outint (datarecord, field (datarecord));
Outtext (" of file : "); Outtext (Filename (datastream) . Strip);
END; Outimage; GOTO endjob;
END;
IF NOT optionwanted [7] THEN INSPECT inputdata DO
BEGIN IF optionwanted [6] THEN printdata; reset; END;
stackpointer:= tablelimit + 1; programcounter:= modelinstructions + 1;
END *** checkmodel ***;

PROCEDURE checkinput;
BEGIN	INTEGER i; REAL size, check;
IF NOT optionwanted [7] THEN INSPECT inputdata DO
IF nextnumber NE size - 1 THEN INSPECT printstream DO
BEGIN	Eject (Line + 5); Outtext ("Not all input data (in this");
Outtext (" record) is processed    (message)"); Outimage;
END;
size:= leftdata [1] . size; IF size EQ 0 THEN error2 (70, size, 0);
FOR i:= 2 STEP 1 UNTIL leftcounter DO
BEGIN check:= leftdata [i] . size; IF size NE check THEN error2 (71, size, check); END;
size:= rightdata [1] . size; IF size EQ 0 THEN error2 (72, size, 0);
FOR i:= 2 STEP 1 UNTIL rightcounter DO
BEGIN check:= rightdata [i] . size; IF size NE check THEN error2 (73, size, check); END;
stack [stackpointer]:= size; stackpointer:= stackpointer + 1;
repeatsize:= 0; weightsize:= 0; programcounter:= 1; casepointer:= size + 1;
END *** checkinput ***;

COMMENT ************************************* execution ******************************************;

PROCEDURE execute (instructions); INTEGER instructions;
BEGIN	INTEGER i, j, instruction, address; LONG REAL F, G; BOOLEAN missingvalue, missingcase;
PROCEDURE checkmissing;
FOR j:= stacklength + 1 STEP 1 UNTIL stackspace DO IF stack [j] EQ G THEN missingvalue:= TRUE;
SWITCH macro:=	stk, add, sub, mul, div, idi, ttp, neg,
eql, rep, dep, res, tst, ini, frs, fls, smc;
SWITCH macro2:= tcv, tav, rav, sav, lav, rtn, jmp, fix, skp, fun;
next: instruction:= program [programcounter]; address:= instruction // 100;
programcounter:= programcounter + 1;
GOTO IF instruction LT 100 THEN macro  [instruction]
ELSE macro2 [Mod (instruction, 100)];
stk: IF stackpointer GT stacklength THEN overflow ("Stack space");
stack [stackpointer]:= F; stackpointer:= stackpointer + 1; GOTO exit;
add: stackpointer:= stackpointer - 1; F:= stack [stackpointer] + F; GOTO exit;
sub: stackpointer:= stackpointer - 1; F:= stack [stackpointer] - F; GOTO exit;
mul: stackpointer:= stackpointer - 1; F:= stack [stackpointer] * F; GOTO exit;
div: IF F EQ 0 THEN error2 (50, termpointer, casepointer);
stackpointer:= stackpointer - 1; F:= stack [stackpointer] / F; GOTO exit;
idi: IF F EQ 0 THEN error2 (51, termpointer, casepointer);
stackpointer:= stackpointer - 1; F:= stack [stackpointer] // F; GOTO exit;
ttp: stackpointer:= stackpointer - 1; G:= stack [stackpointer]; I:= F;
IF G EQ 0 AND F LT 0 THEN error2 (54, F, casepointer);
IF G LT 0 AND F NE I THEN error2 (55, F, casepointer);
F:= IF F EQ I THEN G ** I ELSE G ** F; GOTO exit;
neg: F:= - F; GOTO exit;
eql: G:= inputdata . nextnumber;
IF Abs (F - G) GT precision * Abs (F) THEN error2 (74, F, G); GOTO exit;
rep: IF Abs (F) GT maxobs THEN error2 (52, 0, casepointer);
INSPECT designdata [designcounter] DO
BEGIN IF F LT min THEN min:= F; IF F GT max THEN max:= F; END; repeats:= repeats + 1;
repeatsum:= repeatsum + F; repeatsquares:= repeatsquares + F * F; GOTO exit;
dep: IF F LE precision THEN error2 (56, F, casepointer);
weight:= F; weightroot:= Sqrt (weight); repeatroot:= Sqrt (repeats);
repeatsize:= repeatsize + repeats; weightsize:= weightsize + repeats * weight;
designdata [weightcounter] . store (repeatroot * weightroot);
INSPECT designdata [designcounter] DO
BEGIN	store (repeatsum / repeatroot * weightroot);
sum:= sum + repeatsum * weight; squares:= squares + repeatsquares * weight;
END; termpointer:= 1; GOTO exit;
res: IF Abs (F) GT maxobs THEN error2 (53, termpointer, casepointer);
INSPECT designdata [termpointer] DO
BEGIN	store (F * repeatroot * weightroot);
IF F LT min THEN min:= F; IF F GT max THEN max:= F;
G:= F * repeats * weight; sum:= sum + G; squares:= squares + F * G;
END; termpointer:= termpointer + 1; GOTO exit;
tst: F:= leftdata [1] . lastnumber; FOR i:= 2 STEP 1 UNTIL leftcounter DO
BEGIN G:= leftdata [i] . lastnumber; IF F NE G THEN error2 (75, F, G); END;
ini: stack [stackpointer]:= F; stackpointer:= stackpointer + 1;
repeats:= 0; repeatsum:= repeatsquares:= 0; casepointer:= casepointer - 1;
frs: missingvalue:= FALSE; FOR i:= 1 STEP 1 UNTIL rightcounter DO INSPECT rightdata [i] DO
BEGIN G:= stack [address]:= lastnumber; IF optionwanted [10] THEN checkmissing; END;
missingcase:= missingvalue; GOTO exit;
fls: missingvalue:= FALSE; FOR i:= 1 STEP 1 UNTIL leftcounter DO INSPECT leftdata [i] DO
BEGIN G:= stack [address]:= lastnumber; IF optionwanted [10] THEN checkmissing; END;
IF missingcase OR missingvalue THEN programcounter:= repeatinstruction; GOTO exit;
smc: IF missingcase OR repeats EQ 0 THEN programcounter:= modelinstructions; GOTO exit;
tav: F:= stack [address]; GOTO exit;
rav: rightdata [stack [address]] . store (inputdata . nextnumber); GOTO exit;
sav: stack [address]:= inputdata . nextnumber; GOTO exit;
lav: leftdata [stack [address]] . store (inputdata . nextnumber); GOTO exit;
rtn: F:= stack [stackpointer - 1]; IF F LT 1.5 THEN stackpointer:= stackpointer - 1 ELSE
BEGIN stack [stackpointer - 1]:= F - 1; programcounter:= address; END; GOTO exit;
jmp: IF Abs (F - Entier (F + precision * Abs (F))) GT 2 * precision * Abs (F) THEN error2
(76, F, IF stackpointer EQ tablelimit + 1 THEN 0 ELSE stack [stackpointer - 1] - 1);
IF F LT 0.5 THEN programcounter:= address; GOTO exit;
fix: FOR i:= programcounter - 2 STEP - 1 UNTIL address DO
IF Mod (program [i], 100) EQ 5 THEN leftdata [stack [program [i] // 100]] . fixrepeats; GOTO exit;
skp: inputdata . skipnumbers (F * address); GOTO exit;
fun: BEGIN	SWITCH function:= fabs, fsign, fsqrt, fsin, fcos, ftan, fln, flog, fexp,
fent, frnd, fmod, fmin, fmax, fasin, facos, fatan, fsinh, fcosh, ftanh, find;
fabs: F:= Abs (F); GOTO exit;
fsign: F:= Sign (F); GOTO exit;
fsqrt: IF F LT 0 THEN error2 (60, F, casepointer);
F:= Sqrt (F); GOTO exit;
fsin: F:= Sin (F); GOTO exit;
fcos: F:= Cos (F); GOTO exit;
ftan: F:= Tan (F); GOTO exit;
fln: IF F LE 0 THEN error2 (61, F, casepointer);
F:= Ln (F); GOTO exit;
flog: IF F LE 0 THEN error2 (62, F, casepointer);
F:= Ln (F) / Ln10; GOTO exit;
fexp: IF F GT maxexp THEN error2 (63, F, casepointer);
F:= IF F LT - maxexp THEN 0 ELSE Exp (F); GOTO exit;
fent: F:= Entier (F); GOTO exit;
frnd: F:= Entier (F + 0.5); GOTO exit;
fmod: stackpointer:= stackpointer - 1;
F:= Mod (stack [stackpointer], F); GOTO exit;
fmin: stackpointer:= stackpointer - 1;
F:= Lmin (stack [stackpointer], F); GOTO exit;
fmax: stackpointer:= stackpointer - 1;
F:= Lmax (stack [stackpointer], F); GOTO exit;
fasin: IF Abs (F) GT 1 THEN error2 (64, F, casepointer);
F:= Arcsin (F); GOTO exit;
facos: IF Abs (F) GT 1 THEN error2 (65, F, casepointer);
F:= Arccos (F); GOTO exit;
fatan: F:= Arctan (F); GOTO exit;
fsinh: IF Abs (F) GT maxexp THEN error2 (66, F, casepointer);
F:= Sinh (F); GOTO exit;
fcosh: IF Abs (F) GT maxexp THEN error2 (67, F, casepointer);
F:= Cosh (F); GOTO exit;
ftanh: F:= Tanh (F); GOTO exit;
find: stackpointer:= stackpointer - 2; G:= stack [stackpointer + 1];
F:= IF stack [stackpointer] - G LE precision * Abs (G) AND
G - F LE precision * Abs (G) THEN 1 ELSE 0; GOTO exit;
END *** fun ***;
exit: IF programcounter LE instructions THEN GOTO next;
END *** execute ***;

PROCEDURE initializestorage;
BEGIN	INTEGER i;
FOR i:= 1 STEP 1 UNTIL weightcounter DO designdata [i]:- NEW designstorage;
FOR i:= 1 STEP 1 UNTIL  rightcounter DO  rightdata [i]:- NEW  rightstorage;
FOR i:= 1 STEP 1 UNTIL   leftcounter DO   leftdata [i]:- NEW   leftstorage;
END *** initializestorage ***;

initializestorage;
checkmodel; execute (inputinstructions);
checkinput; execute (modelinstructions);
designsize:= designdata [1] . size;

END *** execution ***;

BEGIN	INTEGER ARRAY  pivot [1 : termcounter];
LONG REAL ARRAY diag [1 : termcounter], resprod [1 : designcounter],
data [1 : designsize, 1 : weightcounter], depvar [1 : designsize];

COMMENT ************************************* print layout ***************************************;

INTEGER low, upp, max, pagenr, partnr, matrixwidth;

PROCEDURE printname (col); INTEGER col;
INSPECT printstream DO
IF col EQ weightcounter THEN
BEGIN	IF repeatsize EQ weightsize THEN
Outtext (" repeats") ELSE Outtext ("  weight");
END ELSE IF col EQ designcounter THEN Outtext ("dep.var.")
ELSE printitem (atomhash . hashtable [designdata [col] . address]);

PROCEDURE underline (width); INTEGER width;
INSPECT printstream DO BEGIN WHILE Pos LE width DO Outchar ('='); Outimage; END;

PROCEDURE dashedline;
INSPECT printstream DO
BEGIN	Eject (Line + 1); Outtext ("----------------------------------");
Outtext ("----------------------------------------------------");
Outtext ("-------------------------"); Outimage; Eject (Line + 1);
END *** dashedline ***;

INSPECT printstream DO
BEGIN	INTEGER j; Setpos (14);
FOR j:= low STEP 1 UNTIL upp DO
BEGIN printname (j); Setpos (Pos // 12 * 12 + 14); END;
Outimage; Eject (Line + 1);

PROCEDURE newpage1;
INSPECT printstream DO
BEGIN	Eject (1); Outtext ("Transformed data matrix");
IF designsize GT 50 OR designcounter GT 10 THEN
BEGIN	pagenr:= pagenr + 1;
Outtext ("  -  page"); Outint (pagenr, 3);
END; Outimage; underline (23); Eject (Line + 1);
END *** newpage1 ***;

PROCEDURE newpart (string); NAME string; TEXT string;
INSPECT printstream DO
BEGIN	IF Line GT 49 - termcounter + low
THEN Eject (1) ELSE Eject (Line + 5);
Outtext ("Correlation matrix of the "); Outtext (string);
IF termcounter GT 10 AND designcounter GT 10 THEN
BEGIN	partnr:= partnr + 1;
Outtext ("  -  part"); Outint (partnr, 3);
END; Outimage; underline (35); Eject (Line + 1); heading;
END *** newpart ***;

PROCEDURE newpage2;
INSPECT printstream DO
BEGIN	Eject (1); Outtext ("Residual analysis");
IF designsize GT 50 THEN
BEGIN	pagenr:= pagenr + 1;
Outtext ("  -  page"); Outint (pagenr, 3);
END; Outimage; underline (17); Setpos (95);
Outtext ("standardized         studentized"); Outimage;
Outtext ("obs.no.        observation        fitted value");
Outtext ("    standard deviation          residual");
Outtext ("            residual            residual");
Outimage; Eject (Line + 1);
END *** newpage2 ***;

COMMENT ************************************* Dekker matrix routines *****************************;

LONG REAL PROCEDURE vecvec (l, u, a, b);
INTEGER l, u; LONG REAL ARRAY a, b;
BEGIN	INTEGER k; LONG REAL s; s:= 0;
FOR k:= l STEP 1 UNTIL u DO s:= a [k] * b [k] + s;
vecvec:= s;
END *** vecvec ***;

LONG REAL PROCEDURE matvec (l, u, i, a, b);
INTEGER l, u, i; LONG REAL ARRAY a, b;
BEGIN	INTEGER k; LONG REAL s; s:= 0;
FOR k:= l STEP 1 UNTIL u DO s:= a [i,k] * b [k] + s;
matvec:= s;
END *** matvec ***;

LONG REAL PROCEDURE tamvec (l, u, i, a, b);
INTEGER l, u, i; LONG REAL ARRAY a, b;
BEGIN	INTEGER k; LONG REAL s; s:= 0;
FOR k:= l STEP 1 UNTIL u DO s:= a [k,i] * b [k] + s;
tamvec:= s;
END *** tamvec ***;

LONG REAL PROCEDURE tammat (l, u, i, j, a, b);
INTEGER l, u, i, j; LONG REAL ARRAY a, b;
BEGIN	INTEGER k; LONG REAL s; s:= 0;
FOR k:= l STEP 1 UNTIL u DO s:= a [k,i] * b [k,j] + s;
tammat:= s;
END *** tammat ***;

PROCEDURE elmcol (l, u, i, j, a, b, x);
INTEGER l, u, i, j; LONG REAL x; LONG REAL ARRAY a, b;
FOR l:= l STEP 1 UNTIL u DO a [l,i]:= a [l,i] + b [l,j] * x;

PROCEDURE elmveccol (l, u, i, a, b, x);
INTEGER l, u, i; LONG REAL x; LONG REAL ARRAY a, b;
FOR l:= l STEP 1 UNTIL u DO a [l]:= a [l] + b [l,i] * x;

PROCEDURE ichcol (l, u, i, j, a);
INTEGER l, u, i, j; LONG REAL ARRAY a;
BEGIN	LONG REAL r;
FOR l:= l STEP 1 UNTIL u DO
BEGIN r:= a [l,i]; a [l,i]:= a [l,j]; a [l,j]:= r; END;
END *** ichcol ***;

PROCEDURE ichrow (l, u, i, j, a);
INTEGER l, u, i, j; LONG REAL ARRAY a;
BEGIN	LONG REAL r;
FOR l:= l STEP 1 UNTIL u DO
BEGIN r:= a [i,l]; a [i,l]:= a [j,l]; a [j,l]:= r; END;
END *** ichrow ***;

PROCEDURE ichrowcol (l, u, i, j, a);
INTEGER l, u, i, j; LONG REAL ARRAY a;
BEGIN	LONG REAL r;
FOR l:= l STEP 1 UNTIL u DO
BEGIN r:= a [i,l]; a [i,l]:= a [l,j]; a [l,j]:= r; END;
END *** ichrowcol ***;

PROCEDURE chlinv2 (a, n); INTEGER n; LONG REAL ARRAY a;
BEGIN	LONG REAL r; INTEGER i, j, i1;
LONG REAL ARRAY u [1 : n];
FOR i:= n STEP - 1 UNTIL 1 DO
BEGIN	r:= 1 / a [i,i]; i1:= i + 1;
FOR j:= i1 STEP 1 UNTIL n DO u [j]:= a [i,j];
FOR j:= n STEP - 1 UNTIL i1 DO
a [i,j]:= - (tamvec (i1, j, j, a, u) + matvec (j + 1, n, j, a, u)) * r;
a [i,i]:= (r - matvec (i1, n, i, a, u)) * r;
END;
END *** chlinv2 ***;

INTEGER PROCEDURE lsqdec (a, n, m, tol, aid, ci);
INTEGER n, m; LONG REAL tol; LONG REAL ARRAY a, aid; INTEGER ARRAY ci;
BEGIN	INTEGER j, k, kpiv;
LONG REAL beta, sigma, norm, w, eps, akk, aidk;
LONG REAL ARRAY sum [1 : m];
norm:= 0; lsqdec:= m;
FOR k:= 1 STEP 1 UNTIL m DO
BEGIN	w:= sum [k]:= tammat (1, n, k, k, a, a);
IF w GT norm THEN norm:= w;
END;
w:= Sqrt (norm); eps:= tol * w;
FOR k:= 1 STEP 1 UNTIL m DO
BEGIN	sigma:= sum [k]; kpiv:= k;
FOR j:= k + 1 STEP 1 UNTIL m DO
IF sum [j] GT sigma THEN
BEGIN sigma:= sum [j]; kpiv:= j; END;
IF kpiv NE k THEN
BEGIN sum [kpiv]:= sum [k]; ichcol (1, n, k, kpiv, a); END;
ci [k]:= kpiv; akk:= a [k,k];
sigma:= tammat (k, n, k, k, a, a); w:= Sqrt (sigma);
aidk:= aid [k]:= IF akk LT 0 THEN w ELSE - w;
IF w LT eps THEN BEGIN lsqdec:= k - 1; GOTO exit; END;
beta:= 1 / (sigma - akk * aidk); a [k,k]:= akk - aidk;
FOR j:= k + 1 STEP 1 UNTIL m DO
BEGIN	elmcol (k, n, j, k, a, a, - beta * tammat (k, n, k, j, a, a));
sum [j]:= sum [j] - a [k,j] ** 2;
END;
END;	exit:
END *** lsqdec ***;

PROCEDURE lsqsol (a, n, m, aid, ci, b);
INTEGER n, m; LONG REAL ARRAY a, aid, b; INTEGER ARRAY ci;
BEGIN	INTEGER k, cik; LONG REAL w;
FOR k:= 1 STEP 1 UNTIL m DO
elmveccol (k, n, k, b, a, tamvec (k, n, k, a, b) / (aid [k] * a [k,k]));
FOR k:= m STEP - 1 UNTIL 1 DO
b [k]:= (b [k] - matvec (k + 1, m, k, a, b)) / aid [k];
FOR k:= m STEP - 1 UNTIL 1 DO
BEGIN	cik:= ci [k]; IF cik NE k THEN
BEGIN w:= b [k]; b [k]:= b [cik]; b [cik]:= w; END;
END;
END *** lsqsol ***;

PROCEDURE lsqinv (a, m, aid, c);
INTEGER m; LONG REAL ARRAY a, aid; INTEGER ARRAY c;
BEGIN	INTEGER i, ci; LONG REAL w;
FOR i:= 1 STEP 1 UNTIL m DO a [i,i]:= aid [i];
chlinv2 (a, m);
FOR i:= m STEP - 1 UNTIL 1 DO
BEGIN	ci:= c [i]; IF ci NE i THEN
BEGIN	ichcol (1, i - 1, i, ci, a);
ichrowcol (i + 1, ci - 1, i, ci, a);
ichrow (ci + 1, m, i, ci, a);
w:= a [i,i]; a [i,i]:= a [ci,ci]; a [ci,ci]:= w;
END;
END;
END *** lsqinv ***;

COMMENT ************************************* distribution functions *****************************;

REAL PROCEDURE Phi (x); REAL x;
BEGIN	REAL y, z, w; y:= Abs (x) / 2;
IF y GE 3 THEN z:= 1 ELSE IF y LT 1 THEN
BEGIN	w:= y * y;
z:= (((((((( 0.000124818987 * w
- 0.001075204047) * w + 0.005198775019) * w
- 0.019198292004) * w + 0.059054035642) * w
- 0.151968751364) * w + 0.319152932694) * w
- 0.531923007300) * w + 0.797884560593) * y * 2
END ELSE
BEGIN	y:= y - 2;
z:= ((((((((((((( - 0.000045255659 * y
+ 0.000152529290) * y - 0.000019538132) * y
- 0.000676904986) * y + 0.001390604284) * y
- 0.000794620820) * y - 0.002034254874) * y
+ 0.006549791214) * y - 0.010557625006) * y
+ 0.011630447319) * y - 0.009279453341) * y
+ 0.005353579108) * y - 0.002141268741) * y
+ 0.000535310849) * y + 0.999936657524
END;
Phi:= IF x GT 0 THEN (z + 1) / 2 ELSE (1 - z) / 2
END *** Phi ***;

REAL PROCEDURE Fisher (x, df1, df2); REAL x; INTEGER df1, df2;
IF x LE 0 THEN Fisher:= 0 ELSE
IF x GE 10000 THEN Fisher:= 1 ELSE
IF df1 GT 1000 AND df2 GT 1000 THEN
Fisher:= Phi ((x ** (1/3) * (1 - 2/df2/9) + 2/df1/9 - 1) /
Sqrt (x ** (2/3) *      2/df2/9  + 2/df1/9)) ELSE
BEGIN	INTEGER a, b, i, j; REAL w, y, z, zk, d, p;
a:= df1 // 2 * 2 - df1 + 2; b:= df2 // 2 * 2 - df2 + 2;
w:= x * df1 / df2; z:= 1 / (1 + w);
IF a EQ 1 THEN
BEGIN	IF b EQ 1 THEN
BEGIN	p:= Sqrt (w); y:= 0.318309886184;
d:= y * z / p; p:= Arctan (p) * 2 * y;
END ELSE
BEGIN p:= Sqrt (w * z); d:= 0.5 * p * z / w; END;
END ELSE
IF b EQ 1 THEN
BEGIN p:= Sqrt (z); d:= 0.5 * z * p; p:= 1 - p; END
ELSE BEGIN d:= z * z; p:= w * z; END;
y:= 2 * w / z;
IF a EQ 1 THEN
BEGIN	FOR j:= b + 2 STEP 2 UNTIL df2 DO
BEGIN	d:= (a / (j - 2) + 1) * d * z;
p:= d * y / (j - 1) + p;
END;
END ELSE
BEGIN	zk:= z ** ((df2 - 1) // 2); d:= d * zk * df2 / b;
p:= p * zk + w * z * (zk - 1) / (z - 1);
END;
y:= w * z; z:= 2 / z; b:= df2 - 2;
FOR i:= a + 2 STEP 2 UNTIL df1 DO
BEGIN	j:= i + b; d:= y * d * j / (i - 2);
p:= p - z * d / j;
END;
Fisher:= IF p LT 0 THEN 0 ELSE IF p GT 1 THEN 1 ELSE p;
END *** Fisher ***;

COMMENT ************************************* transformed data matrix ****************************;

PROCEDURE transformed_data_matrix;
IF optionwanted [1] AND submodel EQ 0 THEN
BEGIN	INTEGER i, j, datacounter; datacounter:=
IF designsize EQ weightsize THEN designcounter ELSE weightcounter;
INSPECT outputstream DO
BEGIN	Outint (designsize, 12); Outint (datacounter, 12); Outimage;
FOR i:= 1 STEP 1 UNTIL designsize DO
BEGIN	weightroot:= data [i, weightcounter];
FOR j:= 1 STEP 1 UNTIL datacounter DO IF j NE weightcounter
THEN Outreal (data [i,j] / weightroot, 16, 22)
ELSE Outreal (weightroot * weightroot, 16, 22); Outimage;
END;
END;
pagenr:= 0; upp:= 0;
INSPECT printstream DO
FOR low:= upp + 1 WHILE upp LT datacounter DO
BEGIN	upp:= IF datacounter - low GT matrixwidth
THEN low + matrixwidth ELSE datacounter; newpage1;
FOR i:= 1 STEP 1 UNTIL designsize DO
BEGIN	Outint (i, 6); Setpos (10);
weightroot:= data [i, weightcounter];
FOR j:= low STEP 1 UNTIL upp DO IF j NE weightcounter
THEN Outfix (data [i,j] / weightroot, 3, 12)
ELSE Outfix (weightroot * weightroot, 3, 12); Outimage;
IF Mod (i, 10) EQ 0 AND i LT designsize THEN
BEGIN IF Mod (i, 50) EQ 0 THEN newpage1 ELSE Eject (Line + 1); END;
END;
END;
END *** transformed_data_matrix ***;

COMMENT ************************************* control information ********************************;

INTEGER minterm, constterm; LONG REAL standdev;
BOOLEAN constant, depvarconstant, oneconstant, moreconstants;

PROCEDURE control_information;
INSPECT printstream DO
BEGIN	INTEGER j; constterm:= 0;
constant:= depvarconstant:= oneconstant:= moreconstants:= FALSE;
Eject (1); Outtext ("Control information");
IF submodel GT 0 THEN
BEGIN Outtext ("  -  submodel"); Outint (submodel, 3); END;
Outimage; underline (19); Eject (Line + 1);
Outtext ("transformed variable"); Outimage;
Outtext ("denoted by parameter          mean           ");
Outtext ("standard deviation                minimum    ");
Outtext ("              maximum"); Outimage; Eject (Line + 1);
FOR j:= termcounter + 1 STEP 1 UNTIL designcounter - 1 DO
BEGIN	Setpos (2); printname (j); Setpos (14);
Outtext ("omitted"); Outimage;
END;
IF submodel GT 0 THEN Eject (Line + 1);
FOR j:= 1 STEP 1 UNTIL termcounter, designcounter DO
INSPECT designdata [j] DO
BEGIN	Setpos (2); printname(j); Setpos (12);
standdev:= squares - sum * sum / weightsize;
standdev:= IF standdev LT precision * squares OR weightsize LE 1
THEN 0 ELSE Sqrt (standdev / (weightsize - 1));
Outfix (sum / weightsize, 6, 25); Outfix (standdev, 6, 25);
Outfix (min, 6, 25); Outfix (max, 6, 25); Outimage;
IF standdev EQ 0 THEN
BEGIN	depvarconstant:= j EQ designcounter;
oneconstant:= termcounter EQ 1; constterm:= j;
IF NOT depvarconstant THEN
BEGIN moreconstants:= constant; constant:= TRUE; END;
END;
END;
Eject (Line + 1); Outtext ("Number of observations :  ");
Outint (designsize, Imax (field (designsize), 3)); Outimage;
minterm:= IF constterm EQ 1 THEN 2 ELSE 1;
IF NOT (constant OR optionwanted [4]) THEN
BEGIN	Eject (Line + 5); Outtext ("There is no constant independent");
Outtext (" variable in the transformed (sub)model    (message)"); Outimage;
END;
END *** control_information ***;

COMMENT ************************************* correlation matrix *********************************;

PROCEDURE correlation_matrix;
IF optionwanted [2] AND submodel EQ 0 THEN
BEGIN	INTEGER i, j, pj; LONG REAL resdevj, correlation;
LONG REAL ARRAY resdev [1 : designcounter],
correl [1 : (designcounter + 1) * designcounter // 2];
FOR j:= 1 STEP 1 UNTIL designcounter DO
INSPECT designdata [j] DO
BEGIN	resdevj:= squares - sum * sum / weightsize;
resdevj:= resdev [j]:= IF resdevj LT precision * squares
THEN 0 ELSE Sqrt (resdevj * weightsize);
pj:= (j - 1) * j // 2; correl [pj + j]:= 1;
FOR i:= j - 1 STEP - 1 UNTIL 1 DO correl [pj + i]:=
IF resdev [i] EQ 0 OR resdevj EQ 0 THEN Maxreal ELSE
(tammat (1, designsize, i, j, data, data) * weightsize
- designdata [i] . sum * sum) / (resdev [i] * resdevj);
END;
INSPECT outputstream DO
BEGIN	Outint (designcounter, 12); Outimage;
FOR j:= 1 STEP 1 UNTIL designcounter DO
BEGIN	pj:= (j - 1) * j // 2;
FOR i:= 1 STEP 1 UNTIL j DO
Outreal (correl [pj + i], 16, 22);
END; Outimage;
END;
partnr:= 0; upp:= 0;
INSPECT printstream DO
FOR low:= upp + 1 WHILE upp LT designcounter DO
BEGIN	upp:= IF designcounter - low GT matrixwidth
THEN low + matrixwidth ELSE designcounter;
newpart ("variables");
FOR j:= low STEP 1 UNTIL designcounter DO
BEGIN	Setpos (2); printname (j); Setpos (10);
pj:= (j - 1) * j // 2; max:= Imin (j, upp);
FOR i:= low STEP 1 UNTIL max DO
BEGIN	correlation:= correl [pj + i];
IF Abs (correlation) LT 1.0000005
THEN Outfix (correlation, 6, 12)
ELSE Outtext ("        *   ");
END; Outimage;
END;
END;
END *** correlation_matrix ***;

COMMENT ************************************* no regression analysis *****************************;

PROCEDURE no_regression_analysis;
IF optionwanted [4] OR designsize LE termcounter
OR depvarconstant OR oneconstant OR moreconstants THEN
INSPECT printstream DO
BEGIN	Eject (Line + 5); Outtext ("No regression analysis");
IF NOT optionwanted [4] THEN Outtext (" possible") ELSE
BEGIN Outtext (" by option"); Outimage; GOTO endjob; END; Outimage;
IF depvarconstant THEN
BEGIN	Eject (Line + 2); Outtext ("The dependent");
Outtext (" variable in the transformed model is constant"); Outimage;
END;
IF oneconstant THEN
BEGIN	Eject (Line + 2); Outtext ("The only independent");
Outtext (" variable in the transformed model is constant"); Outimage;
END;
IF moreconstants THEN
BEGIN	Eject (Line + 2); Outtext ("There are several constant");
Outtext (" independent variables in the transformed model"); Outimage;
END;
IF designsize LE termcounter THEN
BEGIN	Eject (Line + 2); Outtext ("The number of respondents");
Outtext (" is less than or equal to the number of");
Outtext (" independent variables in the transformed model"); Outimage;
END;
erroneous:= TRUE; GOTO endjob;
END *** no_regression_analysis ***;

COMMENT ************************************* least squares **************************************;

INTEGER dftotal, dfmean, dfreg, dfres, dflack, dferror, dfred;
LONG REAL sstotal, ssmean, ssreg, ssres, sslack, sserror, ssred, explained, adjusted,
msreg, msres, mslack, mserror, msred, frpar, frmean, frreg, frlack, frred;
BOOLEAN	perfectfit, lackoffit, reduction;

PROCEDURE least_squares;
BEGIN	INTEGER rank;
rank:= lsqdec (data, designsize, termcounter, precision, diag, pivot);
IF rank LT termcounter THEN INSPECT printstream DO
BEGIN	Eject (Line + 5); Outtext ("No regression analysis possible"); Outimage;
Eject (Line + 2); Outtext ("The rank of the transformed data matrix is ");
Outint (termcounter - rank, field (termcounter - rank));
Outtext (" less than the number of independent variables in the transformed model");
Outimage; erroneous:= TRUE; GOTO endjob;
END;
lsqsol (data, designsize, termcounter, diag, pivot, depvar);
lsqinv (data, termcounter, diag, pivot);
perfectfit:= lackoffit:= FALSE;
dftotal:= repeatsize; sstotal:= designdata [designcounter] . squares;
dfreg:= termcounter; ssreg:= vecvec (1, termcounter, depvar, resprod);
dfres:= dftotal - dfreg; ssres:= sstotal - ssreg;
IF ssres LE precision * sstotal THEN BEGIN ssres:= 0; perfectfit:= TRUE; END;
IF constant THEN
BEGIN	ssmean:= designdata [designcounter] . sum ** 2 / weightsize;
dfmean:= 1; dfreg:= dfreg - dfmean; ssreg:= ssreg - ssmean;
END ELSE BEGIN dfmean:= 0; ssmean:= 0; END;
explained:= 1 - ssres / (sstotal - ssmean);
adjusted:= 1 - ssres / (sstotal - ssmean) * (dftotal - dfmean) / dfres;
msreg:= ssreg / dfreg; msres:= ssres / dfres;
frmean:= IF perfectfit THEN Maxreal ELSE ssmean / msres;
frreg:= IF perfectfit THEN Maxreal ELSE msreg / msres;
dferror:= repeatsize - designsize;
IF dferror GT 0 THEN
BEGIN	sserror:= sstotal - resprod [designcounter];
IF sserror LE precision * sstotal THEN BEGIN sserror:= 0; lackoffit:= TRUE; END;
sslack:= ssres - sserror; IF sslack LE precision * ssres THEN sslack:= 0;
dflack:= dfres - dferror; mslack:= sslack / dflack; mserror:= sserror / dferror;
frlack:= IF perfectfit OR lackoffit THEN 0 ELSE mslack / mserror;
END;
IF submodel EQ 0 AND (optionwanted [5] OR optionwanted [8]) THEN
BEGIN	savemodel:= TRUE; savefit:= perfectfit;
savedf:= dfres; savess:= ssres; savems:= msres;
END;
reduction:= IF (submodel GT 0 OR optionwanted [9]) AND savemodel
THEN dfres GT savedf AND ssres GT savess ELSE FALSE;
IF reduction THEN
BEGIN	dfred:= dfres - savedf; ssred:= ssres - savess; msred:= ssred / dfred;
frred:= IF savefit THEN Maxreal ELSE msred / savems;
END;
END *** least_squares ***;

COMMENT ************************************* multiple correlation *******************************;

PROCEDURE multiple_correlation;
INSPECT printstream DO
BEGIN	Eject (Line + 5); Outtext ("Multiple correlation coefficient");
Outfix (Sqrt (explained), 6, 13); Outtext ("    (adjusted");
Outfix (Sqrt ( adjusted), 6, 11); Outchar (')'); Outimage; underline (32);
END *** multiple_correlation ***;

COMMENT ************************************* explained variation ********************************;

PROCEDURE explained_variation;
INSPECT printstream DO
BEGIN	Eject (Line + 5); Outtext ("Proportion of variation explained");
Outfix (explained, 6, 12); Outtext ("    (adjusted");
Outfix ( adjusted, 6, 11); Outchar (')'); Outimage; underline (33);
END *** explained_variation ***;

COMMENT ************************************* standard error deviation ***************************;

PROCEDURE standard_error_deviation;
INSPECT printstream DO
BEGIN	Eject (Line + 5); Outtext ("Standard deviation of the error term");
Outfix (Sqrt (msres), 6, 22); Outimage; underline (36);
END *** standard_error_deviation ***;

COMMENT ************************************* regression parameters ******************************;

PROCEDURE regression_parameters;
INSPECT printstream DO
BEGIN	INTEGER j;
INSPECT outputstream DO
BEGIN	IF submodel EQ 0 THEN
BEGIN Outint (subselection + 1, 12); Outimage; END;
IF submodel EQ 0 OR subselection GT 0 THEN
BEGIN	Outint (termcounter, 12); Outimage;
FOR j:= 1 STEP 1 UNTIL termcounter DO
Outreal (depvar [j], 16, 22); Outimage;
END;
END;
Eject (1); Outtext ("Regression parameters"); Outimage;
underline (21); Setpos (101); Outtext ("right  tail"); Outimage;
Outtext ("parameter                 estimate           ");
Outtext ("standard deviation              F - ratio    ");
Outtext ("          probability"); Outimage; Eject (Line + 1);
FOR j:= 1 STEP 1 UNTIL termcounter DO
BEGIN	Setpos (2); printname (j); Setpos (12); standdev:= Sqrt (data [j,j] * msres);
frpar:= IF perfectfit THEN Maxreal ELSE depvar [j] ** 2 / standdev / standdev;
Outfix (depvar [j], 10, 25); Outfix (standdev, 10, 25); Outfix (frpar, 6, 25);
Outfix (1 - Fisher (frpar, 1, dfres), 6, 25); Outimage; Eject (Line + 1);
END;
END *** regression_parameters ***;

COMMENT ************************************* covariance matrix **********************************;

PROCEDURE covariance_matrix;
IF optionwanted [2] AND (submodel EQ 0 OR subselection GT 0) THEN
BEGIN	INTEGER i, j; LONG REAL datajj;
INSPECT outputstream DO
BEGIN	Outint (termcounter, 12); Outimage;
FOR j:= 1 STEP 1 UNTIL termcounter DO
FOR i:= 1 STEP 1 UNTIL j DO
Outreal (data [i,j], 16, 22); Outimage;
END;
partnr:= 0; upp:= 0;
INSPECT printstream DO
FOR low:= upp + 1 WHILE upp LT termcounter DO
BEGIN	upp:= IF termcounter - low GT matrixwidth
THEN low + matrixwidth ELSE termcounter;
newpart ("estimates");
FOR j:= low STEP 1 UNTIL termcounter DO
BEGIN	Setpos (2); printname (j); Setpos (10);
datajj:= data [j,j]; max:= Imin (j, upp);
FOR i:= low STEP 1 UNTIL max DO
Outfix (data [i,j] / Sqrt (data [i,i] * datajj), 6, 12);
Outimage;
END;
END;
END *** covariance_matrix ***;

COMMENT ************************************* analysis of variance *******************************;

PROCEDURE analysis_of_variance;
INSPECT printstream DO
BEGIN	Eject (Line + 5); Outtext ("Analysis of variance"); Outimage;
underline (20); Eject (Line + 1); Outtext ("source of");
Setpos (101); Outtext ("right  tail"); Outimage;
Outtext ("variation          df          sum of squares");
Outtext ("           mean square             F - ratio ");
Outtext ("          probability"); Outimage; dashedline;
Outtext ("total"); Setpos (16); Outint (dftotal, 6);
Outfix (sstotal, 6, 24); Outimage; dashedline;
IF constant THEN
BEGIN	Outtext ("mean"); Setpos (16); Outint (dfmean, 6);
Outfix (ssmean, 6, 24); Outfix (ssmean, 6, 22); Outfix (frmean, 6, 22);
Outfix (1 - Fisher (frmean, dfmean, dfres), 6, 22); Outimage;
END;
Outtext ("regression"); Setpos (16); Outint (dfreg, 6);
Outfix (ssreg, 6, 24); Outfix (msreg, 6, 22); Outfix (frreg, 6, 22);
Outfix (1 - Fisher (frreg, dfreg, dfres), 6, 22); Outimage;
Outtext ("residual"); Setpos (16); Outint (dfres, 6);
Outfix (ssres, 6, 24); Outfix (msres, 6, 22); Outimage; dashedline;
IF dferror GT 0 THEN
BEGIN	Outtext (" lack of fit"); Setpos (16); Outint (dflack, 6);
Outfix (sslack, 6, 24); Outfix (mslack, 6, 22); Outfix (frlack, 6, 22);
Outfix (1 - Fisher (frlack, dflack, dferror), 6, 22); Outimage;
Outtext (" pure error"); Setpos (16); Outint (dferror, 6);
Outfix (sserror, 6, 24); Outfix (mserror, 6, 22); Outimage; dashedline;
END;
IF reduction THEN
BEGIN	Outtext (" reduction"); Setpos (16); Outint (dfred, 6);
Outfix (ssred, 6, 24); Outfix (msred, 6, 22); Outfix (frred, 6, 22);
Outfix (1 - Fisher (frred, dfred, savedf), 6, 22); Outimage; dashedline;
END;
END *** analysis_of_variance ***;

COMMENT ************************************* null hypotheses ************************************;

PROCEDURE null_hypotheses;
INSPECT printstream DO
BEGIN	INTEGER j;
Eject (Line + 1); Outtext ("regression null hypothesis :  ");
FOR j:= 1 STEP 1 UNTIL termcounter DO
BEGIN	IF j NE constterm THEN
BEGIN printname (j); Outtext (" = "); END;
IF Mod (j, 8) EQ 0 AND j LT termcounter THEN
BEGIN Outimage; Setpos (31); END;
END;
Outchar ('0'); IF submodel GT 0 THEN Outtext ("  (in the reduced model)"); Outimage;
IF submodel GT 0 AND reduction THEN
BEGIN	Eject (Line + 1); Outtext (" reduction null hypothesis :  ");
FOR j:= termcounter + 1 STEP 1 UNTIL designcounter - 1 DO
BEGIN	printname (j); Outtext (" = ");
IF Mod (j - termcounter, 8) EQ 0 AND j LT designcounter - 1 THEN
BEGIN Outimage; Setpos (31); END;
END;
Outtext ("0  (in the original model)"); Outimage;
END;
IF (submodel GT 0 OR optionwanted [9]) AND NOT reduction THEN
BEGIN	Eject (Line + 5); IF NOT savemodel
THEN Outtext ("No reduction test possible, the original model was not saved")
ELSE Outtext ("No increase in residual degrees of freedom or sum of squares");
Outtext ("    (message)"); Outimage;
END;
END *** null_hypotheses ***;

COMMENT ************************************* residual analysis **********************************;

PROCEDURE residual_analysis;
IF optionwanted [3] AND (submodel EQ 0 OR subselection GT 0) THEN
INSPECT printstream DO
BEGIN	INTEGER i, j, studnum; LONG REAL ARRAY resrow, resdev [1 : termcounter];
LONG REAL observation, fittedvalue, residual, residualsum,
standres, studdev, studres, maxstud, dfstud, frstud, studprob, sqrtmsres;
resetdesigndata; pagenr:= 0; newpage2;
studnum:= 0; residualsum:= 0; maxstud:= - Maxreal;
sqrtmsres:= Sqrt (msres * (designsize - termcounter) / designsize);
FOR i:= 2 STEP 1 UNTIL termcounter DO
FOR j:= i - 1 STEP - 1 UNTIL 1 DO data [i,j]:= data [j,i];
INSPECT outputstream DO BEGIN Outint (designsize, 12); Outimage; END;
FOR i:= 1 STEP 1 UNTIL designsize DO
BEGIN	FOR j:= 1 STEP 1 UNTIL termcounter DO
resrow [j]:= designdata [j] . lastnumber;
weightroot := data [i, weightcounter];
observation:= data [i, designcounter] / weightroot;
fittedvalue:= vecvec (1, termcounter, resrow, depvar) / weightroot;
FOR j:= 1 STEP 1 UNTIL termcounter DO
resdev [j]:= tamvec (1, termcounter, j, data, resrow);
standdev:= Sqrt (vecvec (1, termcounter, resrow, resdev) * msres) / weightroot;
residual:= observation - fittedvalue;
residualsum:= residualsum + residual * weightroot * weightroot;
standres:= IF perfectfit THEN 0 ELSE residual / sqrtmsres;
studdev:= msres - standdev * standdev;
studres:= IF studdev LE precision * msres THEN 0 ELSE residual / Sqrt (studdev);
IF studres * studres GT maxstud THEN
BEGIN studnum:= i; maxstud:= studres * studres; END;
Outint (i, 6); Outfix (observation, 6, 20); Outfix (fittedvalue, 6, 20);
Outfix (standdev, 6, 20); Outfix (residual, 6, 20);
Outfix (standres, 6, 20); Outfix (studres, 6, 20); Outimage;
IF Mod (i, 10) EQ 0 AND i LT designsize THEN
BEGIN IF Mod (i, 50) EQ 0 THEN newpage2 ELSE Eject (Line + 1); END;
INSPECT outputstream DO
BEGIN	Outreal (observation, 16, 22); Outreal (fittedvalue, 16, 22);
Outreal (standdev, 16, 22); Outreal (residual, 16, 22);
Outreal (standres, 16, 22); Outreal (studres, 16, 22); Outimage;
END;
END;
dfstud:= designsize - termcounter - 1; frstud:= dfstud + 1 - maxstud;
studprob:= IF dfstud EQ 0 THEN 1 ELSE IF frstud LE precision THEN 0
ELSE designsize * (1 - Fisher (maxstud * dfstud / frstud, 1, dfstud));
Eject (Line + 1); Setpos (56);
Outtext ("sum of residuals :"); Outfix (residualsum, 6, 13); Outimage;
Eject (Line + 1); Setpos (20 - field (studnum));
Outtext ("Upper bound for the right tail probability of ");
Outtext ("the largest absolute studentized residual (no. ");
Outint (studnum, field (studnum)); Outtext (") :"); Setpos (118);
Outfix (Lmin (studprob, 1&&0), 6, 9); Outimage;
END *** residual_analysis ***;

COMMENT ************************************* process submodels **********************************;

INTEGER subselection, submodel, notprocessed;

PROCEDURE resetdesigndata;
BEGIN	INTEGER j;
FOR j:= 1 STEP 1 UNTIL termcounter DO designdata [j] . reset;
END *** resetdesigndata ***;

PROCEDURE process_submodels;
BEGIN	IF optionwanted [5] THEN
BEGIN	IF subselection EQ 0 THEN
BEGIN	submodel:= submodel + 1;
termcounter:= termcounter - 1;
IF termcounter GE minterm THEN
BEGIN resetdesigndata; GOTO process; END;
notprocessed:= minterm - 1;
END ELSE
FOR programcounter:= programcounter + 1
WHILE programcounter LE optioninstructions DO
BEGIN	submodel:= program [programcounter];
termcounter:= designcounter - 1 - submodel;
IF termcounter GE minterm THEN
BEGIN resetdesigndata; GOTO process; END;
notprocessed:= notprocessed + 1;
END;
END;
IF notprocessed GT 0 THEN INSPECT printstream DO
BEGIN	Eject (Line + 5); Outtext ("Transformed submodels with as only independent");
Outtext (" variable a constant, are not processed    (message)"); Outimage;
Eject (Line + 5); Outtext ("Number of submodels not processed :");
Outint (notprocessed, 3); Outimage;
END;
END *** process_submodels ***;

COMMENT ************************************* regression analysis ********************************;

PROCEDURE initialize_regression;
BEGIN	INTEGER i, j;
FOR i:= 1 STEP 1 UNTIL designsize DO
BEGIN	FOR j:= IF submodel EQ 0 THEN weightcounter ELSE termcounter
STEP - 1 UNTIL 1 DO data [i,j]:= designdata [j] . lastnumber;
depvar [i]:= data [i, designcounter];
END;
IF submodel EQ 0 THEN
FOR j:= 1 STEP 1 UNTIL designcounter DO
resprod [j]:= tamvec (1, designsize, j, data, depvar);
END *** initialize_regression ***;

resetindex:= designdata [1] . index;
matrixwidth:= (imagelength - 21) // 12;
subselection:= optioninstructions - inputinstructions;
programcounter:= inputinstructions; submodel:= notprocessed:= 0;

process : initialize_regression;	transformed_data_matrix;	control_information;
correlation_matrix;		no_regression_analysis;		least_squares;
multiple_correlation;		explained_variation;		standard_error_deviation;
regression_parameters;	covariance_matrix;		analysis_of_variance;
null_hypotheses;		residual_analysis;		process_submodels;

END *** regression ***;

END *** designdata ***;

END *** translator ***;

COMMENT ************************************* final messages *************************************;

endjob : INSPECT outputstream DO
BEGIN Outtext ("""Eor"" of job :"); Outint (jobnumber, 3); Outimage; END;

IF printstream =/= Sysout THEN
BEGIN	IF erroneous THEN BEGIN Outtext ("Erroneous"); Outimage; END;
Outtext ("End of job :"); Outint (jobnumber, 3); Outimage;
END;

INSPECT printstream DO
BEGIN	Eject (Line + 4); Outimage; Dotypeout (Sysout);
Outtext ("End of job :"); Outint (jobnumber, 3); Outimage;
IF printstream =/= Sysout THEN Eject (1); erroneous:= FALSE; GOTO job;
END;

exit : Eject (Line + 1); Dotypeout (Sysout);
Outtext ("End of Multiple Linear Regression Analysis"); Outimage; closefiles;

END *** environment ***;

END *** inputsystem ***;
```