Trailing-Edge
-
PDP-10 Archives
-
decuslib20-05
-
decus/20-0149/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 ***;
PROCEDURE readinpspec (filespec); TEXT filespec;
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;
END *** readinpspec ***;
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;
END *** readprispec ***;
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;
END *** readdatspec ***;
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));
END *** readoutspec ***;
PROCEDURE readfilespecs;
BEGIN runable:= TRUE; IF lastchar EQ '/' AND lastline . More THEN
BEGIN TEXT scanline;
lastline:- Copy (Rest (lastline)); scanline:- Scanto (lastline, '=');
readprispec (Scanto (scanline, ';')); readoutspec (Scanto (scanline, ';'));
readinpspec (Scanto (lastline, ';')); readdatspec (Scanto (lastline, ';'));
IF NOT runable THEN readline (inputstream);
END;
END *** readfilespecs ***;
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 ***;
Link CLASS databuffer;
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;
pointer:- NONE; base:- NEW Head;
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;
BEGIN INTEGER address;
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 *******************************;
BOOLEAN endtext, endline, dataread, dataskipped, datafileread;
INTEGER datarecord; TEXT lastline, lastitem; CHARACTER lastchar;
REF (textstorage) currenttext, modeltext, inputtext, optiontext; REF (inputstorage) inputdata;
PROCEDURE readtext (sourcetext);
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;
END *** readtext ***;
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;
END *** readline ***;
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;
END; readline (stream);
END; dataskipped:= erroneous;
dataread:= size GT 0 AND NOT dataskipped;
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;
END *** readdata ***;
PROCEDURE readdatafile;
INSPECT datastream DO
BEGIN readline (datastream);
WHILE (endtext OR endline) AND NOT Endfile DO
BEGIN IF endtext THEN
BEGIN datarecord:= datarecord + 1; endtext:= FALSE;
IF readkeyword EQ "EO" THEN readline (datastream);
END; IF endline THEN readline (datastream);
END; readdata (datastream); datarecord:= datarecord + 1;
END *** readdatafile ***;
TEXT PROCEDURE readkeyword;
BEGIN readkeyword:- Upcase (Copy (nextitem (2)));
WHILE nextchar NE '"' AND NOT endline DO;
WHILE nextchar LT ' ' AND NOT endline DO;
END *** readkeyword ***;
PROCEDURE readjob;
BEGIN TEXT keyword; readline (inputstream);
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;
END *** readjob ***;
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;
readline (inputstream);
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;
dataread:= dataskipped:= datafileread:= savemodel:= FALSE;
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];
PROCEDURE load (instruction, address); INTEGER instruction, address;
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
BEGIN load (tst, 0); load (fls, 0); lefthandpart; load (dep, 0);
IF lastatom == equal THEN nextatom ELSE error (11);
IF lastatom == addsign THEN nextatom; righthandpart; load (rtn, 1);
END OTHERWISE error (10);
PROCEDURE lefthandpart;
BEGIN idtype:= 2; expression; load (rep, 0); load (rtn, 2);
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;
WHILE lastatom == addsign DO
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;
IF lastatom == addsign THEN
BEGIN load (stk, 0); nextatom; term; load (add, 0);
nextterm;
END ELSE
IF lastatom == subsign THEN
BEGIN load (stk, 0); nextatom; term; load (sub, 0);
nextterm;
END *** nextterm ***;
PROCEDURE term; BEGIN factor; nextfactor; END;
PROCEDURE nextfactor;
IF lastatom == mulsign THEN
BEGIN load (stk, 0); nextatom; factor; load (mul, 0);
nextfactor;
END ELSE
IF lastatom == divsign THEN
BEGIN load (stk, 0); nextatom; factor; load (div, 0);
nextfactor;
END ELSE
IF lastatom == idisign THEN
BEGIN load (stk, 0); nextatom; factor; load (idi, 0);
nextfactor;
END *** nextfactor ***;
PROCEDURE factor; BEGIN primary; nextprimary; END;
PROCEDURE nextprimary;
WHILE lastatom == ttpsign DO
BEGIN load (stk, 0); nextatom; primary; load (ttp, 0); END;
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
BEGIN INTEGER address; address:= programcounter;
load (0, 0); load (stk, 0); nextatom; description (TRUE);
program [address]:= programcounter * 100 + jmp;
END ELSE load (eql, 0);
PROCEDURE description (repeat); BOOLEAN repeat;
BEGIN INTEGER address; address:= programcounter;
IF lastatom == lparent THEN
BEGIN nextatom; inputstatement;
IF lastatom == rparent THEN nextatom ELSE error (23);
IF repeat THEN load (rtn, address);
END ELSE
IF lastatom == lbracket THEN
BEGIN three:= four:= five:= FALSE;
nextatom; variablelist;
IF lastatom == rbracket THEN nextatom ELSE error (24);
IF repeat THEN mixed (address)
ELSE IF five THEN load (fix, address);
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);
PROCEDURE mixed (address); INTEGER address;
IF NOT (three OR five) THEN skip (address) ELSE
IF NOT (three OR four) THEN BEGIN load (rtn, address); load (fix, address); END
ELSE BEGIN IF five THEN load (fix, address); load (rtn, address); END;
PROCEDURE skip (address); INTEGER address;
BEGIN INTEGER count; count:= programcounter - address;
programcounter:= address - 1; load (skp, count);
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
BEGIN leftcounter:= programcounter:= 1; load (tcv, 1); load (ini, 0); END;
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;
designdata [term] . address:= address;
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;
rightdata [right] . address:= address;
END ELSE
IF type EQ 5 AND stack [address] EQ 0 THEN
BEGIN left:= left + 1; stack [address]:= left;
leftdata [left] . address:= address;
END ELSE
IF type EQ 2 THEN error1 (43, term)
WHEN number DO IF entry . Pos EQ 1 THEN
BEGIN stack [address]:= Scanreal (entry);
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 OR datafileread) THEN readdatafile;
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;
tcv: F:= address; 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;
GOTO function [address];
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 ***;
PROCEDURE heading;
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);
END *** heading ***;
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);
Outtext ("obs.no."); heading;
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;
IF adjusted GT precision THEN
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 ***;