Trailing-Edge
-
PDP-10 Archives
-
decus_20tap5_198111
-
decus/20-0137/stepr/stepr.for
There is 1 other file named stepr.for in the archive. Click here to see a list.
C WESTERN MICHIGAN UNIVERSITY
C STEPR.F4 (FILE NAME ON LIBRARY DECTAPE)
C STEPR, 1.3.2 (CALLING NAME, SUBLST #)
C STEPWISE REGRESSION
C ORIGINAL PROGRAM CAME FROM WAYNE STATE UNIVERSITY. LATER
C IT WAS MODIFIED BY SAM ANEMA AND R.R. BARR.
C LIBRARY DECTAPE PROG. USED: USAGE.MAC
C FORWMU PROGS. USED: TTYPTY, EXISTS, DEVCHG, PRINTS, XPRODH
C APLIB PROGS. USED: IO, GETFOR, FISHER
C INTERNAL SUBR. USED: LANG2, TRANS, LANG1, REGGR, TNUM
C ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
C
C
DIMENSION X(70),SX(70),NN(70),SXY(70,70),XMEAN(70),IDEN(16)
DIMENSION FMD(70),IFT(96),IT(70),KTRAN(50)
DIMENSION KVAR1(50),KVAR2(50),KVAR3(50),CONST(50),IFOPT(15)
DIMENSION VAR(70),IRROR(2)
INP=2
IOUT=3
IND=-4
NOUTD = -1
WRITE(NOUTD,2021)
2021 FORMAT('1',20X,'WMU STEPWISE REGRESSION'//)
C CALL USAGE('STEPR ')
C---------------TTYPTY RETURNS ZERO - TTY JOB, MINUS ONE - BATCH JOB;
C--------------- ICODE USED AS ARG. IN CALL IO.
CALL TTYPTY(ICODE)
C---------------IDV1, IDV2 ARE RETURNED. OTHER ARGS. ARE INPUT.
C--------------- 1 MEANS OUTPUT? PRINTS, 0 MEANS INPUT? PRINTS.
CALL IO(IOUT,IDV1,NOUTD,IND,1,ICODE)
1111 NTR=0
C---------------ITYPE REFERS TO MISSING DATA REPLACEMENT METHOD
C--------------- SEE ST. 110.
ITYPE=3
CALL IO(INP,IDV2,NOUTD,IND,0,ICODE)
C---------------IFT, ISTD ARE RETURNED. 96=NO. OF OBJ. TIME FORMAT
C---------------WORDS. (6 LINES) 2 MEANS F-TYPE FORMAT ONLY.
CALL GETFOR(NOUTD,IND,IFT,ISTD,96,2)
IF(ISTD.EQ.1)IFT(1)='(10F)'
1 WRITE(NOUTD,2)
2 FORMAT(' ENTER NUMBER OF VARIABLES'/)
READ(IND,3)NVAR
3 FORMAT(I)
IF(NVAR.LT.71.AND.NVAR.GT.1)GO TO 5
CALL TTYPTY(ICC)
IF(ICC.EQ.-1)CALL EXIT
GO TO 1
5 WRITE(NOUTD,4)
4 FORMAT(' ENTER IDENTIFICATION IF DESIRED OTHERWISE RETURN'/)
READ(IND,11)(IDEN(I),I=1,16)
11 FORMAT(16A5)
REWIND 1
121 WRITE(NOUTD,6)
6 FORMAT(' ENTER OPTIONS OR TYPE "HELP"'/)
READ(IND,7)(IT(I),I=1,70)
7 FORMAT(70A1)
IF(IT(1).NE.'H'.OR.IT(2).NE.'E'.OR.IT(3).NE.'L')GO TO 91
WRITE(NOUTD,9991)
9991 FORMAT(/,' OPTION TABLE',//' CODE DISCRIPTION',/,
1 ' ---- -----------',/,' MIS* MISSING DATA',/,
2 ' TRA* TRANSFORMATION OF VARIABLES',/,
3 ' MVS MEANS STANDARD DEV., AND VARIABLES',/,
4 ' XPR RAW SUMS OF SQUARES AND CROSS PRODUCTS',/,
5 ' COV VARIANCE - COVARIANCE MATRIX',/,
6 ' COR CORRELATION, FISHER Z, AND T MATRICES',/,
7 ' FVL F - VALUES',/,
8 ' FOR FORCE VARIABLES',/,
9 ' ELM ELIMINATE VARIABLES',/,
A ' RES RESIDUALS',/,
B ' ANA ANALYSIS OF VARIANCE AT FINAL STAGE',/,
C ' DUR DURBIN - WATSON STATISTICS',/,
D ' ZER FORCE REGRESSION THROUGH ZERO',//,
E ' * READ PERTINATE REMARKS IN SECTION 3.0',/,
F ' OF LIBRARY PROGRAM WRITEUP #1.3.2 BEFORE USING.',//)
GO TO 121
91 DO 8 I=1,15
8 IFOPT(I)=0
CALL LANG1(IT,IFOPT)
C---------------IFOPT(I)=1 MEANS USER CHOOSES THIS OPTION.
C--------------- IFOPT(1), ... , IFOPT(15) MEANS MVS, COR, COV,
C--------------- RES, ANA, DUR, ZER, TRA, MIS, OBJ, FVL, FOR, ELM,
C--------------- XPR, TOL RESPECTIVELY - OBJ AND TOL DO NOT APPEAR
C--------------- IN WRITE UP IN HELP LIST.
IF(IFOPT(9).NE.1)GO TO 120
WRITE(NOUTD,102)
102 FORMAT(' IS THERE MORE THAN ONE MISSING DATA SYMBOL ?'/)
READ(IND,103)MDS
103 FORMAT(A3)
IF(MDS.EQ.'YES') GO TO 104
WRITE(NOUTD,106)
106 FORMAT(' ENTER MISSING DATA SYMBOL'/)
READ(IND,107)FMD(1)
107 FORMAT(10F)
DO 108 I=2,NVAR
108 FMD(I)=FMD(1)
GO TO 111
104 WRITE(NOUTD,109)
109 FORMAT(' ENTER MISSING DATA SYMBOLS'/)
READ(IND,107)(FMD(I),I=1,NVAR)
111 WRITE(NOUTD,110)
110 FORMAT(' HOW WOULD YOU LIKE TO COMPENSATE FOR MISSING DATA? TYPE:'
1/' 1 - TO REPLACE MISSING DATA BY MEAN VALUE'/
2' 2 - TO DELETE THE OBSERVATION'/)
READ(IND,3)ITYPE
IF(ITYPE .EQ.1.OR.ITYPE.EQ.2)GO TO 120
CALL TTYPTY(ICC)
IF(ICC.EQ.-1)CALL EXIT
GO TO 111
120 MAX=0
C---------------IFOPT(8)=1 MEANS CHOOSE TRANS
IF(IFOPT(8).NE.1)GO TO 30
C---------------INDICES OF VARS. I, J, K AS PRESCIBED BY USER
C--------------- (SEE WRITE UP) FOR THE TRANSFORMATIONS ARE RETURNED
C--------------- IN KVAR1, KVAR2, KVAR3 RESPECTIVELY FOR LATER
C--------------- USE BY SUBR. TRANS.
C---------------KTRAN(I)=ONE OF 15 POSSIBLE TRANSF. LISTED IN SECTION
C--------------- 3.0 OF WRITE UP UNDER OPTION 2; CONST=A AS
C--------------- DESCRIBED IN THE SAME SECTION.
CALL LANG2(KTRAN,CONST,KVAR1,KVAR2,KVAR3,IT,NTR)
C---------------NTR=NO. OF TRANSF. AND IS RETURNED BY SUBR. LANG2.
DO 9 I=1,NTR
IF(KVAR1(I).GT.MAX)MAX=KVAR1(I)
9 CONTINUE
C---------------NVAR=NO. OF VARS, MAX=MAXIMUM INDEX OF VARS.
C--------------- GENERATED BY TRANSF.
30 NVTR=MAX0(NVAR,MAX)
40 REWIND 20
REWIND 1
DO 299 J=1,NVTR
SX(J)=0.0
NN(J)=0
DO 300 I=1,NVTR
SXY(I,J)=0.0
300 CONTINUE
299 CONTINUE
K=0
C---------------IDV2 IS RETURNED BY IO.
IF(IDV2.EQ.'TTY')GO TO 13
WRITE(NOUTD,14)
14 FORMAT(' YOUR DATA IS BEING READ'/)
GO TO 17
13 IF(ISTD.EQ.1)WRITE(NOUTD,15)
15 FORMAT(' ENTER DATA (AT MOST 10 PER LINE)'/)
IF(ISTD.EQ.0)WRITE(NOUTD,16)
16 FORMAT(' ENTER DATA'/)
17 IRROR(1)=0
IRROR(2)=0
GO TO (301,351,391),ITYPE
C---------------HERE IF WE REPLACE MISSING DATA BY MEAN VALUE,
C--------------- SEE ST. 110+1.
301 READ(INP,IFT,END=340)(X(I),I=1,NVAR)
303 K=K+1
DO 304 I=1,NVAR
C---------------FIND(I)=MISSING DATA SYMBOLS SEE ST. 109+1.
IF(X(I).EQ.FMD(I)) GO TO 304
SX(I)=SX(I)+X(I)
NN(I)=NN(I)+1
304 CONTINUE
WRITE(20)(X(I),I=1,NVAR)
GO TO 301
340 REWIND 20
DO 341 I=1,NVAR
XMEAN(I)=SX(I)/FLOAT(NN(I))
341 SX(I)=0.0
C---------------K=NO. OF ABS., SEE ST. 303.
DO 343 J=1,K
READ(20)(X(I),I=1,NVAR)
DO 342 I=1,NVAR
IF(X(I).EQ.FMD(I))X(I)=XMEAN(I)
342 CONTINUE
CALL TRANS(KTRAN,CONST,KVAR1,KVAR2,KVAR3,X,NTR,IRROR,KAR)
C---------------KAR RETURNED BY SUBR. TRANS AND KAR INDICATES
C--------------- ERRORS DETECTED BY SUBR. TRANS.
C---------------X(I) CONTAINS NEW VALUES GENERATED BY SUBR. TRANS.
IF(KAR.EQ.1)GO TO 343
WRITE(1)(X(I),I=1,NVTR)
CALL XPRODH(X,SX,SXY,NVTR,70)
343 CONTINUE
K=K-IRROR(1)-IRROR(2)
GO TO 2020
C---------------HERE WE DELETE ENTIRE OBS. IF ONE VAR. HAS
C--------------- MISSING DATA.
351 READ(INP,IFT,END=2020)(X(I),I=1,NVAR)
360 DO 354 I=1,NVAR
IF(X(I).EQ.FMD(I))GO TO 351
354 CONTINUE
CALL TRANS(KTRAN,CONST,KVAR1,KVAR2,KVAR3,X,NTR,IRROR,KAR)
IF(KAR.EQ.1)GO TO 351
K=K+1
WRITE(1)(X(I),I=1,NVTR)
CALL XPRODH(X,SX,SXY,NVTR,70)
GO TO 351
C---------------HERE IF NO MISSING DATA
391 READ(INP,IFT,END=2020)(X(I),I=1,NVAR)
395 CALL TRANS(KTRAN,CONST,KVAR1,KVAR2,KVAR3,X,NTR,IRROR,KAR)
IF(KAR.EQ.1)GO TO 391
K=K+1
WRITE(1)(X(I),I=1,NVTR)
CALL XPRODH(X,SX,SXY,NVTR,70)
GO TO 391
2020 IF(K.GT.1)GO TO 20210
IF(IRROR(1).GT.0)WRITE(NOUTD,7654)IRROR(1)
IF(IRROR(2).GT.0)WRITE(NOUTD,7655)IRROR(2)
WRITE(NOUTD,7656)K
7656 FORMAT(' ONLY ',I2,' OBSERVATIONS REMAIN',/)
IF(ICODE.EQ.-1)CALL EXIT
GO TO 1111
C---------------VAR(I)=VARIANCE OF ITH VAR., SX(I)=MEAN OF ITH VAR.
C---------------SEE ST. 2023+3.
20210 DO 2022 I=1,NVTR
SX(I)=SX(I)/FLOAT(K)
C---------------THIS IS COMP. FORMULA FOR VARIANCE TAKING NOTE
C--------------- OF THE FACT THAT ON THE LINE ABOVE
C--------------- SX(I)=SX(I)/FLOAT(K). SEE STEPR FOLDER.
VAR(I)=(SXY(I,I)-FLOAT(K)*SX(I)**2)/FLOAT(K-1)
IF(ABS(VAR(I)).LT..0000001)GO TO 19
2022 CONTINUE
NVAR=NVTR
WRITE(IOUT,12)(IDEN(I),I=1,16)
12 FORMAT('1',16A5/)
WRITE(IOUT,2026)K
2026 FORMAT(/' THE NUMBER OF OBSERVATIONS IS',I5//)
IF(IRROR(1).GT.0)WRITE(IOUT,7654)IRROR(1)
7654 FORMAT(' NUMBER OF OBS. REJECTED BECAUSE OF ILLEGAL TRANS. LOG IS'
1,I5/)
IF(IRROR(2).GT.0)WRITE(IOUT,7655)IRROR(2)
7655 FORMAT(' NUMBER OF OBS. REJECTED BECAUSE OF ILLEGAL TRANS.',
1' DIVISION IS',I5/)
IF(IFOPT(1).NE.1)GO TO 50
WRITE(IOUT,2023)
2023 FORMAT(4X,'VARIABLE',4X,'MEAN',8X,'VARIANCE',6X,'STD. DEV.'/)
DO 2024 I=1,NVAR
SD=SQRT(VAR(I))
WRITE(IOUT,2025)I,SX(I),VAR(I),SD
2025 FORMAT(6X,I2,2X,3G)
2024 CONTINUE
50 IF(IFOPT(14).NE.1)GO TO 55
WRITE(IOUT,2030)
2030 FORMAT('1RAW SUMS OF SQUARES AND CROSS PRODUCTS'/)
KSIZE=5
C---------------RETURNED BY TTYPTY. ICODE=-1 MEANS BATCH JOB
IF(ICODE.EQ.-1)KSIZE=10
KF=1
2032 KL=MIN0(KF+KSIZE-1,NVAR)
WRITE(IOUT,2043)(I,I=KF,KL)
DO 2034 I=KF,NVAR
KLL=MIN0(KL,I)
C---------------RETURNED BY XPRODH
WRITE(IOUT,2035)I,(SXY(I,J),J=KF,KLL)
2035 FORMAT(/I3,1X,12G12.5)
2034 CONTINUE
IF(NVAR.EQ.KL)GO TO 55
KF=KF+KSIZE
WRITE(IOUT,2066)
GO TO 2032
C---------------IFOPT(7)=1 MEANS FORCE EQ. THRU ZERO
55 IF(IFOPT(7).EQ.1)GO TO 3220
C---------------CALCULATE AND WRITE OUT COVARIANCE MATRIX.
DO 2040 I=1,NVAR
DO 2040 J=1,I
C---------------SEE COMMENTS WITH ST. 20210+2 ABOVE
2040 SXY(I,J)=(SXY(I,J)-FLOAT(K)*SX(I)*SX(J))/FLOAT(K-1)
C---------------IFOPT(3)=1 MEANS PRINT COV. MATRIX
IF(IFOPT(3).NE.1)GO TO 60
WRITE(IOUT,2041)
2041 FORMAT('1COVARIANCES'/)
KSIZE=5
IF(ICODE.EQ.-1)KSIZE=10
KF=1
2042 KL=MIN0(KF+KSIZE-1,NVAR)
WRITE(IOUT,2043)(I,I=KF,KL)
2043 FORMAT(4X,12(5X,I2,5X))
DO 2044 I=KF,NVAR
KLL=MIN0(KL,I)
WRITE(IOUT,2045)I,(SXY(I,J),J=KF,KLL)
2045 FORMAT(/I3,1X,12F12.5)
2044 CONTINUE
IF(NVAR.EQ.KL)GO TO 60
KF=KF+KSIZE
WRITE(IOUT,2066)
GO TO 2042
C---------------HERE FROM ST. 55 I.E. USER WANTS FORCE THRU ZERO;
C--------------- MATRIX SUBMITTED TO SUBROUTINE REGGR IS SIMILAR TO
C--------------- CORRELATION MATRIX BUT MEANS ARE NOT SUBTRACTED
C--------------- I.E., WE USE RAW SUMS OF SQUARES AND CROSS PRODUCTS
C--------------- SEE STEPR FOLDER.
3220 REWIND 20
DO 3221 I=1,NVAR
3221 X(I)=SQRT(SXY(I,I))
DO 3222 I=1,NVAR
DO 3222 J=1,I
SXY(I,J)=SXY(I,J)/(X(I)*X(J))
3222 SXY(J,I)=SXY(I,J)
WRITE(20)(SX(I),I=1,NVAR)
WRITE(20)(X(I),I=1,NVAR)
DO 3223 I=1,NVAR
3223 WRITE(20)(SXY(I,J),J=I,NVAR)
C---------------CALL REGGR
GO TO 70
60 REWIND 20
DO 2050 I=1,NVAR
X(I)=SQRT(VAR(I)*FLOAT(K-1))
C---------------FORM CORRELATION MATRIX
DO 2050 J=1,I
SXY(I,J)=SXY(I,J)/SQRT(VAR(I)*VAR(J))
IF(I.EQ.J)SXY(I,J)=1.0
2050 SXY(J,I)=SXY(I,J)
WRITE(20)(SX(I),I=1,NVAR)
WRITE(20)(X(I),I=1,NVAR)
DO 2049 I=1,NVAR
2049 WRITE(20)(SXY(I,J),J=I,NVAR)
C---------------IFOPT(2)=1 MEANS PRINT CORRELATION MATRIX
IF(IFOPT(2).NE.1)GO TO 70
WRITE(IOUT,2051)
2051 FORMAT('1CORRELATIONS'/)
KSIZE=7
IF(ICODE.EQ.-1)KSIZE=14
KF=1
2052 KL=MIN0(KF+KSIZE-1,NVAR)
WRITE(IOUT,2053)(I,I=KF,KL)
2053 FORMAT(4X,14(4X,I2,3X))
DO 2054 I=KF,NVAR
KLL=MIN0(KL,I)
WRITE(IOUT,2055)I,(SXY(I,J),J=KF,KLL)
2055 FORMAT(/I3,1X,14F9.5)
2054 CONTINUE
IF(NVAR.EQ.KL)GO TO 2058
KF=KF+KSIZE
WRITE(IOUT,2066)
GO TO 2052
2058 WRITE(IOUT,2059)
2059 FORMAT('1T-VALUES'/)
KF=1
N=NVAR
2060 KL=MIN0(KF+KSIZE-1,N)
WRITE(IOUT,2061)(I,I=KF,KL)
2061 FORMAT(4X,14(3X,I2,4X))
DO 2062 I=KF,N
KLL=MIN0(KL,I)
DO 2063 J=KF,KLL
IF(I.EQ.J)GO TO 2064
SXY2=1.0-SXY(I,J)**2
IF(SXY2.LE.0.0)GO TO 2064
X(J)=SXY(I,J)*SQRT(FLOAT(K-2)/SXY2)
GO TO 2063
2064 X(J)=0.0
2063 CONTINUE
WRITE(IOUT,2065)I,(X(J),J=KF,KLL)
2065 FORMAT(/I3,1X,14F9.4)
2062 CONTINUE
IF(NVAR.EQ.KL)GO TO 2068
KF=KF+KSIZE
WRITE(IOUT,2066)
2066 FORMAT(///)
GO TO 2060
2068 WRITE(IOUT,2069)
2069 FORMAT('1Z-VALUES'/)
KF=1
2070 KL=MIN0(KF+KSIZE-1,N)
WRITE(IOUT,2061)(I,I=KF,KL)
DO 2072 I=KF,N
KLL=MIN0(KL,I)
DO 2073 J=KF,KLL
IF(I.EQ.J)GO TO 2074
IF(1.0-SXY(I,J))2071,2074,2071
2071 SXY2=(1.0+SXY(I,J))/(1.0-SXY(I,J))
IF(SXY2.LE.0.0)GO TO 2074
X(J)=SQRT(FLOAT(K-3))*.5*ALOG(SXY2)
GO TO 2073
2074 X(J)=0.0
2073 CONTINUE
WRITE(IOUT,2065)I,(X(J),J=KF,KLL)
2072 CONTINUE
IF(NVAR.EQ.KL)GO TO 70
KF=KF+KSIZE
WRITE(IOUT,2066)
GO TO 2070
C---------------K=NO. OF OBS. (SEE ST. 303)
70 CALL REGGR(SXY,K,VAR,SX,NVAR,IFOPT)
WRITE(NOUTD,71)
71 FORMAT(' DO YOU WISH TO REANALYZE THE SAME DATA?'/)
READ(IND,18)KWHER
18 FORMAT(A3)
IF(KWHER.EQ.'YES')GO TO 70
C---------------ST. 1111 IS ON FIRST PAGE. NEXT PROMPTING
C--------------- IS INPUT?
GO TO 1111
19 WRITE(NOUTD,21)I
21 FORMAT(' VARIABLE ',I3,' HAS ZERO VARIANCE'/)
IF(ICODE.EQ.-1)CALL EXIT
GO TO 1111
END
C---------------KTRAN, CONST, KVAR1, KVAR2, KVAR3, NTR ARE
C--------------- INPUT. X, KAR ARE OUTPUT.
SUBROUTINE TRANS(KTRAN,CONST,KVAR1,KVAR2,KVAR3,X,NTR,IRROR,KAR)
DIMENSION X(70),IRROR(2)
DIMENSION KTRAN(50),KVAR1(50),KVAR2(50),KVAR3(50),CONST(50)
KAR=0
IF(NTR.EQ.0)RETURN
DO 10 I=1,NTR
GO TO (20,21,22,23,24,25,26,27,28,29,30,31,32,33,34),KTRAN(I)
C---------------SEE OPTION 2 PAGE 5 OF WRITE UP
20 X(KVAR1(I))=-X(KVAR2(I))
GO TO 10
21 X(KVAR1(I))=X(KVAR2(I))
GO TO 10
22 X(KVAR1(I))=X(KVAR2(I))**CONST(I)
GO TO 10
23 X(KVAR1(I))=X(KVAR2(I))+X(KVAR3(I))
GO TO 10
24 X(KVAR1(I))=X(KVAR2(I))*X(KVAR3(I))
GO TO 10
25 IF(X(KVAR2(I)).GT.0.0)GO TO 40
IRROR(1)=IRROR(1)+1
KAR=1
RETURN
40 X(KVAR1(I))=ALOG(X(KVAR2(I)))
GO TO 10
26 IF(X(KVAR2(I)).GT.0.0)GO TO 41
IRROR(1)=IRROR(1)+1
KAR=1
RETURN
41 X(KVAR1(I))=ALOG10(X(KVAR2(I)))
GO TO 10
27 X(KVAR1(I))=EXP(X(KVAR2(I)))
GO TO 10
28 X(KVAR1(I))=10.0**X(KVAR2(I))
GO TO 10
29 X(KVAR1(I))=X(KVAR2(I))+CONST(I)
GO TO 10
30 X(KVAR1(I))=X(KVAR2(I))*CONST(I)
GO TO 10
31 X(KVAR1(I))=SIN(X(KVAR2(I)))
GO TO 10
32 X(KVAR1(I))=COS(X(KVAR2(I)))
GO TO 10
33 IF(X(KVAR3(I)).NE.0.0)GO TO 42
IRROR(2)=IRROR(2)+1
KAR=1
RETURN
42 X(KVAR1(I))=X(KVAR2(I))/X(KVAR3(I))
GO TO 10
34 X(KVAR1(I))=X(KVAR2(I))-X(KVAR3(I))
10 CONTINUE
RETURN
END
C---------------KTRAN, CONST, KVAR1, KVAR2, KVAR3, NTR ARE OUTPUT;
C--------------- IT IS INPUT. NTR=NO. OF TRANSF.,
SUBROUTINE LANG2(KTRAN,CONST,KVAR1,KVAR2,KVAR3,IT,NTR)
DIMENSION KTRAN(50),CONST(50),KVAR1(50),KVAR2(50)
DIMENSION KVAR3(50),IT(70)
NOUTD=-1
IND=-4
NTR=0
WRITE(NOUTD,2001)
2001 FORMAT(' ENTER TRANSFORMATIONS'/)
2000 READ(IND,100)(IT(I),I=1,70)
100 FORMAT(70A1)
IER=0
K=1
DO 1 I=1,70
IF(IT(I).EQ.' ')GO TO 1
IT(K)=IT(I)
K=K+1
1 CONTINUE
NK=K-1
DO 2 I=K,70
2 IT(I)=' '
IF(IT(1).EQ.'E'.AND.IT(2).EQ.'N'.AND.IT(3).EQ.'D')GO TO 3000
NTR=NTR+1
IF(IT(1).EQ.'(')GO TO 500
3 WRITE(NOUTD,1100)
1100 FORMAT(' ERROR IN TRANSFORMATION'/)
CALL TTYPTY(ICC)
IF(ICC.EQ.-1)CALL EXIT
NTR=NTR-1
GO TO 2000
500 I=2
CALL TNUM(I,X,IT,IER)
IF(IER.NE.0) GO TO 3
KVAR1(NTR)=X+.01
I=I+1
IF(IT(I).EQ.'=')GO TO 501
GO TO 3
501 I=I+1
IF(IT(I).EQ.'(')GO TO 300
IF(IT(I).EQ.'E')GO TO 400
IF(IT(I).EQ.'L')GO TO 600
IF(IT(I).EQ.'1')GO TO 700
IF(IT(I).EQ.'S')GO TO 800
IF(IT(I).EQ.'C')GO TO 900
IF(IT(I).EQ.'-')GO TO 1000
GO TO 3
300 I=I+1
CALL TNUM(I,X,IT,IER)
IF(IER.NE.0)GO TO 3
KVAR2(NTR)=X+.01
I=I+1
IF(IT(I).EQ.'*')GO TO 310
IF(IT(I).EQ.'+')GO TO 320
IF(IT(I).EQ.'-')GO TO 330
IF(IT(I).EQ.'/')GO TO 340
IF(IT(I).EQ.' ')GO TO 350
GO TO 3
310 I=I+1
IF(IT(I).EQ.'*')GO TO 304
IF(IT(I).EQ.'(')GO TO 306
IF(IT(I).EQ.'+'.OR.IT(I).EQ.'-'.OR.IT(I).EQ.'.')GO TO 308
IF(IT(I).LE.'9'.AND.IT(I).GE.'0')GO TO 308
GO TO 3
304 I=I+1
CALL TNUM(I,X,IT,IER)
IF(IER.NE.0)GO TO 3
CONST(NTR)=X
KTRAN(NTR)=3
GO TO 2000
306 I=I+1
CALL TNUM(I,X,IT,IER)
IF(IER.NE.0)GO TO 3
KVAR3(NTR)=X+.01
KTRAN(NTR)=5
GO TO 2000
308 CALL TNUM(I,X,IT,IER)
IF(IER.NE.0)GO TO 3
CONST(NTR)=X
KTRAN(NTR)=11
GO TO 2000
320 I=I+1
IF(IT(I).EQ.'(')GO TO 322
IF(IT(I).EQ.'.')GO TO 326
IF(IT(I).GE.'0'.AND.IT(I).LE.'9')GO TO 326
GO TO 3
322 I=I+1
CALL TNUM(I,X,IT,IER)
IF(IER.NE.0)GO TO 3
KVAR3(NTR)=X+.01
KTRAN(NTR)=4
GO TO 2000
326 CALL TNUM(I,X,IT,IER)
IF(IER.NE.0)GO TO 3
CONST(NTR)=X
KTRAN(NTR)=10
GO TO 2000
330 I=I+1
IF(IT(I).EQ.'(')GO TO 332
IF(IT(I).EQ.'.')GO TO 336
IF(IT(I).GE.'0'.AND.IT(I).LE.'9')GO TO 336
GO TO 3
332 I=I+1
CALL TNUM(I,X,IT,IER)
IF(IER.NE.0)GO TO 3
KVAR3(NTR)=X+.01
KTRAN(NTR)=15
GO TO 2000
336 I=I-1
CALL TNUM(I,X,IT,IER)
IF(IER.NE.0)GO TO 3
CONST(NTR)=X
KTRAN(NTR)=10
GO TO 2000
340 I=I+1
IF(IT(I).EQ.'(')GO TO 341
GO TO 3
341 I=I+1
CALL TNUM(I,X,IT,IER)
IF(IER.NE.0)GO TO 3
KVAR3(NTR)=X+.01
KTRAN(NTR)=14
GO TO 2000
350 KTRAN(NTR)=2
GO TO 2000
400 I=I+4
IF(IT(I-3).EQ.'*'.AND.IT(I-2).EQ.'*'.AND.IT(I-1).EQ.'(')GO TO 450
GO TO 3
450 CALL TNUM(I,X,IT,IER)
IF(IER.NE.0)GO TO 3
KVAR2(NTR)=X+.01
KTRAN(NTR)=8
GO TO 2000
600 I=I+1
IF(IT(I).EQ.'N')GO TO 620
IF(IT(I).EQ.'O')GO TO 640
GO TO 3
620 I=I+1
IF(IT(I).EQ. '(')GO TO 625
GO TO 3
625 I=I+1
CALL TNUM(I,X,IT,IER)
KVAR2(NTR)=X+.01
KTRAN(NTR)=6
GO TO 2000
640 I=I+2
IF(IT(I-1).EQ.'G'.AND.IT(I).EQ.'(')GO TO 645
GO TO 3
645 I=I+1
CALL TNUM(I,X,IT,IER)
IF(IER.NE.0)GO TO 3
KVAR2(NTR)=X+.01
KTRAN(NTR)=7
GO TO 2000
700 I=I+4
IF(IT(I-3).EQ.'0'.AND.IT(I-2).EQ.'*'.AND.IT(I-1).EQ.'*'.AND.IT(I).
1 EQ.'(')GO TO 750
GO TO 3
750 I=I+1
CALL TNUM(I,X,IT,IER)
IF(IER.NE.0)GO TO 3
KVAR2(NTR)=X+.01
KTRAN(NTR)=9
GO TO 2000
800 I=I+3
IF(IT(I-2).EQ.'I'.AND.IT(I-1).EQ.'N'.AND.IT(I).EQ.'(')GO TO 850
GO TO 3
850 I=I+1
CALL TNUM(I,X,IT,IER)
IF(IER.NE.0)GO TO 3
KVAR2(NTR)=X+.01
KTRAN(NTR)=12
GO TO 2000
900 I=I+3
IF(IT(I-2).EQ.'O'.AND.IT(I-1).EQ.'S'.AND.IT(I).EQ.'(')GO TO 950
GO TO 3
950 I=I+1
CALL TNUM(I,X,IT,IER)
IF(IER.NE.0)GO TO 3
KVAR2(NTR)=X+.01
KTRAN(NTR)=13
GO TO 2000
1000 I=I+1
IF(IT(I).EQ.'(')GO TO 1050
GO TO 3
1050 I=I+1
CALL TNUM(I,X,IT,IER)
IF(IER.NE.0)GO TO 3
KVAR2(NTR)=X+.01
KTRAN(NTR)=1
GO TO 2000
3000 RETURN
END
C---------------IT INPUT, IFOPT OUTPUT.
SUBROUTINE LANG1(IT,IFOPT)
DIMENSION ICOM(15,3),IT(70),IFOPT(15)
DATA (ICOM(I,1),I=1,10)/'M','C','C','R','A','D','Z','T','M','O'/
DATA (ICOM(I,2),I=1,10)/'V','O','O','E','N','U','E','R','I','B'/
DATA (ICOM(I,3),I=1,10)/'S','R','V','S','A','R','R','A','S','J'/
DATA (ICOM(I,1),I=11,15)/'F','F','E','X','T'/
DATA (ICOM(I,2),I=11,15)/'V','O','L','P','O'/
DATA (ICOM(I,3),I=11,15)/'L','R','M','R','L'/
K=1
C**** WMU AM: 1.3.2-2, MTO, 1-FEB-78
DO 1I=1,69
C**** END = LANG1, STMT #1-4
IF(IT(I).EQ.' ')GO TO 1
IT(K)=IT(I)
K=K+1
1 CONTINUE
NK=K
DO 2 I=K,70
2 IT(I)=' '
C**** WMU AM: 1.3.2-2, MTO, 1-FEB-78
DO 3 I=4,NK
IF(IT(I).NE.','.AND.IT(I).NE.' ')GO TO 3
DO 3 J=1,15
IF(IT(I-3).EQ.ICOM(J,1).AND.IT(I-2).EQ.ICOM(J,2).AND.IT(I-1).EQ.
1ICOM(J,3))IFOPT(J)=1
C**** END = LANG1, STMT #2+5
3 CONTINUE
RETURN
END
C---------------I, X, IER ARE OUTPUT; IT IS INPUT
SUBROUTINE TNUM(I,X,IT,IER)
DIMENSION IT(70),IN(15),ITP(3)
DO 10 J=I,70
IF(IT(J).EQ.'.')GO TO 10
IF(IT(J).GE.'0'.AND.IT(J).LE.'9')GO TO 10
IF(IT(J).EQ.' '.OR.IT(J).EQ.')')GO TO 20
IF(J.NE.I)IER=1
IF(IT(J).EQ.'+')GO TO 10
IF(IT(J).EQ.'-')GO TO 10
IER=1
10 CONTINUE
20 NPT=0
KF=J-1
DO 30 J=I,KF
IF(IT(J).EQ.'.')NPT=NPT+1
30 CONTINUE
IF(NPT.GT.1)IER=1
IF(IER.NE.0)RETURN
NCH=15-(KF-I)
DO 40 J=1,15
40 IN(J)=' '
K=I
DO 50 J=NCH,15
IN(J)=IT(K)
50 K=K+1
ENCODE(15,60,ITP)IN
60 FORMAT(15A1)
DECODE(15,70,ITP)X
70 FORMAT(F15.0)
I=KF+1
RETURN
END
C PDP-10 VERSION OF ESSO STEPWISE REGRESSION ANALYSIS - PART 2
C BLOCK E - READ MATRIX, ELIMINATE VARIABLES
C---------------ALL ARGS. ARE INPUT. ARRAY=EITHER CORRELATION
C--------------- MATRIX OR A MATRIX DERIVED FROM CORRELATION COEFFS.
C--------------- WHERE MEANS ARE NOT SUBTRACTED. SEE ST. 3220 AND
C--------------- ST. 2049+1 IN MAIN PROG.; NDATA=NO. OF OBS.;
C--------------- SIGMA=VARIANCES; AVG=MEANS, NTERS=NO. OF VAR.,
C--------------- IFOPT=USER SELECTED OPTIONS.
SUBROUTINE REGGR(ARRAY,NDATA,SIGMA,AVG,NTERS,IFOPT)
C**L.M.#1.3.2-1 RRB 26SEP77
C**ADDED ARRAY "IFORCE" TO REPLACE "INDEX" BEFORE DO LOOP 104
C**SO LIST OF FORCED VARIABLES IS NOT OVER-WRITTEN.
DIMENSION COEN(70),SIGCO(70),INDEX(70),SIGMA(70),STUDT(70),
1DATA(70),AVG(70),BETA(70),ARRAY(70,70),INDX2(70),NAME(70),
2LAME(70),KAME(70),KFR(70),IFOPT(15),IFORCE(70)
IND=-4
NOUTD=-1
IOUT=3
TOL=.0001
REWIND 20
REWIND 1
INP=5
COE=0.0
EFIN=0.0
EFOUT=0.0
NTERM=NTERS
NOELM=0
NOFRC=0
OBSER=FLOAT(NDATA)
WRITE(NOUTD,2000)
2000 FORMAT(///' WHICH IS THE DEPENDENT VARIABLE?'/)
READ(IND,2001)NDVAR
2001 FORMAT(I)
IF(IFOPT(15).NE.1)GO TO 2019
WRITE(NOUTD,2022)
2022 FORMAT(' ENTER TOLERANCE'/)
READ(IND,2003)TOL
2019 IF(IFOPT(11).NE.1)GO TO 3862
WRITE(NOUTD,2002)
2002 FORMAT(' ','ENTER F-VALUE FOR ENTERING A VARIABLE.'/)
READ(IND,2003)EFIN
2003 FORMAT(F)
WRITE(NOUTD,2004)
2004 FORMAT(' ','ENTER F-VALUE FOR OMITTING A VARIABLE.'/)
READ(IND,2003)EFOUT
3862 IF(IFOPT(13).NE.1)GO TO 2006
WRITE(NOUTD,2005)
2005 FORMAT(' ','HOW MANY VARIABLES WOULD YOU LIKE TO ELIMINATE?'/)
READ(IND,2001)NOELM
IF(NOELM)2006,2006,2007
2007 WRITE(NOUTD,2008)
2008 FORMAT(' ','WHICH ARE THEY?(MAX: 20 PER LINE)'/)
READ(IND,2009)(INDX2(I),I=1,NOELM)
2009 FORMAT(20I)
2006 NOVAR=NDVAR
IF(IFOPT(12).NE.1)GO TO 715
WRITE(NOUTD,2010)
2010 FORMAT(' ','ENTER NUMBER OF VARIABLES TO BE FORCED INTO THE ',
1'REGRESSION.'/)
READ(IND,2001)NOFRC
IF(NOFRC)2011,715,2011
2011 WRITE(NOUTD,2008)
C**L.M.#1.3.2-1 RRB 26SEP77 - CHANGED "INDEX" TO "IFORCE"
READ(IND,2009)(IFORCE(I),I=1,NOFRC)
715 LDVAR = NDVAR
ISTRE = 0
KDV1=0
IDW=1
SDEV2=0.0
IDIF=0
KTR = 0
SDIF2 = 0.0
IF(NOELM)28,28,302
28 NOELM=0
GO TO 70
302 IF(NOELM)994,994,304
994 NPCRD = 2
NCOL=6
995 WRITE(IOUT,10)NPCRD,NCOL
10 FORMAT(24H ERROR ON PARAMETERCARD,I3,7X,4HCOL.,I2)
GO TO 1004
304 IF(NOELM - 1)994,303,72
72 L15=NOELM - 1
DO 451 K8 = 1,L15
K9=K8+1
C---------------SORT VARS. TO BE ELEMINATED. INDX2(I) ARE
C--------------- INDICES OF VARS. TO BE ELIM.
DO 451 K10=K9,NOELM
IF(INDX2(K8) - INDX2(K10))451,451,450
450 LN=INDX2(K8)
INDX2(K8)=INDX2(K10)
INDX2(K10)=LN
451 CONTINUE
303 L=NOELM+1
INDX2(L)=0
C---------------OBSER=NO. OF OBS. SEE ST. 2000-2.
70 IF(OBSER)71,73,73
71 IDW=0
IFDEV=0
OBSER= -OBSER
C---------------NTERM=NO. OF VARS. SEE ST. 2000-5. COEN(I)=
C--------------- SQRT. OF ROW SUMS OF SQUARES. (SEE ST. 3222+2
C--------------- AND 3221 IN MAIN PROG.) IF USER CHOOSES FORCE
C--------------- THRU ZERO OPTION; OTHERWISE COEN(I)=SQRT. OF
C--------------- QUANTITY VARIANCE TIMES NO. OF OBS. MINUS 1 (SEE
C--------------- ST. 60+2 IN MAIN PROG.) IF USER DOES NOT CHOOSE FORCE
C--------------- THRU ZERO OPTION. (ALSO SEE ST. 2050+2 IN MAIN PROG.)
73 READ(20)(AVG(I),I=1,NTERM)
READ(20)(COEN(I),I=1,NTERM)
LTERM=NTERM
DO 1050 K=1,NTERM
1050 NAME(K)=K
518 WRITE(IOUT,502)NDVAR
502 FORMAT(13H/VARIABLE NO.,I4,14H IS DEPENDENT///)
IF(NOELM) 700,703,700
700 WRITE(IOUT,701)(INDX2(I),I=1,NOELM)
701 FORMAT(1H0,18X,'ELIMINATE ', 30I3/(30X,30I3))
GO TO 1703
703 WRITE(IOUT,1701)NTERM
1701 FORMAT(1H ,I4,29H VARIABLES NONE ELIMINATED)
1703 J=1
L = 1
N=INDX2(1)
C---------------COMPRESS VARS. NOT TO BE ELIMINATED. KFR(I) WILL
C---------------GIVE INDICES FOR NON-ELIMINATED VARS. IN SUCH A WAY THAT
C--------------- THE KTH POSITION GIVES ROW AND COL. NO. OF MATRIX ARRAY
C--------------- FOR THE KTH VARIABLE.
DO 87 I=1,NTERM
IF(NOELM)994,85,86
86 IF(I-N)88,89,88
88 AVG(J)=AVG(I)
85 SIGMA(J)=COEN(I)
KAME(J) = NAME(I)
K6 = NAME(I)
KFR(K6) = J
IF(I - NDVAR)851,850,851
850 KDV = J
851 J=J+1
GO TO 87
89 L=L+1
N=INDX2(L)
87 CONTINUE
852 J1 = J - 1
869 IF(KDV - J1)870,871,3380
C---------------PROG. SHOULD NEVER GET TO 3380 SINCE BY ST. 850
C--------------- KDV IS ONE OF THE NON-ELIMINATED VARS.
3380 WRITE(IOUT,3381)
3381 FORMAT(66H NO. OF DEPENDENT VARIABLE IS GREATER THAN THE NUMBER OF
1 VARIABLES)
C---------------GO TO RETURN
GO TO 1004
870 R1=AVG(KDV)
R2=SIGMA(KDV)
N1=KAME(KDV)
C---------------PUT VALUES FOR DEP. VAR. IN LAST POSITION OF
C---------------NON-ELIMINATED VARS.
DO 872 I=KDV,J1-1
AVG(I)=AVG(I+1)
SIGMA(I)=SIGMA(I+1)
KAME(I)=KAME(I+1)
872 CONTINUE
AVG(J1)=R1
SIGMA(J1)=R2
KAME(J1)=N1
DO 875 I=1,NTERM
IF(KFR(I).EQ.KDV)GO TO 873
875 CONTINUE
873 DO 874 J=I+1,NTERM
874 KFR(J)=KFR(J)-1
871 M = 1
LSTRT = 1
DO 331 I=1,NTERM
READ(20)(COEN(J),J=I,NTERM)
IF(NOELM)994,333,332
332 L = LSTRT
N=INDX2(L)
IF(I - N)339,330,339
330 LSTRT = LSTRT+1
GO TO 331
339 K=M
L=LSTRT
N=INDX2(L)
DO 336 J=I,NTERM
IF(J-N) 335,337,335
335 ARRAY(M,K)=COEN(J)
3352 K=K+1
GO TO 336
337 L=L+1
N=INDX2(L)
336 CONTINUE
M=M+1
GO TO 331
333 DO 338 K=I,NTERM
338 ARRAY(I,K)=COEN(K)
331 CONTINUE
KDV1=KDV
DEFR=OBSER-1.0
C---------------NOVAR=NO. OF TERMS NOT TO BE ELIMINATED
NOVAR = NTERM - NOELM
NOVMI=NOVAR-1
DO 300 I=2,NOVAR
IM1=I-1
DO 300 J=1,IM1
300 ARRAY(I,J)=ARRAY(J,I)
IF(NOELM) 994,3312,3310
C---------------KDV1=INDEX OF DEP. VAR. SEE ST. 331+1.
3310 IF(KDV1 - NOVAR)3311,3049,3380
3311 NDVAR = KDV1
NTERM = NOVAR
GO TO 3382
3312 IF(NDVAR - NTERM)3382,3049,3380
3382 DO 3000 I=1,NTERM
T2 = ARRAY(I,NDVAR)
DO 3001 J=NDVAR,NTERM-1
3001 ARRAY(I,J)=ARRAY(I,J+1)
3000 ARRAY(I,NTERM)=T2
DO 3002 I=1,NTERM
T2=ARRAY(NDVAR,I)
DO 3003 J=NDVAR,NTERM-1
3003 ARRAY(J,I)=ARRAY(J+1,I)
3002 ARRAY(NTERM,I)=T2
C---------------NOFRC=NO. OF VARS. TO BE FORCED.
3049 IF(NOFRC)305,350,305
C**L.M.#1.3.2-1 RRB 26SEP77 - CHANGED "INDEX" TO "IFORCE"
305 WRITE(IOUT,702)(IFORCE(I),I=1,NOFRC)
702 FORMAT(23X, 'FORCE ', 30I3/(30X,30I3))
C
C BLOCK F - SOLVE SYSTEM OF EQUATIONS FOR COEFFICIENTS
313 DO 306 L=1,NOFRC
C**L.M.#1.3.2-1 RRB 26SEP77 - CHANGED "INDEX" TO "IFORCE"
K=IFORCE(L)
K= KFR(K)
DIAG=ARRAY(K,K)
DO 307 I=1,NOVAR
IF(I-K)308,307,308
308 FACTR=ARRAY(I,K)/DIAG
DO 309 J=1,NOVAR
IF(J-K)310,309,310
310 ARRAY(I,J)=ARRAY(I,J)-(FACTR*ARRAY(K,J))
309 CONTINUE
307 CONTINUE
DO 311 I=1,NOVAR
311 ARRAY(I,K)=-ARRAY(I,K)/DIAG
DO 312 J=1,NOVAR
312 ARRAY(K,J)=ARRAY(K,J)/DIAG
306 ARRAY(K,K)=1.0/DIAG
FACTR = NOFRC
DEFR = DEFR-FACTR
FLEVL=(ARRAY(K,NOVAR)**2/DIAG)*DEFR/ARRAY(NOVAR,NOVAR)
NOENT = K
NOSTP=NOFRC-1
6004 GO TO 1006
350 NOSTP=-1
6005 GO TO 1006
C
C BLOCK G-SELECT VARIABLES BY ANAL OF VARIANCE AND SOLVE EQUATIONS
1006 NOSTP=NOSTP+1
IF(ARRAY(NOVAR,NOVAR))147,14,101
147 WRITE(IOUT,1014)
1014 FORMAT(34H LAST DIAGONAL ELEMENT IS NEGATIVE)
GO TO 1004
14 WRITE(IOUT,15)
15 FORMAT(30H LAST DIAGONAL ELEMENT IS ZERO)
GO TO 1004
101 IF(ARRAY(NOVAR,NOVAR).LT.0.0)WRITE(NOUTD,2222)
2222 FORMAT(' 1R')
SIGY=SIGMA(NOVAR)*SQRT( ARRAY(NOVAR,NOVAR)/DEFR )
DEFR=DEFR-1.0
IF(DEFR) 16,103,103
16 WRITE(IOUT,17)DEFR
17 FORMAT(22H DEGREES OF FREEDOM IS,F8.2)
GO TO 1004
103 VMIN=-10000.0
VMAX=0.0
NOIN=0
DO 104 I=1,NOVMI
DIAG=ARRAY(I,I)
C**L.M.#1.3.2-1 RRB 26SEP77
C**EXCLUDE FORCED VARIABLES FROM ENTER/REMOVE TESTING
IF(DIAG.EQ.0)GO TO 104
IF(NOFRC.EQ.0)GO TO 105
DO 1991 JJ=1,NOFRC
1991 IF(KAME(I).EQ.IFORCE(JJ))GO TO 1992
GO TO 105
1992 SIG=ABS(SIGMA(I))
IF(SIG.EQ.0)GO TO 18
GO TO 134
C**END OF L.M.#1.3.2-1
105 IF(DIAG)107,104,108
107 IF(DIAG+TOL)109,104,104
109 WRITE(IOUT,12)KAME(I)
12 FORMAT( 43H NEGATIVE DIAGONAL ELEMENT BELOW TOLERANCE,,8HVARIABLE,
1I4)
GO TO 1004
108 IF(DIAG-TOL)104,110,110
110 VAR=ARRAY(I,NOVAR)*ARRAY(NOVAR,I)/DIAG
SIG=SIGMA(I)
IF(VAR)106,104,112
106 IF(SIG)129,18,122
18 WRITE(IOUT,19)SIG,KAME(I)
19 FORMAT(22H STANDARD DEVIATION IS,F14.7,13HFOR VARIABLE ,I3)
GO TO 1004
129 SIG=-SIG
GO TO 134
122 IF(VAR-VMIN)134,134,113
113 VMIN=VAR
C---------------INDEX OF VAR. CORRESPONDING TO MIN. VAR.
NOMIN=I
C---------------NOIN IS NO. OF EQS. IN REG. EQ. (SEE ST. 133+2)
134 NOIN=NOIN+1
INDEX(NOIN)=I
COEN(NOIN)=ARRAY(I,NOVAR)*SIGMA(NOVAR)/SIG
IF(DIAG.LT.0.0)WRITE(NOUTD,2223)
2223 FORMAT(' 2R')
SIGCO(NOIN)=(SIGY/SIG)*SQRT(DIAG)
1341 LAME(NOIN)=KAME(I)
IF(KAME(I)-LDVAR)104,1342,104
1342 LAME(NOIN)=J1
GO TO 104
112 IF(SIG) 18,18,146
146 IF(VAR-VMAX)104,104,114
114 VMAX=VAR
C---------------NOMAX=INDEX OF VAR. CORRESPONDING TO MAX. VAR.
NOMAX=I
104 CONTINUE
IF(NOIN) 99,115,116
115 KTR=KTR+1
WRITE(IOUT,400)SIGY
400 FORMAT(25H0 STANDARD ERROR OF Y = G)
IF(KTR - 1)123,123,455
455 WRITE(IOUT,456)
456 FORMAT(38H ERROR, NO VARIABLE CAN ENTER EQUATION)
RETURN
116 IF(IFOPT(7).NE.1)GO TO 117
118 CNST=0.0
SALPHA=0.0
TALPHA=0.0
C**** WMU AM: 1.3.2-2, MTO, 30-JAN-78
IDFALP=NDATA-NOIN
DO 6031 I=1,NOIN
6031 SIGCO(I)=SIGCO(I)*SQRT(FLOAT(IDFALP-1)/IDFALP)
GO TO 120
117 CNST=AVG(NOVAR)
SALPHA=1.0/FLOAT(NDATA)
DO 119 J=1,NOIN
L=INDEX(J)
DO 8022 I=1,NOIN
LL=INDEX(I)
SALPHA=SALPHA+AVG(L)*AVG(LL)*ARRAY(L,LL)/(SIGMA(L)*SIGMA(LL))
8022 CONTINUE
119 CNST=CNST-(COEN(J)*AVG(L))
SALPHA=SQRT(SALPHA)*SIGY
TALPHA=CNST/SALPHA
IDFALP=NDATA-NOIN-1
120 FLEV=FLEVL
IF (IFOPT(7).EQ.1) FLEV=IDFALP*FLEV/(IDFALP-1)
NOEN=NOENT
K1=K
IF(VMIN+10000.0)121,123,121
121 FLEVL=-VMIN*(DEFR+1.0)/ARRAY(NOVAR,NOVAR)
IF(EFOUT-FLEVL)123,123,124
124 K=NOMIN
DEFR=DEFR+2.0
NOENT=0
GO TO 126
123 FLEVL=VMAX*DEFR/(ARRAY(NOVAR,NOVAR)-VMAX)
IF(EFIN - FLEVL)125,1007,1007
125 K=NOMAX
NOENT=K
126 CONTINUE
COEFDT=1.0-ARRAY(NOVAR,NOVAR)
COD=COEFDT-COE
COE=COEFDT
IF(COEFDT.LT.0.0)WRITE(NOUTD,2224)
2224 FORMAT(' 3R')
RRR=SQRT(COEFDT)
127 IF(NOIN)99,128,130
130 IF(NOEN)131,131,132
132 WRITE(IOUT,401)NOSTP,KAME(K1)
401 FORMAT(///9H STEP NO. I3/20H VARIABLE ENTERING I6)
ISTRE = K1
GO TO 133
131 IF(K1 - ISTRE)8037,8038,8037
8038 IF(K1)8033,8037,8033
8037 WRITE(IOUT,402)NOSTP,KAME(K1)
402 FORMAT(///9H STEP NO. I3/19H VARIABLE REMOVED I6)
133 FPRB=FISHER(1,IDFALP,FLEV)
C****** END = REGGR, STMT #133
C**L.M.#1.3.2-1 RRB 26SEP77 - ADD TPROB TO STEPS
TPRB=FISHER(1,IDFALP,TALPHA*TALPHA)
WRITE(IOUT,403)FLEV,FPRB,SIGY,COEFDT,RRR,COD,IDFALP,CNST,SALPHA,
1TALPHA,TPRB,(LAME(J),COEN(J),SIGCO(J),J=1,NOIN)
403 FORMAT(3X,'F LEVEL',2X,F12.4,5X,'F-PROB =',2X,F8.5,
1/3X,'STANDARD ERROR OF ESTIMATE = ',
1F12.4/' COEFFICIENT OF DETERMINATION = ',F8.5/
1' COEFFICIENT OF MULTIPLE REGRESSION = ',F8.5/
1 ' INCREASE IN COEFFICIENT OF DETERMINATION = ',F8.5/
13X,'DEGREES OF FREEDOM =',I6,/,
15X,'CONSTANT',5X,'STD. ERR.',13X,'T',9X,'T-PROB'/1X,3G,F8.5/13X,
1'VARIABLE COEFFICIENT STD ERRO',
2'R OF COEF' /(16H X=I3,F15.5,F15.5))
C**END OF L.M.#1.3.2-1
128 DIAG=ARRAY(K,K)
DO 140 I=1,NOVAR
IF(I-K)141,140,141
141 FACTR= ARRAY(I,K)/DIAG
DO 143 J=1,NOVAR
IF(J-K)142,143,142
142 ARRAY(I,J)=ARRAY(I,J)-(FACTR*ARRAY(K,J))
143 CONTINUE
140 CONTINUE
DO 144 I=1,NOVAR
144 ARRAY(I,K)=-ARRAY(I,K)/DIAG
DO 145 J=1,NOVAR
145 ARRAY(K,J)=ARRAY(K,J)/DIAG
ARRAY(K,K)=1.0/DIAG
IF(NOSTP-250)1006,1450,1450
1450 WRITE(IOUT,1451)NOSTP
1451 FORMAT(25H ERROR, NUMBER OF STEPS =,I4,15H CHECK F-VALUES)
GO TO 1004
99 WRITE(IOUT,13)NOIN
13 FORMAT(15H ERROR, NOIN IS,I7)
GO TO 1004
C
C BLOCK H-WRITE FINAL OUTPUT
C COMMON ARRAY,STUDT,BETA,INDEX,NOVAR,DEFR,NOSTP,
C 1 NOIN,NOEN,NUMBR,YPRED,DEV,I,J,K,
C 2 NDATA,NTERM,IFINV,IFDEV,NOELM,IFEOF,INDX2,
C 3 COEN,SIGCO,K1,CNST,SIGY,FLIV
1007 IF( NOIN.EQ.0)GO TO 455
DO 201 I=1,NOIN
J=INDEX(I)
STUDT(I)=COEN(I)/SIGCO(I)
201 BETA(I)=ARRAY(J,NOVAR)
IF(NOEN)997,202,203
997 WRITE(IOUT,22)NOEN
22 FORMAT(14H ERROR, NOEN =,I6)
GO TO 1004
203 WRITE(IOUT,601)NOSTP,KAME(K1)
601 FORMAT(///9H STEP NO. I3/20H VARIABLE ENTERING I6)
ISTRE = K1
GO TO 204
202 IF(K1-ISTRE)8036,8035,8036
8035 IF(K1)8033,8036,8033
8033 WRITE(IOUT,8034)
8034 FORMAT ( 42H VARIABLE REMOVED SAME AS VARIABLE ENTERED )
GO TO 1004
8036 WRITE(IOUT,602)NOSTP,KAME(K1)
602 FORMAT(///9H STEP NO. I3/19H VARIABLE REMOVEDI6)
204 COEFDT=1.0-ARRAY(NOVAR,NOVAR)
COD=COEFDT-COE
IF(COEFDT.LT.0.0)TYPE 2225
2225 FORMAT(' 4R')
RRR=SQRT(COEFDT)
C****** WMU AM: 1.3.2-2, MTO, 30-JAN-78
FPRB=FISHER(1,IDFALP,FLEV)
TPRB=FISHER(1,IDFALP,TALPHA*TALPHA)
WRITE(IOUT,600)FLEV,FPRB,SIGY,COEFDT,RRR,COD,IDFALP,CNST,SALPHA,
1TALPHA,TPRB
600 FORMAT(12H F LEVEL =,F12.4,5X,'F-PROB =',2X,F8.5,
1/' STANDARD ERROR OF ESTIMATE = ',
1F12.4/' COEFFICIENT OF DETERMINATION = ',F8.5/
1' COEFFICIENT OF MULTIPLE REGRESSION = ',F8.5/
1' INCREASE IN COEFFICIENT OF DETERMINATION = ',F8.5/
13X,'DEGREES OF FREEDOM =',I6,/,
15X,'CONSTANT',5X,'STD. ERR.',13X,'T',9X,'T-PROB'/1X,3G,F8.5//
1' VARIABLE',4X,'COEFF. STD ERR OF COEFF.',
#17X,'T-VALUE',3X,'T-PROB',/37X,'STANDARDIZED COEFF.')
DO 6002 J=1,NOIN
TPRB=FISHER(1,IDFALP,STUDT(J)*STUDT(J))
C****** END = REGGR, STMT #6002-1
6002 WRITE(IOUT,6001)LAME(J),COEN(J),SIGCO(J),BETA(J),STUDT(J),TPRB
6001 FORMAT(2X,'X(',I2,')',F13.5,F13.6,F20.6,F9.3,F9.5)
IF(DEFR.EQ.0)WRITE(IOUT,8002)
8002 FORMAT(/,' NO MORE VARIABLES ARE INCLUDED BECAUSE',/,
1 ' DEGREES OF FREEDOM ARE ZERO',/)
8003 IF(IFOPT(4).NE.1.AND.IFOPT(5).NE.1.AND.IFOPT(6).NE.1)GO TO 1000
IF(IFOPT(4).NE.1)GO TO 2027
205 WRITE(IOUT,604)
604 FORMAT(///20X,'PREDICTED VS ACTUAL RESULTS' /
111X,'OBS. NO. ACTUAL PREDICTED DEVIATION' //)
2027 SSS=0.
SST=0.
DO 2070 J=1,NDATA
NUMBR=J
READ(1)(DATA(L),L=1,LTERM)
8 YPRED=CNST
SSS=DATA(LDVAR)+SSS
SST=SST+DATA(LDVAR)**2
11 DO 208I=1,NOIN
K=LAME(I)
208 YPRED=YPRED+COEN(I)*DATA(K)
DEV=DATA(LDVAR)-YPRED
51 SDEV2 = SDEV2 + DEV*DEV
33 IF(IDIF)35,36,35
35 DIF=DEV-DEVP
SDIF2 = SDIF2 + DIF*DIF
GO TO 37
36 IDIF = 1
37 DEVP = DEV
IF(IFOPT(4).NE.1)GO TO 2070
207 WRITE(IOUT,605)NUMBR,DATA(LDVAR),YPRED,DEV
2070 CONTINUE
605 FORMAT(8H I8,3F14.5)
IF(IFOPT(5).NE.1)GO TO 2472
C**** WMU AM: 1.3.2-2, MTO, 30-JAN-78
IF (IFOPT(7).NE.1) SST=SST-SSS*SSS/FLOAT(NDATA)
SSR=SST-SDEV2
NERR=IDFALP
NTOTL=IDFALP+NOIN
ERMS=SDEV2/FLOAT(NERR)
RMS=SSR/FLOAT(NOIN)
F=RMS/ERMS
FPRB=FISHER(NOIN,IDFALP,F)
C**** END = REGGR, STMT #2432-2
WRITE(IOUT,2432)
2432 FORMAT(//20X,'ANALYSIS OF VARIANCE')
WRITE(IOUT,2020)SSR,NOIN,RMS,F,FPRB,SDEV2,NERR,ERMS,SST,NTOTL
2020 FORMAT(/4X,'SOURCE',6X,'SUM OF SQ.',6X,'DF',5X,'MEAN SQ.',
16X,'F',6X,'F-PROB'/' REGRESSION',F17.6,I6,F15.5,F9.3,F9.5/,
16X,'ERROR',F17.6,
2I6,F15.5/6X,'TOTAL',F17.6,I6)
2472 IF(IFOPT(6).NE.1)GO TO 1000
DWS=SDIF2/SDEV2
DW1=4.-DWS
WRITE(IOUT,40)SDIF2,DWS,DW1
40 FORMAT(//14X,38HDURBIN-WATSON TEST FOR AUTOCORRELATION//
110X,28HSUM OF SQUARED DIFFERENCES =,F14.4/10X,5HDWS =,
2F12.8/10X,9H4 - DWS =,F12.8)
1000 RETURN
1004 CALL TTYPTY(IERCOD)
IF(IERCOD.NE.0)CALL EXIT
RETURN
END