Programmer 7500/MAX_PROGRAMMERS.iso/PASCAL/HPSTAT.ZIP/HPSTAT2.PAS

Pascal/Delphi Source File | 1986-04-15 | 57.0 KB | 2,249 lines

Program Statistics; {$C-}
Const
  DC2           = ^O;
Type
  Panel         = Array[0..24,0..79] of
                  record
                    ch : char; attr : byte;
                  end;
  Address       = ^Integer;

  DataPointer   =  ^DataRecord;
  DataRecord    =  record
                     XX : Array[1..10900] of Real;
                   End;
  Regs          =  record
                     AX,BX,CX,DX,BP,SI,DI,DS,ES,Flags: Integer;
                   End;
  SmallRow      =  array[1..31] of Real;
  DataRow       =  array[1..100] of Real;
  SmallArray    =  array[1..31,1..31] of Real;
  IntRow        =  array[1..31] of Integer;
  Message       =  string[255];
  A_Line                    =     string[68];

Var
  graphics_scrn    : Panel absolute $B800:$0000;
  screen           : ^Panel;
  SaveScreen       : panel;
  X_save,Y_save    : integer;

  LabXY                        : array[1..2,1..20] of integer;
  Titles                       : Array[1..7] of A_Line;
  PageBuf                 : Array[1..8] of Integer;
  NDP                     : Array[1..100] of Integer;
  HeapTop                 : ^Integer;
  RR                      : SmallArray;
  FirstData,LastData,
  Data                    : DataPointer;
  CH,Comma,Iff                         : Char;
  Avg,StDev,SS,Diagonal,Vertical       : DataRow;
  T,Berr,Bmin,Bmax,SX,A,RY,B           : SmallRow;


  Mst,Sdr,FX,SY,XDF,TTT,Xnobs,RSUM,SSUM,
  Z,YY,SSB,SSB0,SST,RC1,RC2,S2R,S1R,RCP,
  SSR,Yrange,PSDR,PSDR1,Bmod,U,V,W,X,Y            : Real;

  I,J,K,L,M,N,Xstart,Ystart,Code,MR,MC,Skips,
  IY1,Imod,ICX,IIC,Nobs,Nvar,DF1,DF,FC,Irow,
  IY,Xpos,Ypos,LastX,LastY,Degree,RC,SC,PlotCall  : Integer;


  Result                  : Regs;
  Ansr,Cmd                : Char;
  FileName                : String[20];
  InFile,OutFile,IO       : Text;
  Amount,Stn              : String[10];
  Line                    : String[80];
  IOErr,FileWrite         : Boolean;
  Icol                    : IntRow;
  Len,Decimal             : Array[1..10] of Integer;

Function Sign(A: Real): Real;
Begin
  If A > 0 then Sign := 1.0
  Else Sign := -1.0;
End;

Function Yes: Boolean;
Var
  Ch  : Char;
Begin
  Repeat
    Write(' (Y,N)?');
    Read(Kbd,Ch);
    Writeln;
  Until Ch in['N','n','Y','y'];
  Yes := Ch in['Y','y']
End;

Procedure Pause;
Begin
  Writeln;
  GotoXY(1,24); TextBackground(4);
  Write(' Press Any Key to Continue ');
  TextBackground(0);
  Read(Kbd,Ch);
  GotoXY(1,24); TextBackground(0);
  Write('                           ');
End;

Function Power(X,Y: Real): Real;
Const
  Epsilon    = 0.000001;

Var
  I: Integer;
  P: Real;

Begin
  If Abs(Y - Round(Y)) < Epsilon then
  Begin
    P := 1.0;
    If Y >= 0.0 then
      For I := 1 to Round(Y) do P := P*X
    Else if X = 0.0 then Writeln('Negative Power of 0.0')
    Else
      For I := -1 downto Round(Y) do
        P := P/X;
    Power := P;
  End
  Else if X > 0.0 then
    Power := Exp(Y*Ln(X))
  Else if Abs(X) < Epsilon then
    Power := 0.0
  Else
    Writeln('Attempt to take negative number to non-integer power')
End;

Procedure IOcheck;
Const
 IoVal   : Integer = 0;
var
  Ch                   : Char;
begin
  IOVal := IOresult;
  IOErr := (IOVal <> 0);
  GotoXY(1,23); ClrEol;        { Clear error line in any case }
  if IOErr then begin
    Write(Chr(7));
    case IOVal of
      $01  :  Writeln('File does not exist');
      $05  :  Writeln('Can''t read from this file');
      $06  :  Writeln('Can''t Write to this file');
      $10  :  Writeln('Error in numeric format');
      $91  :  Writeln('Seek beyond end of file');
      $99  :  Writeln('Unexpected end of file');
      $F0  :  Writeln('Disk Write error');
      $F1  :  Writeln('Directory is full');
      $F2  :  Writeln('File size overflow');
      $FF  :  Writeln('File disappeared')
    else      Writeln('Unknown I/O error:  ');
    end;
    Pause;
  end
end; { of proc IOCheck }

Function  GetX(I,J: Integer): Real;
  Var Index : Integer;
Begin
  Index := (J-1)*MR + I;
  GetX := Data^.XX[Index];
End;

Procedure PutX(Number: Real; I,J: Integer);
  Var Index : Integer;
Begin
  Index := (J-1)*MR + I;
  Data^.XX[Index] := Number;
End;

Procedure InitializeArray;
Begin
  For M := 1 to 10900 do
    Data^.XX[M] := -2.0E35;
End;

Procedure CheckSkips;
Begin
  ClrScr;
  Skips := 0;
  Write(^J' Do you wish to skip any observations for this analysis');
  If Yes then
    Begin
      Write(^J' How many observations do you wish to skip ? ');
      Readln(Skips);
      For J := 1 to Skips do
        Begin
          GotoXY(1,20);
          Write(' Enter an observation number to skip : '); ClrEol;
          Readln(Irow);
          PutX(1.0,Irow,SC);
        End;
    End;
End;

Procedure UnsKip;
Begin
  For J := 1 to Nobs do
    Begin
      If GetX(J,SC) > 0.0 then
        PutX(-2.0E35,J,SC);
    End;
End;

Procedure MatrixMultiply(  C: SmallRow;
                           B: SmallArray;
                           N: Integer;
                       Var A: SmallRow);
{ THIS SUBROUTINE MULTIPLIES TWO MATRICES  }
Var
  Temp   : Real;
Begin
  For I := 1 to N do
    Begin
      Temp := 0.00;
      For J:= 1 to N do
        Temp := Temp + B[I,J]*C[J];
      A[I] := Temp;
    End;
End;

Procedure ChooseColumns(Var IY,ICX: Integer;
                        Var Icol  : IntRow ; MinY: Integer);
Var
  OldLen          : Array[1..31] of Integer;
  Col             : Array[1..3] of Integer;
  X,Y,Index       : Integer;
  Amount          : String[5];

Procedure Initialize;
Begin
  For I := 1 to 31 do Icol[I] := 0 ;
  Amount := '';
End;


Procedure Convert(Var Len: Integer);
Var
  Value              :Integer;
Begin;
  Len := Length(Amount);
  Val(AMOUNT,Value,Code);
    Icol[Index] := Value
End;

Procedure WriteNum(A,B,C: Integer);
Begin
  GotoXY(A,B);
  If OldLen[Index] > C then
    Begin
      For L := 1 to OldLen[Index] do
        Write(' ');
      GotoXY(A,B);
    End;
  If Icol[Index] <> 0 then Write(Icol[Index]);
  OldLen[Index] := C;
End;


Procedure CheckNumber;
Var
  C  : Integer;
Begin
  C := 10;
  If(Amount <> '') then
    Convert(C);
  GotoXY(X,Y);
  WriteNum(X,Y,C);
  IY := 0;
  Amount := '';
End;

Begin
  Col[1] := 15; Col[2] := 40; Col[3] := 65;
  For J := 1 to 31 do OldLen[J] := 0;
  Initialize;
{
       READ INPUT DATA
}
  ClrScr;
  WriteLn(' Enter Column Numbers for the Regression Variables below. Use');
  Writeln(' the GotoXY Keys to move arount on the entry field. For the ');
  Writeln(' Correlation Matrix, NO Y VARIABLE IS NEEDED ');
  TextBackground(15); TextColor(0);
  Writeln(' TYPE X WHEN FINISHED ');
  TextBackground(0);  TextColor(15);
  If MinY = 1 then
    Begin
      GotoXY(1,6);Write(' Y - Axis :');
    End;
  For I := 1 to 3 do
    Begin
      For J := 1 to 10 do
        Begin
          K := (I-1)*25 + 1;
          GotoXY(K,8+J);
          K := (I-1)*10 + J ;
          Write(' X - ',K:2,' : ');
        End;
    End;
  I := 1; J := MinY;
  GotoXY(Col[I],6+(J-1)*3); Index := 1; IY := 0;
  Repeat
  Read(Kbd,CH);
  If(CH = ^[) then
    Begin
      Read(Kbd,Ch);
      CheckNumber;
      Case CH of
      'P':  Begin
              J := J + 1;
              If (I = 1) and (J > 11) then
                J := MinY
              Else If (I > 1) and  (J > 10) then
                J := 1  ;
            End;
      'H':  Begin
              J := J - 1;
              If(J < 1) and (I > 1) then J := 1;
              If(J < MinY) and (I = 1) then J := MinY;
            End;
      'G':  Begin
              I := 1;
              J := MinY;
            End;
      'M':  Begin
              I := I + 1;
              If (I > 3) then I := 3 ;
              If (J > 1) and (I = 2) then J := J - 1;
            End;
      'K':  Begin
              I := I - 1;
              If (I < 1) then I := 1;
              If (I = 1) and (MinY = 2) then J := J + 1;
            End;
      End;
    End;
    X := Col[I];
    Y := 8 + J;
    Index := (I-1)*10 + J;
    If I = 1 then
      Begin
        Index := Index - 1;
        Y := Y - 1;
        If J = 1 then Y := 6;
      End;
    If Index = 0 then Index := 31;
    GotoXY(X,Y);
  If(CH in ['E','e','-','.','0'..'9']) then
    Begin
      IY:= IY + 1;
      If (IY > 5) then IY := 5;
      Insert(Ch,AMOUNT,IY);
      GotoXY(X,Y); Write(Amount);GotoXY(X+IY,Y);
    End;
  If (CH = ^H) then
    Begin
      Delete(Amount,IY,1);
      GotoXY(X,Y);
      Write(Amount,' ');
      GotoXY(X+IY-1,Y);
      IY := IY - 1;
      If IY < 0 then IY := 0;
    End;
   If (Ch = 'X') or (Ch = 'x') then
     CheckNumber;
  Until (CH = 'X') or (CH = 'x');
  ClrScr; TextBackground(0);
  J := 1; ICX := 0 ;
  While (Icol[J] > 0) and (J <= 30) do
    Begin
      ICX := ICX + 1;
      J := J + 1;
    End;
  If MinY = 1 then
    Begin
      Icol[Icx+1] := Icol[31];
      IY := Icol[Icx+1]
    End;
End;



Procedure MatrixInvert(Var A:SmallArray; N : Integer );
Var
  KROW,IROW,JJ,KP1  : Integer;
  Inter             : Array[0..50,0..2] of Integer;
{
      THIS SUBROUTINE CALCULATES THE INVERSE OF A MATRIX USING
      PARTIAL PIVOTING, WHICH AVOIDS PROBLEMS WITH VERY SMALL
      PIVOT POINTS
 }

Begin
  For K := 1 to N do
    Begin
{
      SEARCH FOR LARGEST PIVOT ELEMENT
}
      JJ := K;
      IF(K <> N) Then
        Begin
          KP1 := K+1;
          Z := ABS(A[K,K]);
          For I := KP1 to N do
            Begin
              V := ABS(A[I,K]);
              IF(Z < V) Then
                Begin
                  Z := V;
                  JJ := I;
                End;
            End;
        End;
{
      STORE NUMBERS OF ROWS INTERCHANGED
}
      INTER[K,1] := K ;
      INTER[K,2] := JJ;
      IF(JJ <> K) Then
        Begin
{
      ROW INTERCHANGE
}
          For J :=  1 to N do
            Begin
              W := A[JJ,J];
              A[JJ,J] := A[K,J];
              A[K,J] := W;
            End;
        End;
{
      CALCULATE NEW ELEMENTS OF PIVOT ROW EXCEPT FOR PIVOT ELEMENT
}
      For J := 1 to N do
        Begin
          IF(J <> K) Then
             A[K,J] := A[K,J]/A[K,K];
        End;
{
      CALCULATE NEW ELEMENT REPLACING PIVOT ELEMENT
}
      A[K,K] := 1.0/A[K,K];
{
      CALCULATE NEW ELEMENTS NOT IN PIVOT ROW OR COLUMN
}
      For L := 1 to N do
        Begin
          IF(L <> K) Then
            Begin
              For M := 1 to N do
                Begin
                  IF(M <> K) Then
                  A[L,M] := A[L,M]-A[K,M]*A[L,K];
                End;
            End;
        End;
{
      CALCULATE NEW ELEMENTS FOR PIVOT COLUMN, EXCEPT FOR PIVOT ELEMENT
}
      For I := 1 to N do
        Begin
          IF(I <> K) Then
            A[I,K] := -A[I,K]*A[K,K];
        End;
    End;
{
      REARRANGE COLUMNS OF FINAL MATRIX
}
    For L := 1 to N do
      Begin
        K := N - L+1;
        KROW := INTER[K,1];
        IROW := INTER[K,2];
        IF(KROW <> IROW) Then
          Begin
            For I := 1 to N do
              Begin
                W := A[I,KROW];
                A[I,KROW] := A[I,IROW];
                A[I,IROW] := W;
              End;
          End;
      End;
End;


Procedure StandardDeviation(A,B: Integer;Var StDev: Datarow);
{
      CALCULATES STANDARD DEVIATIONS FOR ALL VARIABLES IN AN ARRAY
}
Var
  V,S1,S2          : Real;
Begin
  Z := Int(Nobs - Skips);
  For K := 1 to B do
    Begin
      S1 := 0.0; S2 := 0.0;
      For J := 1 to A do
        Begin
          If(GetX(J,SC) < 0.0) then
            Begin
              S1 := S1 + GetX(J,K);
              S2 := S2 + GetX(J,K)*GetX(J,K);
            End;
        End;
      V  := (S2-(S1*S1)/Z)/(Z-1.0);
      StDev[K] := Sqrt(V);
    End;
End;


Procedure Average(Var Avg: DataRow; A,B: Integer);
{
      CALCULATES THE AVERAGES FOR ALL VARIABLES IN AN ARRAY
}
Begin
  Z := Int(Nobs - Skips);
  For I := 1 to B do
    Begin
      Avg[I]:= 0.0;
      For J := 1 to A do
        Begin
          If GetX(J,SC) < 0.0 then
            Avg[I] := Avg[I] + GetX(J,I)
        End;
      Avg[I] := Avg[I]/Z;
    End;
End;


Procedure Scorr(Rows: IntRow; Nobs,Ind : Integer;
                Var RR: SmallArray);
{
      CALCULATES THE CROSS CORRELATION COEFFICIENTS FOR SELECTED VARIABLES
}
Var
  D1,D2,X1,X2,X3  : Real;
  I,J,K           : Integer;
Begin
  Z := Int(Nobs - Skips);
  ClrScr; GotoXY(10,12);
  TextBackground(4);
  Writeln(' WAIT A SECOND, I''M THINKING ');
  TextBackground(0);
  FillChar(RR,5400,Chr(0));
  For I := 1 to Nobs do
    Begin
      If GetX(I,SC) < 0.0 then
        Begin
          For J := 1 to Ind do
            Begin
              D1 := GetX(I,Rows[J]) - Avg[Rows[J]];
              For K := 1 to J do
                Begin
                  D2 := GetX(I,Rows[K])-AVG[Rows[K]];
                  RR[J,K] := RR[J,K] + D1*D2;
                  If J <> K then RR[K,J] := RR[J,K];
                End;
            End;
        End;
    End;
  For J := 1 to Ind do
    Begin
      Diagonal[J] := RR[J,J];
      Vertical[J] := RR[J,Ind];
    End;
  For J := 1 to Ind do
    Begin
      For K := 1 to J do
        Begin
          X2 := RR[J,K]/Z;
          X1 := SQRT(Diagonal[J]/Z);
          X3 := SQRT(Diagonal[K]/Z);
          If X1*X3 <> 0.0 then
            RR[J,K] := X2/(X1*X3)
          Else
            RR[J,K] := 1.0;
          If(K = J) Then RR[J,K] := 1.0000;
          If K <> J then RR[K,J] := RR[J,K];
        End;
    End;
    ClrScr;
End;

Procedure InitCorr;
Begin
  Average(Avg,Nobs,Nvar);
  StandardDeviation(Nobs,Nvar,StDev);
End;

Procedure DifSave;

{                                                        }
{     PROGRAM TO CONVERT RPLOT DATA FILES TO DIF FORMAT  }
{     FOR USE IN LOTUS, VISICALC, ETC.                   }
{                                                        }
Begin
  Repeat
   Writeln(^J' ENTER A FILENAME. A  ''.DIF'' SUFFIX WILL AUTOMATICALLY ');
   Write(' BE APPENDED TO THE FILE NAME...');
   Readln(FileName);
   If Pos('.',FileName) = 0 then
     FileName := FileName + '.DIF';
{$I-} Assign(OutFile,FileName);
     Rewrite(OutFile); IOCheck; {$I+}
  Until Not IOerr;
  Writeln(OutFile,'TABLE');
  Writeln(OutFile,'0,1');
  Writeln(OutFile,'""');
  Writeln(OutFile,'VECTORS');
  Writeln(OutFile,'0,',Nvar);
  Writeln(OutFile,'""');
  Writeln(OutFile,'TUPLES');
  Writeln(OutFile,'0,',Nobs);
  Writeln(OutFile,'""');
  Writeln(OutFile,'DATA');
  Writeln(OutFile,'0,0');
  Writeln(OutFile,'""');
  Writeln(OutFile,'-1,0');
    For I := 1 to Nobs do
    Begin
      Writeln(OutFile,'BOT');
      For J := 1 to Nvar do
      Begin
        Writeln(OutFile,'0,',GetX(I,J));
        Writeln(OutFile,'V');
      End;
      Writeln(OutFile,'-1,0');
    End;
  Writeln(OutFile,'EOD');
  Close(OutFile);
End;

Procedure SetPageSize;
Begin
 If Nvar*Nobs > 10900.0 then
   Begin
     ClrScr;GotoXY(1,12);
     TextBackground(4);
     Write(' TOO MANY DATA POINTS!! ');
     TextBackground(0);
   End
 Else
   Begin
     Release(HeapTop);
     Mark(HeapTop); New(Data);
     MC := Nvar + 8;
     MR := 10900 div MC;
     GotoXY(1,5);
     Writeln(' SUPERSTAT  has sized  the data area to ',MC-3,' columns  ');
     Write(' and ',MR,' rows.  Do you wish to change this');
     If Yes then
       Begin
         Writeln(' To change this configuration, enter the number of ');
         Writeln(' Columns you wish the data area to have ');
         Writeln(' (Remember that Columns*Rows can not be greater than 10000) ');
         Write(' Enter the number of Columns '); Readln(MC);
         MC := MC + 3; MR := 10900 div MC;
       End;
     FC := MC -1; RC := MC ; SC := MC - 2;
     InitializeArray;
   End;
End;


Procedure DFread ;

{ PROCEDURE TO READ DIF FILE }

Const
  Vectors  = 'VEC';
  Tuples   = 'TUP';
  Bot      = 'BOT';
Var
   KK,S,N                   :     Integer;
   Alpha                    :     string[3];
   V                        :     real;
   NS                       :     string[25];
   Title                    :     string[72];
   AA                       :     string[2];
   SS                       :     string[1];
begin;
ClrScr;
  for I := 1 to 4 do
  begin
    Readln(Infile,Alpha);
    Readln(Infile,SS,Comma,N);
    Readln(Infile,Comma);
    If(Alpha = Vectors) then  Nvar := N;
    If(Alpha = Tuples) then   Nobs := N;
  end;
 SetPageSize;
 KK := 0;
 For I := 1 to Nobs do
  Begin
    Readln(InFile,AA);
    If(AA <> '-1') then Halt;
    Readln(InFile,Alpha);
    If(Alpha <> Bot) then Halt;
      For J := 1 to Nvar do
      Begin
       Readln(InFile,SS,Comma,V);
       Readln(InFile,Title);
       Val(SS,S,Code);
       If(S = 0) then  PutX(V,I,J)
       else
       KK := KK + 1;
      End;
  End;
Close(InFile)
End;

Procedure PRN_READ ;
begin;
  ClrScr;
  Write(^J^J' ENTER THE NUMBER OF VARIABLES IN THE PRN FILE...');
  Readln(Nvar);
  Write(^J' ENTER THE NUMBER OF OBSERVATIONS PER VARIABLE ...');
  Readln(Nobs);
  SetPageSize;
  For I := 1 to Nobs do
    Begin
      For J := 1 to Nvar do
        Begin
          Read(InFile,V);
          PutX(V,I,J);
        End;
    End;
  Close(InFile)
End;

Procedure SpreadSheet;

Procedure ConvertToReal;
  Var
    Len, NDX              :Integer;
  Begin;
    Len := 0; NDX := 0;
    Repeat
      NDX := NDX + 1;
      If(COPY(AMOUNT,NDX,1) = ' ') then  Len := NDX - 1
    Until (Len > 0);
    Val(Copy(AMOUNT,1,Len),W,Code);
    PutX(W,I,J);
   End;

Procedure WriteData(Row,Col : Integer);
  Begin
    If(GetX(Row,Col) > -1.0E35) then
       Write(GetX(Row,Col):10:NDP[Col])
    Else Write('          ');
  End;

Procedure HiLite;
  Begin
    TextBackground(15);TextColor(0); WriteData(I,J);
    TextBackground(0); TextColor(15);
  End;

Procedure Xrow(Start : Integer);
 Var NDX : Integer;
 Begin
   For M := 0 to 5 do
     Begin
       NDX := M*12 + 13;
       GotoXY(NDX,4);
       Write(Start+M:2);
     End;
 End;

Procedure Yrow(Start : Integer);
Begin
  For M := 0 to 14 do
    Begin
      GotoXY(1,M+6);
      Write(Start+M:4);
    End;
End;

Procedure UnLite;
  Begin
    GotoXY(XPOS,YPOS);
    TextBackground(0) ;WriteData(I,J);
  End;


Procedure EntString;
Begin
  GotoXY(1,22);Write('Value:  ');ClrEol;
End;

Procedure Where(A,B: Integer);
Begin
  GotoXY(1,3);
  Write('Cell[',A:3,',',B:2,']');
End;


Procedure CheckAmount;
Begin
  If(Copy(Amount,1,10) <> '          ') then
    ConvertToReal;
  GotoXY(XPOS,YPOS);
  WriteData(I,J);
  IY := 0;
  Amount := '          ';
End;

Procedure CountData;
Begin
  Nobs := 1; Nvar := 1;
  While GetX(Nobs,1) > -1.0E35 do
    Nobs := Nobs + 1;
  While GetX(1,Nvar) > -1.0E35 do
    Nvar := Nvar + 1;
  Nobs := Nobs - 1; Nvar := Nvar - 1;
End;

Procedure ClearPage;
Var
  Clear : Integer;
Begin
  ClrScr;
End;

Procedure FillScreen(NewScreen: Panel);
Var
   XX                   : Real;
   Str3                 : String[3];
   I,J                  : Integer;
Begin
   Window(1,1,80,25);
   For I  := 0 to 14 do
     Begin
       Str(I+Ystart:3,Str3);
       For J := 1 to 3 do NewScreen[I+5][J].ch := Str3[J];
       For K := 0 to 5 do
         Begin
           L  := K*12 + 5;
           M  := Xstart + K;
           XX  := GetX(Ystart+I,M);
           If(XX > -1.0E35) then
             Str(XX:10:NDP[M],Stn)
           Else Stn := '          ';
           For J := 1 to 10 do NewScreen[I+5][L+J].ch := Stn[J];
         End;
     End;
     Screen^ := NewScreen;
End;

Procedure Move;
  Begin
    If Xstart = LastX then
      Begin
        If(J < Xstart) then Xstart := J;
        If(J > Xstart + 5) then Xstart := J-5;
        If (Xstart < 1) then Xstart := 1;
      End;
    If Xstart <> LastX then
      Begin
        SaveScreen := Screen^
        FillScreen(SaveScreen);
        Xrow(Xstart);
      End;
    If Ystart = LastY then
      Begin
        If(I > (Ystart+14)) then Ystart := I - 14;
        If(I < Ystart)  then Ystart := I ;
        If(I < 2) then Ystart := 1;
      End;
    If Ystart <> LastY then
      Begin
        SaveScreen := Screen^
        FillScreen(SaveScreen);
      End;
    XPOS := 7 + 12*(J - Xstart);
    YPOS := 6 + I - Ystart;
    Where(I,J);
    GotoXY(XPOS,YPOS);Hilite;EntString;
    LastX := Xstart; LastY := Ystart;
  End;


Procedure ScreenEntry;
Begin
  ClearPage;  Screen := addr(graphics_scrn);
  XPOS := 7;  YPOS := 6;  AMOUNT := '          ';
  Ystart := 1; Xstart := 1; I := 1; J:= 1;
  LastX := Xstart; LastY := Ystart;
  ClrScr;
  TextBackground(15); TextColor(0);
  Write('SUPERSTAT  Statistical Analysis Program ||');
  Writeln(' Type X when  finished ');
  TextBackground(0); TextColor(15);
  GotoXY(1,3); TextBackground(15); TextColor(0);
  Write('            ');
  TextBackground(0); TextColor(15);
  GotoXY(8,5);Write(Line);Write('==========');
  Move;
  Xrow(Xstart);Yrow(Ystart);GotoXY(Xpos,Ypos);HiLite;
  Where(I,J);
  SaveScreen := Screen^
  EntString;
  IY := 0;
  Repeat
  Read(Kbd,CH);
  If(CH = ^[) then
    Read(Kbd,CH);
  If Ch in['M','K','P','H',DC2,'Q','I','G',^I] then
    Begin
      CheckAmount;
      Unlite;
    End;
  Case CH of
    'P':  Begin
            I := I + 1;
            If(I > MR) then I := MR;
          End;
    'H':  Begin
            I  := I - 1;
            If(I < 1) then I := 1;
          End;
    'M':  Begin
            J := J + 1;
            If(J > MC-3) then J := MC-3;
            End;
    'K':  Begin
            J := J - 1;
            If(J < 1) then J := 1;
          End;


    DC2 :  Begin
            Xstart := Xstart - 6;
            J := J - 6;
            If J < 1 then J := 1;
            If Xstart < 1 then Xstart := 1;
          End;
    'Q':  Begin
            Ystart := Ystart + 14;
            I := I + 14;
            If Ystart > MR - 14 then Ystart := MR - 14;
            If I > MR then I := MR;
         End;
    'I':  Begin
            Ystart := Ystart  - 14;
            I := I - 14;
            If(Ystart < 1) then Ystart := 1;
            If(I < 1) then I := 1;
          End;
    'G':  Begin
            J := 1; I := 1;
          End;
     ^I:  Begin
            Xstart := Xstart + 6;
            If Xstart > MC - 9 then Xstart := MC - 9;
            J := Xstart;
          End;
     End;
  If Ch in['M','K','P','H',DC2,'Q','I','G',^I] then
       Move;
   If(CH in ['E','e','-','.','0'..'9']) then
     Begin
       IY:= IY + 1;
       If (IY > 10) then IY := 10;
       AMOUNT[IY] := CH;
       GotoXY(9,22); Write(Amount);ClrEol;GotoXY(9+IY,22);
     End;
   If( CH = ^M) then
     Begin
       CheckAmount;
       EntString;
     End;
   If (CH = ^H) then
     Begin
       Amount[IY] := ' ';
       GotoXY(9,22);Write(Amount);ClrEol;GotoXY(8+IY,22);IY := IY - 1;
       If (IY < 0) then IY := 0;
     End;
   If (CH = '/') then
     Begin
       GotoXY(1,22);ClrEol;
       Write(' Commands: (F)ormat, (S)ave ');
       Read(Kbd,Ch);Writeln;
       Case Ch of
       'F','f' : Begin
                 SaveScreen := Screen^
                 GotoXY(1,22);ClrEol;
                 Write('Enter the number of Decimal places for this column: ');
                 Readln(NDP[J]);FillScreen(SaveScreen);
                 End;
       'S','s' : Begin
                   CountData;
                   SaveScreen := Screen^
                   ClrScr;
                   GotoXY(1,15);
                   DifSave;
                   FillScreen(SaveScreen);
                 End;
        End;
        EntString;
      End;
  Until (CH = 'X') or (CH = 'x');
End;

Begin
  FillChar(Line,81,'=');
  For I := 1 to MC do NDP[I] := 4;
  FirstData := Nil;
  ScreenEntry;
  CountData;
  ClrScr;
  Writeln(' You entered ',Nvar,' variables and ',Nobs,' Observations');
  Pause;
End;

{Print out matrix of correlation coefficients }

Procedure CorrelationMatrix;
Begin
  ClrScr;
  GotoXY(1,2);
  Writeln(' This section prints out cross correlation coefficients ');
  Writeln(' for up to 30 variables at a time.  You must enter the  ');
  Writeln(' number of variables you wish to see in the Matrix, and ');
  Writeln(' then you may enter the column numbers of the variables ');
  Writeln(' in Question '^J);
  Pause;
  ChooseColumns(IY,ICX,Icol,2);
  Scorr(Icol,Nobs,ICX,RR);
  ClrScr;
  Writeln(IO,^J^J^J'Correlation Matrix':33);
  For J := 1 to ICX do
    Begin
      Write(IO,'     ');
      For K := 1 to J do
            Write(IO,Icol[K]:10); Writeln(IO);
      Write(IO,^J'  X[',Icol[J]:2,']');
      For K := 1 to J do
            Write(IO,RR[J,K]:10:4); Writeln(IO);
    End;
  Writeln(IO,^J^J^J);
  Pause;
End;


Procedure ColumnSum(Var A: Real; B,M: Integer);
Begin
  A := 0.0;
  For I := 1 to M do
    Begin
      If GetX(I,SC) < 0.0 then
        A := A + GetX(I,B)
    End;
End;

Procedure DotProduct(Var A: Real; B,C,M: Integer);
Begin
  A := 0.0;
  For I := 1 to M do
    Begin
      If GetX(I,SC) < 0.0 then
        A := A + GetX(I,B)*GetX(I,C)
    End;
End;

Function Minimum(A,B: Integer): Real;
 Var Min : Real;
Begin
  Min := 1.0E30;
  For K := 1 to  B do
    Begin
      If GetX(K,SC) < 0.0 then
         Begin
          If GetX(K,A) < Min then
            Min := GetX(K,A)
         End;
    End;
  Minimum := Min;
End;

Function Maximum(Var A,B: Integer): Real;
 Var Max : Real;
Begin
  Max := -1.0E30;
  For K := 1 to  B do
    Begin
      If GetX(K,SC) < 0.0 then
        Begin
         If GetX(K,A) > Max then
           Max := GetX(K,A);
        End;
    End;
  Maximum := Max;
End;

Procedure Analysis(L: Integer);
Var
  X_P                 : Real;
Begin
{
       CALCULATE PREDICTED Y VALUES AND RESIDUALS
}
  SY := 0.0;
  RSUM := 0.0;
  SSUM := 0.0;
  For I := 1 to L do
    Begin
      PutX(B[1],I,FC);
      SY := SY+(GetX(I,IY) - Avg[IY])*(GetX(I,IY) - Avg[IY]);
      For J := 2 to Degree+1 do
        Begin
          X_P := Int(J - 1);
          W := GetX(I,FC) + B[J]*Power(GetX(I,ICX),X_P);
          PutX(W,I,FC);
        End;
      PutX(GetX(I,IY) - GetX(I,FC),I,RC);
      RSUM :=RSUM+(GetX(I,FC)- Avg[IY])*(GetX(I,FC)-Avg[IY]);
      SSUM :=SSUM+(GetX(I,IY)- Avg[IY])*(GetX(I,IY)-Avg[IY]);
    End;
  RCP := RSUM/SSUM;
  RCP := RCP*100.00;
End;

{$I B:Simplex.Pas}

Overlay Procedure PolyRegr;
Var
  NP,NP2,KX,IJ,JK,Index    : Integer;
  PAVG,Sum,Ysum,PS3R,PS4R  : Real;

Procedure ColumnMultiply;
Begin
  For K := 1 to Nobs do
      PutX(GetX(K,ICX)*GetX(K,Nvar+1),K,Nvar+1);
End;

{
CHOLESKY'S METHOD FOR SOLUTION OF SIMULTANEOUS LINEAR ALGEBRAIC
EQUATIONS
}
Procedure Cholesky(A: SmallArray; M,N: Integer; Var C: SmallRow);
Var

M1,IP1,JM1,II,JJ  : Integer;
J,K,IM1,NN,IL     : Integer;
Sum               : Real;

{
CALCULATE FIRST ROW OF UPPER TRIANGULAR MATRIX
}
Begin;
  M1 := M+1;
  For J := 1 to M1 do
      A[1,J] := A[1,J]/A[1,1];
{
CALCULATE OTHER ELEMENTS OF U AND L MATRICES
}
  For IL := 2 to M do
    Begin
      J := IL;
      For II := J to M do
        Begin
          SUM := 0.0;
          JM1 := J-1;
          For K := 1 to JM1 do
            SUM := SUM+A[II,K]*A[K,J];
          A[II,J] := A[II,J]-SUM;
        End;
      IP1 :=IL + 1;
      For JJ := IP1 to M1 do
        Begin
          SUM := 0.0;
          IM1 := IL - 1;
          For K := 1 to IM1 do
            SUM := SUM+A[IL,K]*A[K,JJ];
          A[IL,JJ] := (A[IL,JJ]-SUM)/A[IL,IL];
        End;
    End;
  C[M] := A[M,M+1];
  L := M-1;
  For NN := 1 to L do
    Begin
      SUM := 0.0;
      IL := M-NN;
      IP1 := IL + 1;
      For J := IP1 to M do
      SUM := SUM+A[IL,J]*C[J];
      C[IL] := A[IL,M1]-SUM;
    End;
End;

Procedure PSXC(Var Z,Sum : Real);
Begin
  Sum := 0.0;
  For N := 1 to Nobs do
    Sum := Sum + Sqr(GetX(N,Nvar+1) - Z);
End;

Begin
  ClrScr;
  Writeln(' This section of SUPERSTAT determines polynomial fits');
  Writeln(' Up to a degree of 20.');
  Write(^J^J' Enter the degree of the polynomial for this problem: ');
  Readln(Degree);
  Write(' Which column contains the X data? '); Readln(ICX);
  Write(' Which column contains the Y data? '); Readln(IY);
  Ysum := 0.0;
  For J := 1 to Nobs do
    Begin
      PutX(1.0,J,Nvar+1);
      Ysum := Ysum + GetX(J,IY);
    End;
    SSB0 := Sqr(Ysum)/XNobs;
    NP := Degree+1;
{
      FILL IN CORRELATION MATRIX
}

    RR[1,1] := Nobs;
    RR[1,NP+1] := Ysum;
    NP2 := 2*Degree;
    For J := 1 to NP2 do
      Begin
        ColumnMultiply;
        If (J < NP + 1) then  DotProduct(Y,Nvar+1,IY,Nobs);
        ColumnSum(W,Nvar+1,Nobs);
        PAVG := W/Xnobs;
        IF(J < NP) Then  PSXC(PAVG,SX[J]);
        IJ := J+1 ;
        IF(IJ < NP) Then
          RR[IJ,1] := W
        Else
          RR[NP,IJ-Degree] := W;
          IF(J <= Degree) Then RR[J+1,NP+1] := Y;
        End;
      For K := 1 to Degree do
        For J := 1 to Degree do
          RR[J,K+1] := RR[J+1,K];
      KX := NP+1;
      For I := 1 to 21 do B[I] := 0.00;
{
      SOLVE THE CORRELATION MATRIX BY CHOLESKY'S METHOD
}
      For JK := 1 to NP do
       T[JK] := RR[JK,NP+1];
       Cholesky(RR,NP,KX,B);
       SSB := 0.0;
       For IJ := 1 to NP do
         SSB := SSB+B[IJ]*T[IJ];
       SST := SSB-SSB0;
       MST := SST/Degree;

      Analysis(Nobs);

      SSR :=SY-SST;
      S2R :=SSR/(NOBS-NP);
      IF(S2R <> 0.0) Then
        FX:=MST/S2R
      Else
        FX := 1.0E+08;
      S1R := 0.0010;
      IF(S2R > 0.0) Then S1R := SQRT(S2R);
      For J := 1 to Degree do
        Begin
          BERR[J] := S1R/SQRT(SX[J]);
          IF(BERR[J] = 0.0) THEN
            T[J] := 1.0E+08
          ELSE
            T[J] := B[J+1]/BERR[J];
        End;
      PS3R:=100.*S1R/Avg[IY];
      Yrange := Maximum(IY,Nobs) - Minimum(IY,Nobs);
      PS4R := 100.0*S1R/YRANGE;
      IF(RCP > 100.0) Then RCP := 100.0 ;
  ClrScr;Writeln(IO,' The  fit is for a polynomial of degree ',Degree:2);
  Writeln(IO,' You may wish to try a higher degree. Remember');
  Writeln(IO,' that a lower degree is better if the regression');
  Writeln(IO,' equation is suitable'^J);
  Writeln(IO,' The polynomial correlation coefficients are :'^J);
  Writeln(IO,' Intercept         ',B[1]:15);
  Write(IO,'Correlation':31,'Std Error':19);
  Writeln(IO,'T-Statistic':21);
  Writeln(IO,'Coefficient':31,'of Coeff.':19);
  Write(IO,'-----------          --------':49);
  Writeln(IO,'-----------':21);
  For J := 2 to NP do
    Begin
      K := J-1;
      Write(IO,' For X**',K:2,'        ',B[J]:15,'    ',Berr[K]:15);
      Writeln(IO,'         ',T[K]:7:3);
    End;
  Writeln(IO,^J' Calculated F-Test : ',FX:15,^J);
  Writeln(IO,' Standard Error of Y estimate :',S1R:10:6);
  Writeln(IO,' Standard Error of Y as a ');
  Writeln(IO,'  Percent of Y average  : ',PS3R:8:4);
  Writeln(IO,' Standard Error of Y as a ');
  Writeln(IO,'  Percent of the Range of Y values  : ',PS4R:8:4);
{
      CALCULATE CORRELATION COEFFICIENT
}
  RCP := RSUM/SSUM;
  RCP := RCP*100.00;
  IF(RCP > 100.0) Then RCP := 100.0;
  Writeln(IO,' The Correlation Coefficient (R-squared) is: ',RCP:8:4,'%'^J);
  For J := 1 to Nobs do
    Begin
      Index := (Nvar*MR) + J;
      Data^.XX[Index] := -2.0E35;
    End;
  Pause;
End;

Overlay Procedure Orthogonal;
Var
  Regs                   : Array[0..50] of Real;
  Temp,Temp1,Temp2,Temp3 : Real;
  Flag0,Flag1,Flag2,Flag3: Boolean;
  MaxDeg                 : Integer;
  X0,XN,A,BB,SumY,Xdiv   : Real;

Function Divisor(K: Integer): Real;
Begin
  If K <> 0 then
    Begin
      L := K + Nobs + 2;
      W := 1.0;
      While (L <> Nobs + 1) do
        Begin
          L := L - 1;
          W := W*L;
        End;
      L := Nobs - K;
      While L <> Nobs do
        Begin
          L := L + 1;
          W := W/L;
        End;
      Divisor := W/(K*2+1);
    End
  Else
      Divisor := Nobs+1;
End;

Procedure Two;
Begin
  J := 11;
  For I := 0 to MaxDeg do
    Begin
      Regs[J] := Regs[J]/Divisor(I);
      J := J + 1;
    End;
  MaxDeg := MaxDeg + 1;
End;

Procedure NewDegree;
Var
 X1,X2,X3,Index : Integer;
Begin
  GotoXY(1,7);
  Write(' Enter the Degree for this Problem: '^['J');Readln(Degree);
   X1 := MaxDeg + 12;
   X2 := MaxDeg*2 + X1 + 3;
   For J := X1 to X2 do Regs[J] := 0.0;
   J := 11;
   I := 11 + MaxDeg + 1;
   K:= I + MaxDeg + 1;
   Regs[I] := 1.0;
   I := I + 1;
   Regs[I] := 1.0;
   Regs[K] := Regs[J];
   Regs[K+1] := Regs[J+1];
   Regs[1] := 1;
   Flag0 := False;
   M := 12;
End;

Procedure Clear(Var TF: Boolean);
Begin
  TF := False;
End;

Procedure SetFlag(Var TF: Boolean);
Begin
  TF := True;
End;

Procedure CheckFlag;
Begin
  Case Flag0 of
    True: Clear(Flag0);
    False: SetFlag(Flag0);
  End;
End;

Procedure Five;
Var Done : Boolean; X: Real;
Begin
Done := False;
Repeat
  M := M + 1;
  I := MaxDeg + 11 + 1;
  If Flag0 = True then I := I + 1;
  Temp1 := (Nobs+Regs[1]+1.0)*Regs[1]/(Nobs - Regs[1]);
  X := -Temp1/(Regs[1] + 1);
  While Round(Regs[I]) <> 0 do
    Begin
      Regs[I] := Regs[I]*X;
      I := I + 2;
    End;
  I := MaxDeg + 11;
  If Flag0 = False then I := I + 1;
  X := (Regs[1]*2 +1)*Nobs/(Regs[1]+1.0)/(Nobs - Regs[1]);
  While Round(Regs[I+1]) <> 0 do
    Begin
      I := I + 1;
      Temp1 := X*Regs[I];
      I := I + 1;
      Regs[I] := Regs[I] + Temp1;
    End;
  Regs[1] := Regs[1] + 1.0;
  I := MaxDeg + 12;
  If Flag0 = True then I := I + 1;
  K := I + MaxDeg + 1;
  CheckFlag;
  Temp1 := K - MaxDeg*2.0 - Regs[1];
  While 15 > Temp1 do
    Begin
      Regs[K] := Regs[K] + Regs[M]*Regs[I];
      K := K + 2;
      I := I + 2;
      Temp1 := K - MaxDeg*2.0 - Regs[1];
    End;
Until Regs[1] = Degree;
End;

Function Factr(X: Real) : Real;
Begin
  If X > 0 then
    Factr := X*Factr(X-1)
  Else
    Factr := 1;
End;

Procedure Binomial(Var X: Real);
  Var Y,Z,T : Real;
Begin
  If (Regs[3] = 0) or (Regs[2] = 0) or (Regs[2] = Regs[3]) then
      X := 1
  Else
      Begin
        Y := Factr(Regs[3]); Z := Factr(Regs[2]);
        Temp3 := Regs[3] - Regs[2];
        T := Factr(Temp3);
        X := Y/(Z*T);
      End;
End;

Procedure Coefficients;
Var
  N1,N2,N3: Integer; X,Z: Real;
Begin
  N1 := MaxDeg*2 + 13 ;
  N2 := N1 + Degree;
  Regs[7] := N1;
  N  := MaxDeg+12;
  Regs[2] := 0;
  Regs[3] := 0;
  For J := N1 to N2 do
    Begin
      Regs[10] := 1;
      X := 0.0;
      For I := J to N2 do
        Begin
          Binomial(Z);
          X := X + Z*Regs[I]*Regs[10];
          Regs[10] := Regs[10]*BB;
          Regs[3] := Regs[3] + 1;
        End;
      Regs[2] := Regs[2] + 1;
      Regs[3] := Regs[2];
      Regs[N] := X;
      N := N + 1;
    End;
    N3 := N2 - N1;
    N1 := MaxDeg  + 12;
    N2 := N1 + N3 ;
    X := 1.0;
    For J := N1 to N2 do
      Begin
        Regs[J] := Regs[J]*X;
        X := A*X;
      End;
    B[1] := Regs[N1];
    For J := N1+1 to N2 do
        B[J+1-N1] := Regs[J];
End;

Procedure Initialize;
Var
   X,Z: Real;
Begin
Clear(Flag0);
ClrScr; For J := 0 to 50 do Regs[J] := 0.0;
Writeln(' This section of SUPERSTAT is for the fitting of Orthogonal ');
Writeln(' Polynomials.  These routines are for data  which the     ');
Writeln(' Independent variables are '^['&dB EVENLY SPACED '^['&d@');
Write(' Which column contains your X values? ');Readln(ICX);
Nobs := Nobs - 1;
Write(' Which column contains your Y values? ');Readln(IY);
X0 := GetX(1,ICX);  XN := GetX(Nobs+1,ICX);
A := -2.0/(XN - X0); BB := (XN + X0)/(XN - X0);
Write(' Enter the Maximum degree of the Polynomial: ');Readln(MaxDeg);
For I := 0 to Nobs do
  Begin
    Regs[11] := Regs[11] + GetX(I+1,IY);
    X  := -2.0*I/Nobs  + 1;
    Z := 1.0;
    Regs[12] := Regs[12] + X*GetX(I+1,IY);
    For J := 1 to MaxDeg - 1 do
      Begin
        Temp := X;
        Temp1   := -(Nobs + J + 1)*J*Z;
        Temp2   := Temp1 +(Nobs - 2*I)*(J*2+1)*X;
        X := (Temp2/(J+1))/(Nobs-J);
        Z := Temp;
        Regs[13+J-1] :=  Regs[13+J-1] + GetX(I+1,IY)*X;
      End;
   End;
End;

Begin
  Initialize;
  Two;
  Repeat
    NewDegree;
    Five;
    Coefficients;
    Analysis(Nobs+1);
    ClrScr;
    Writeln(IO,IFF,' The  fit is for a polynomial of degree ',Degree:2);
    Writeln(IO,' You may wish to try a higher degree. Remember');
    Writeln(IO,' that a lower degree is better if the regression');
    Writeln(IO,' equation is suitable'^J);
    Writeln(IO,' The polynomial correlation coefficients are :'^J);
    Writeln(IO,' Intercept         ',B[1]:15);
    For J := 2 to Degree+1 do
      Begin
        K := J-1;
        Writeln(IO,' A',K,'                ',B[J]:15);
      End;
{
      CALCULATE CORRELATION COEFFICIENT
}
    IF(RCP > 100.0) Then RCP := 100.0;
    Writeln(IO,' The Correlation Coefficient (R-squared) is: ',RCP:8:4,'%'^J);
      GotoXY(1,24);ClrEol;
      Write(' Try Another Degree');
    Until Not Yes;
    Nobs := Nobs + 1;
End;

Overlay Procedure LinearRegression;
Var
  ZZ : Real;

Procedure Exponential;
Begin
  For J := 1 to Nobs do
    Begin
      If (GetX(J,IY) > 0) then
        PutX(Ln(GetX(J,IY)),J,Nvar+1)
      Else
        PutX(0.0,J,Nvar+1);
    End;
  IY := Nvar + 1;
  Nvar := Nvar + 1;
  InitCorr;
  Nvar := Nvar - 1;
End;

Procedure PowerCurve;
Begin
  For J := 1 to Nobs do
    Begin
      If (GetX(J,IY) > 0) then
        PutX(Ln(GetX(J,IY)),J,Nvar+1)
      Else
        PutX(0.0,J,Nvar+1);
      If (GetX(J,IIC) > 0) then
        PutX(Ln(GetX(J,IIC)),J,Nvar+2)
      Else
        PutX(0.0,J,Nvar+2);
    End;
  IY := Nvar + 1;
  ICOL[1] := Nvar + 2;
  Nvar := Nvar + 2;
  InitCorr;
  Nvar := Nvar - 2;
End;

Procedure Logarithmic;
Begin
  For J := 1 to Nobs do
    Begin
      If (GetX(J,IIC) > 0.0) then
        PutX(Ln(GetX(J,IIC)),J,Nvar+1)
      Else
        PutX(0.0,J,Nvar+1);
    End;
  Icol[1] := Nvar + 1;
  Nvar := Nvar + 1;
  InitCorr;
  Nvar := Nvar - 1;
End;


Begin
  ClrScr;
  Writeln(' Multiple Linear Regression Section '^@);
  Writeln(' Any variable may be correlated against as many as 30');
  Writeln(' Other variables in your data set');
  Pause;
  ChooseColumns(IY,ICX,Icol,1);
  IY1 := IY;
  CheckSkips;
  Imod := 1;
  DF1  := NOBS-Skips-ICX-1;
  XDF  := 1.000/Int(DF1);
  TTT  := (1.93237+1.45140*XDF)/(1.0-0.73411*XDF) ;
  If(ICX = 1) Then
    Begin
      IIC := Icol[1];
      ClrScr;
      Writeln(' You have entered a simple 2 variable problem ');
      Writeln('  Y= f(X). Choose one of the following models...');
      Writeln('  1 - Linear      :  Y = a + b*X   (default)  ');
      Writeln('  2 - Exponential :  Y = a*EXP(b*X)   ');
      Writeln('  3 - Power       :  Y = a*X**b       ');
      Writeln('  4 - Logarithmic :  Y = a + b*LN(X)'^J^J);
      Writeln('  Enter your choice ... '); Read(Kbd,Ch);
      Writeln; Imod := Ord(Ch) - 48;
      Case Imod of
        2: Exponential;
        3: PowerCurve;
        4: Logarithmic;
      End;
    End;
{
      MOVE APPROPRIATE CORRELATION COEFFICIENTS INTO CORRELATION
      MATRIX FOR THE PROBLEM AS STATED
}
    If (Skips > 0) and (Imod < 2) then
      InitCorr;
    Scorr(Icol,Nobs,ICX+1,RR);
{
      CALCULATE CORRELATION COEFFICIENTS
}

  ZZ:= Int(Nobs - Skips);
  For I := 1 to ICX do
    RY[I] := RR[I,ICX+1];

  MatrixInvert(RR,ICX);
{
      MULTIPLY THE INVERSE OF THE X CORRELATION MATRIX
      BY THE Y CORRELATION MATRIX TO GET THE A MATRIX
}
  MatrixMultiply(RY,RR,ICX,A);
  SY := Diagonal[ICX+1];
  For I := 1 to ICX do
    Begin
      SX[I] := Diagonal[I];
      SS[I] := Vertical[I];
    End;
{
      CALCULATE CORRELATION COEFFICIENTS
}
  For I := 1 to ICX do
    B[I+1] := A[I]*SQRT(SY/SX[I]);
  B[1] := AVG[IY];
  For I := 1 to ICX do
    Begin
      B[1] := B[1]-B[I+1]*AVG[Icol[I]];
    End;
{
      CALCULATE ANALYSIS OF VARIANCE FIGURES
}
  ColumnSum(YY,IY,Nobs);
  SSB0 := YY*YY/ZZ;
  SSB := B[1]*YY;
  For J := 1 to ICX do
    Begin
      N := Icol[J];
      DotProduct(Z,N,IY,Nobs);
      SSB := SSB + B[J+1]*Z;
    End;
  SST := SSB - SSB0 ;
  SSR := SY - SST ;
  S2R := SSR/Int(Nobs -Skips -ICX - 1);
  S1R := 0.0;
  If(S2R >= 0) then S1R := SQRT(S2R);
  MST := SST/ICX ;

{
      F TEST
}
  If(S2R <> 0.0) then
    FX := MST/S2R
  Else
    FX  :=  1.0E+08;
{
      CALCULATE STD. Writeln OF COEFFICIENTS AND T SCORES
}
  For I := 1 to ICX do
    Begin
      BERR[I] := S1R/SQRT(SX[I]);
      If(BERR[I] = 0) THEN
        T[I]  :=  1.0E+08
      ELSE
        T[I] := B[I+1]/BERR[I];
      BMIN[I] := B[I+1]-TTT*BERR[I];
      BMAX[I] := B[I+1]+TTT*BERR[I];
    End;
{
      CALCULATE R-SQUARED FOR THE REGRESSION
}
  RC1 := SST/SY;
  RC2 := 1.000-(1.000-RC1)*((ZZ-1.0000)/(ZZ-Int(ICX))) ;
  RC1 := 100.0*RC1;
  If(RC1 > 100.0) then RC1 := 100.;
  RC2 := 100.0*RC2 ;
  If(RC2 > 100.0) then RC2 := 100. ;
{
      CALCULATE PREDICTED Y VALUES AND RESIDUALS
}
  If ICX = 1 then
    Begin
      If(IMOD = 2) or (IMOD = 3) then B[1] := Exp(B[1]);
      Icol[1] := IIC;
      IY := IY1;
    End;
  For  J := 1 to NOBS do
    Begin
      If ICX = 1 then
        Begin
          Case Imod of
           1: W := B[1] + B[2]*GetX(J,IIC);
           2: W := B[1]*Exp(B[2]*GetX(J,IIC));
           3: W := B[1]*Power(GetX(J,IIC),B[2]);
           4: W := B[1] + B[2]*Ln(GetX(J,IIC));
          End;
          PutX(W,J,FC);
        End
      Else
        Begin
          PutX(B[1],J,FC);
          For K := 1 to ICX do
            Begin
              N := Icol[K];
              W := GetX(J,FC)+B[K+1]*GetX(J,N);
              PutX(W,J,FC);
            End;
        End;
      W := GetX(J,IY)-GetX(J,FC);
      PutX(W,J,RC);
    End;
  Yrange := Maximum(IY,Nobs) - Minimum(IY,Nobs);
  PSDR1 := 100.0*S1R/Yrange;
  PSDR := 100.0*S1R/AVG[IY];
  If(IMOD = 3) or (IMOD = 4) then ICOL[1] := IIC ;
  Writeln(IO,IFF,'Linear Regression Results':43,^J^J);
  Case Imod of
    2: Writeln(IO,'  Exponential Curve: Y  =  a*EXP(b*X) '^J);
    3: Writeln(IO,'  Power Curve      : Y  =  a*X**b '^J);
    4: Writeln(IO,'  Logarithmic Curve: Y  =  a + b*ALOG(X) '^J);
  End;
  Writeln(IO,'Standard':25,'   Correl   Regression    Std Error of    T');
  Write(IO,'  X     Mean    Deviation   X vs Y   Coefficient   Reg. Coeff.');
  Writeln(IO,'  Value');
  Write(IO,'-----   -----   ----------  ------   -----------   ----------');
  Writeln(IO,'  ------');
  For J := 1 to ICX do
    Begin
      N := ICOL[J];
      Write(IO,N:2,'   ',AVG[N]:8:3,'  ',StDev[N]:8:4,'     ',RY[J]:6:4);
      Writeln(IO,'   ',B[J+1]:11,'   ',Berr[J]:11,'  ',T[J]:6:2);
      Writeln(IO,^J' Limits of Regression       Upper   ',Bmax[J]:11);
      Writeln(IO,' Coefficient (95% Conf.)    Lower   ',Bmin[J]:11,^J);
    End;
    Writeln(IO,^J'Standard':22);
    Writeln(IO,'  Y     Mean    Deviation');
    Writeln(IO,'-----   -----   ----------');
    Writeln(IO,' ',IY:2,'   ',Avg[IY]:8:3,'  ',StDev[IY]:8:4);
    If(IMOD = 2) or (IMOD = 3) then
      Bmod := EXP(B[1])
    Else
      Bmod := B[1];
    Writeln(IO,^J' Y  Intercept',':  ':26,Bmod:15);
    Writeln(IO,^J' R-squared For Correlation          :  ',RC1:8:4);
    Writeln(IO,^J' R-squared, Adjusted for Degrees');
    Writeln(IO,' of Freedom',':  ':28,RC2:8:4);
    Writeln(IO,^J' Std. error of Y estimate           :  ',S1R:8:5);
    Writeln(IO,^J' Std. error as a percent of Y Avg.  :  ',PSDR:8:4);
    Writeln(IO,^J' Std. error as a percent of the ');
    Writeln(IO,' Range of Y values                  :  ',PSDR1:8:4);
    Pause;
    Writeln(IO,^J'         Analysis of Variance Table for the Regression'^J);
    Writeln(IO,'             Degrees      Sum Of       Mean      Calculated');
    Writeln(IO,'   Source   of Freedom    Squares     Squares     F - Value');
    Writeln(IO,'    ----     --------      -----       -----       -------');
    Write(IO,' Regression    ',ICX:2,'        ',SST:7:4,'     ',MST:7:4);
    Writeln(IO,'      ',FX:7:4);
    Writeln(IO,^J' Residual     ',DF1:2,'        ',SSR:7:4,'     ',S2R:7:4);
    DF := NOBS-Skips-1;
    Writeln(IO,^J'  Total,  '^J^M'  Corrected    ',DF:2,'        ',SY:7:4,^J);
  If Skips > 0 then
    Begin
      UnsKip;
      Skips := 0;
      InitCorr;
    End;
  Pause;
End;

(* {$I B:Plot2.Ovl}  *)
Procedure Output;
Begin
  ClrScr;GotoXY(1,8);Iff := ' ';
  Writeln(' Do you wish output to go to :');
  Writeln(' 1 - Your Screen');
  Writeln(' 2 - The Printer, or');
  Writeln(' 3 - A disc File ');
  Read(KBD,Ch);
  If (Ch = '1') then
    Begin
      Assign(IO,'CON:');
      Reset(IO);
    End
  Else  If (Ch = '2') then
    Begin
      Assign(IO,'LST:');
      Reset(IO);
      Iff := ^L;
    End
  Else If (Ch = '3') then
    Begin
      Repeat
        Write(' Enter a filename for output: ');
        Read(FileName);
        Assign(IO,FileName);
{$I-}   Rewrite(IO);  IOCheck; {$I+}
      Until Not IOErr;
      FileWrite := True;
    End;

End;

Procedure ResidualTable;
{
      RESIDUAL ANALYSIS SECTION
}
Begin
  ClrScr;
  Writeln(IO,Iff,'Residual Table:40'^J^J);
  Write(IO,'  Observation ');
  Writeln(IO,'Percent':15);
  Write(IO,'     Number         Y-Actual        Y-Calc');
  Writeln(IO,'        Residual    Deviation'^J^J);
  For J := 1 to Nobs do
    Begin
      If GetX(J,IY) <> 0.0 then
        W := 100.0*(GetX(J,FC)-GetX(J,IY))/GetX(J,IY);
      Write(IO,'      ',J:2,'       ',GetX(J,IY):15,'  ');
      Write(IO,GetX(J,FC):15,' ',GetX(J,RC):15,' ');
      If GetX(J,IY) <> 0.0 then
        Writeln(IO,W:8:3)
      Else
        Writeln(IO,'********');
    End;
    Pause;
End;

Procedure ModifyData;
  Var
    A: Char;  XX : Real; MaxC, NN : Integer;
Begin
  N := Nvar; MaxC := Nvar;
  Repeat
    ClrScr;
    Writeln(^J^J' You may modify colums of data as follows :');
    Writeln(' 1 - Multiply two columns or a constant and a column');
    Writeln(' 2 - Take the Natural log of a column');
    Writeln(' 3 - Divide two columns or a constant and a column');
    Writeln(' 4 - Exponentiate a column ( Col B = Exp(Col A))  ');
    Writeln(' 5 - Raise a column to a power ');
    Writeln(' 6 - Add/Subtract two columns or a constant and a column');
    Writeln(' 7 - Normalize a column of data');
    Writeln(' 8 - Save data to a DIF file');
    Writeln(' 0 - Quit the Data Modification section ');
    Writeln(' Enter your choice : ');
  Read(Kbd,CH);
  NN := Nvar + 1;
  If(Ch <> '0') and (Ch <> '7') then
    Begin
      ClrScr;
      Writeln(' You have ',Nvar:2,' columns of data. You may  ');
      Writeln(' Overwrite data in any of these columns if you wish');
      Writeln(' or store data in the next available empty column, ');
      Writeln(' which is column ',NN:2,^J^J);
    End;
    Case Ch of
{
      MULTIPLY TWO COLUMNS TOGETHER
}
   '1': Begin
         Writeln(' You may:');
         Writeln(' 1 - multiply a column by a constant, or');
         Writeln(' 2 - multiply two columns together');
         Writeln(^J' Enter your choice: ');
         Read(Kbd,A);Writeln;
{
      MULTIPLY A COLUMN BY A CONSTANT
}
         Case A of
        '1': Begin
              Write(' Column C= Column A * Constant X, enter A,X,C :');
              Readln(M,XX,N);
              For I := 1 to Nobs do
               PutX(XX*GetX(I,M),I,N);
             End;
{
      MULTIPLY ONE COLUMN BY ANOTHER
}
      '2': Begin
             Write(' Column C = Column A*Column B, Enter A,B,C : ');
             Readln(L,M,N);
             For I := 1 to Nobs do
               PutX(GetX(I,L)*GetX(I,M),I,N);
           End;
        End;
      End;
{
      TAKE THE LOG OF A COLUMN
}
   '2': Begin
          Write('Column C = Ln(Column A) - Enter A and C');
          Readln(M,N);
          For I := 1 to Nobs do
            PutX(Ln(GetX(I,M)),I,N);
        End;
{
      DIVIDE TWO COLUMNS
}
   '3': Begin
          Writeln(' You may:');
          Writeln(' 1 - divide a column by a constant, or');
          Writeln(' 2 - divide one column by another');
          Writeln(^J' Enter your choice: ');
          Read(Kbd,A);Writeln;
         Case A of
       '1': Begin
              Write(' Column C= Column A / Constant X, enter A,X,C :');
              Readln(M,XX,N);
              For I := 1 to Nobs do
                PutX(GetX(I,M)/XX,I,N);
            End;
       '2': Begin
              Write(' Column C = Column A/Column B, Enter A,B,C : ');
              Readln(L,M,N);
              For I := 1 to Nobs do
                PutX(GetX(I,L)/GetX(I,M),I,N);
            End;
         End;
       End;
{
      EXPONENTIAL FUNCTION
}
   '4': Begin
          Write(' Column C = EXP( Column A) - Enter A, C : ');
          Readln(L,N);
          For I := 1 to Nobs do
            PutX(Exp(GetX(I,L)),I,N);
        End;
{
      RAISE A COLUMN TO A POWER
}
   '5': Begin
          Write(' Column C = Column A**X ,  Enter A, C, X: ');
          Readln(L,N,XX);
          For I := 1 to Nobs do
            PutX(Power(GetX(I,L),XX),I,N);
        End;
{
      ADD TWO COLUMNS
}
   '6': Begin
         Writeln('You may :');
         Writeln(' 1 - Add/Subtract a constant to/from a column ');
         Writeln(' 2 - Add/Subtract any two columns');
         Writeln('  Enter your choice: '); Read(Kbd,A);Writeln;
         Case A of
         '1': Begin
               Writeln(' Column C = Column A + B (Constant). Enter A, B, C ');
               Write(' Enter a Negative constant to subtract !: ');
               Readln(L,XX,N);
               For I := 1 to Nobs do
                  PutX(GetX(I,L) + XX,I,N);
              End;
         '2': Begin
               Writeln(' Column C = Column A + Column B, Enter A, B, C ');
               Write(' To subtract, enter a the negative of B : ');
               Readln(L,M,N);K := Abs(M);
               For I := 1 to Nobs do
                 Begin
                   If M < 0 then
                     PutX(GetX(I,L) + GetX(I,K),I,N)
                   Else
                     PutX(GetX(I,L) - GetX(I,K),I,N);
                 End;
              End;
        End;
        End;
{
}
    '7': Begin
           Writeln(' Enter the column number of the data you wish to');
           Writeln(' normalize, and a column number where the results should');
           Write(' be placed: ');Readln(L,N);
           For I := 1 to Nobs do
             Begin
               XX := (GetX(I,L) - Avg[L])/StDev[L];
               PutX(XX,I,N);
             End;
         End;
    '8': Begin
          ClrScr;
          DifSave;
         End;
      End;
      If N > Nvar then Nvar := Nvar + 1;
      If N > MaxC then Maxc := N;
  Until Ch = '0';
  InitCorr;

End;
Procedure InputData;

{
       This is a simple program to list out the directory of the
       current (logged) drive.
}
type
  Char15arr            = array [ 1..15 ] of Char;
  String20             = string[ 20 ];
  Suffix               = array [1..3] of Char;
Var
  NamR                 : Array[1..20] of String20;
  I,J,Number           : Integer;
  Open                 : Boolean;
  Types                : array[1..2] of Suffix;

Procedure DirList(Descriptor: Suffix; Var Found: Boolean);
var
  DTA                  : array [ 1..43 ] of Byte;
  Mask                 : Char15arr;
  Error, Default       : Integer;
  Drive                : Array[1..2] of Char;
begin { main body of program DirList }
  ClrScr;
  Drive[1] := 'B'; Drive[2] := 'A';
  FillChar(DTA,SizeOf(DTA),0);        { Initialize the DTA buffer }
  FillChar(Mask,SizeOf(Mask),0);      { Initialize the mask }
  FillChar(NamR,SizeOf(NamR),0);      { Initialize the file name }

  Result.AX := $1A00;         { Function used to set the DTA }
  Result.DS := Seg(DTA);      { store the parameter segment in DS }
  Result.DX := Ofs(DTA);      {   "    "      "     offset in DX }
  MSDos(Result);              { Set DTA location }
  Result.AX := $1900;
  MSDos(Result);
  Default := (Result.AX and $FF ) + 1;
  Error := 0;
  WriteLn(Descriptor,' Files on Drive ',Drive[Default]);
  WriteLn;
  Mask := ' :\????????.   ';
  Mask[1] := Drive[Default];
  Mask[13] := Descriptor[1];
  Mask[14] := Descriptor[2];
  Mask[15] := Descriptor[3];
  Result.AX := $4E00;
  Result.DS := Seg(Mask);
  Result.DX := Ofs(Mask);
  Result.CX := 22;
  MSDos(Result);
  Error := Result.AX and $FF;
  I := 1;
  J := 1;
  if (Error = 0) then
    repeat
      NamR[J][I] := Chr(Mem[Seg(DTA):Ofs(DTA)+29+I]);
      I := I + 1;
    until not (NamR[J][I-1] in [' '..'~']) or (I>20);

  NamR[J][0] := Chr(I-1);
  while (Error = 0) do begin
    Error := 0;
    Result.AX := $4F00;
    Result.CX := 22;
    MSDos( Result );
    Error := Result.AX and $FF;
    J := J + 1;
    I := 1;
    repeat
      NamR[J][I] := Chr(Mem[Seg(DTA):Ofs(DTA)+29+I]);
      I := I + 1;
    until not (NamR[J][I-1] in [' '..'~'] ) or (I > 20);
    NamR[J][0] := Chr(I-1);
  end ;
  GotoXY(1,5);
  FileName[1] := Drive[Default];
  FileName[2] := ':';
  FileName[3] := '\';
  Found := False;
  If J-1 > 0 then
    Begin
      For I := 1 to J - 1  do
        Writeln(Chr(64+I),':   ',Namr[I]);
      GotoXY(1,22);
      Write(' Choose the file you want by letter, hit <ESC> to exit ');
      Read(Kbd,Ch);
      Ch := UpCase(Ch);
      Number := Ord(Ch) - 64;
      If Number in[1..J-1] then
        Begin
          Insert(NamR[Number],FileName,4);
          Found := True;
        End;
    End
  Else
    Begin
      Writeln(' No ',Descriptor,' Files Found on Drive ',Drive[Default]);
      Found := False ;
      Pause;
    End;
End;

Begin
  Repeat
    Types[1] := 'DIF'; Types[2] := 'PRN';
    ClrScr;  Skips := 0;
    Writeln(' YOU HAVE THE FOLLOWING OPTIONS : '^J);
    Writeln(' 1 - READ NEW DATA FROM A DIF FILE ');
    Writeln(' 2 - READ NEW DATA FROM A LOTUS PRN FILE ');
    Writeln(' 3 - INPUT NEW DATA VIA THE SPREADSHEET, OR ');
    Writeln(' 4 - MODIFY/ADD TO EXISTING DATA VIA THE SPREADSHEET');
    Writeln(^J' Enter your Choice '); Read(Kbd,CH);
  Until Ch in['1','2','3','4'];
  L := Ord(Ch) - 48;
  Case L of
    1,2: Begin
           FileName := '                    ';
           DirList(Types[L],Open);
           If Open then
             Begin
               GotoXY(1,15);
               Write(' Opening File : ');
               Writeln(FileName);
   {$I-}       Reset(InFile);  IOCheck;  {$I+}
             End;
           If Open then
             If (L = 1) then DFRead else PRN_Read;
         End;
    3  : Begin
           Release(HeapTop);
           Mark(HeapTop); New(Data);
           MC := 45; MR := 200; SC := 43;
           FC := 44; RC := 45;
           InitializeArray;
           SpreadSheet;
           Open := True;
         End;
    4  : Begin
           SpreadSheet;
           Open := True;
         End;
  End;
  If Open then InitCorr;
  Xnobs := Int(Nobs);
End;

Begin
{
      THE DATA ARRAYS ARE DIMENSIONED TO READ UP TO 200 OBSERVATIONS
      OF 50 DIFFERENT VARIABLES, AND ANALYZE THE EFFECTS OF UP TO 20
      OF THESE VARIABLES ON ANY INDEPENDANT VARIABLE (Y)
}
  ClrScr;  PlotCall := 0;
  GotoXY(1,10); TextBackground(0); TextColor(15);
  Writeln(' Welcome to SUPERSTAT, the Statistical analysis program. This ');
  Writeln(' program currently reads data in DIF format files such as Lotus');
  Writeln(' 1,2,3 or Visicalc can create, or Lotus PRN files .');
  Writeln(' You may also enter data interactively via a spreadsheet-like');
  Writeln(' Interface. This data can then be saved to a DIF File. Output');
  Writeln(' is initially directed to the screen. ');
  Assign(IO,'CON:');
  Reset(IO);
  Mark(HeapTop);
  New(Data);
  Pause;
  Repeat
    ClrScr;GotoXY(1,8);
    Writeln(' Choose from the Following Options');
    Writeln(' 1 - Enter/Modify a Data Set');
    Writeln(' 2 - Correlation Matrix');
    Writeln(' 3 - Multiple Linear Regression');
    Writeln(' 4 - Polynomial Regression');
    Writeln(' 5 - Orthogonal Polynomials');
    Writeln(' 6 - Simplex Non-Linear curve fitting');
    Writeln(' 7 - Residual Table');
    Writeln(' 8 - Modify Data');
    Writeln(' 9 - Change Output Device');
{   Writeln(' P - Plot Data');       }
    Writeln(' X - Exit the Program');
    Read(Kbd,Ch);
    Case Ch of
      'X','x'         : Begin
                         Ch := ' ';
                         GotoXY(1,22);ClrEol;
                         Write(' Sure you want to Exit');
                         If Yes then
                           Ch := 'X';
                        End;
      '1'             : InputData;
      '2'             : CorrelationMatrix;
      '3'             : LinearRegression;
      '4'             : PolyRegr;
      '5'             : Orthogonal;
      '6'             : Simplex;
      '7'             : ResidualTable;
      '8'             : ModifyData;
      '9'             : OutPut;
{     'P','p'         : DataPlot;   }
    End;
  Until Ch = 'X';
  GotoXY(1,22);ClrEol;Write(' Do you need to save your data');
  If Yes then DifSave;
  If FileWrite then
    Begin
      Flush(IO);
      Close(IO);
    End;
  Release(HeapTop);
End.