Trailing-Edge
-
PDP-10 Archives
-
decus_20tap5_198111
-
decus/20-0137/lin/lin.for
There are 2 other files named lin.for in the archive. Click here to see a list.
C WESTERN MICHIGAN UNIVERSITY
C LIN.FOR (FILE NAME ON LIBRARY DECTAPE)
C LIN, 1.3.1 (CALLING NAME, SUBLST NO.)
C MULTIPLE LINEAR REGRESSION, PARTIAL CORRELATIONS
C PROGRAMMED BY JACK R. MEAGHER FOR THE 1620, MODIFIED BY
C B. HOUCHARD, NORM GRANT, R.R. BARR
C LIBRARY DECTAPE PROGS. USED: USAGE.MAC
C APLB10 PROGS. USED: IOB, GETFOR, FISHER
C FORWMU PROGS. USED: ZEROH, TTYPTY, ALLCOR, DEVCHG, EXISTS,
C PRINTS, XPRODH, RENAMS, MINVSQ
C INTERNAL SUBR. USED: MAINL
C ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
C
C---------------SPACE(1) IS USED IN CALL ALLCOR, SEE ST. 750+1;
C--------------- ID(16) USED FOR USER'S IDENT OF OUTPUT
C--------------- NOTF PASSED TO SUBR. GETFOR, SEE ST. 200-2,
C--------------- NOTF(48) LIMITS OBJ. TIME FORMAT TO 240 CH.
DIMENSION SPACE(1),NOTF(48),ID(16)
C---------------SUBR. IOB IN APLB10 SHARES COMMON /IOBLK/, COMMON
C--------------- /IOBLKA/, SUBR. MAINL SHARES COMMON
C--------------- /IOBLK/, COMMON /BLOCK1/
COMMON /IOBLK/IDLG,ICC,INP,IOUT,IO2,IO3,ICODE,IBNK,NAMI(2)
COMMON /IOBLKA/NAMO(2),IPJ,IPG,NCOPYS,ITYCH
COMMON/BLOCK1/ISTD,NOTF,ID
C
C***********************************************************************
C DEVICES USED:
C
C IDLG--DEVICE USED TO COMMUNICATE WITH USERS
C IT IS ALWAYS SET TO -4
C ICC---DEVICE USED TO ACCEPT USER'S RESPONSES
C IT IS ALWAYS SET TO -1
C INP---DEVICE USED TO READ DATA
C ITS LOGICAL NUMBER IS DETERMINED BY SUBROUTINE IO
C IOUT--DEVICE USED TO WRITE REPORT
C ITS LOGICAL NUMBER IS DETERMINED BY SUBROUTINE IO
C***********************************************************************
C
IDLG=-1
ICC=-4
INP=2
IOUT=3
C
WRITE(IDLG,180)
180 FORMAT(/' WMU - LINEAR MULTIPLE REGRESSION AND PARTIAL',
1' CORRELATION'//)
C CALL USAGE('LIN ')
C
C***********************************************************************
C DETERMINE IF JOB IS FROM TELETYPE OR PSEUDO-TELETYPE
C
C IF ICODE = 0 JOB IS FROM TELETYPE
C =-1 JOB IS FROM PSEUDO-TELETYPE
C***********************************************************************
C
C---------------ICODE RETURNED
CALL TTYPTY(ICODE)
C
C***********************************************************************
C GATHER INPUT/OUTPUT INFORMATION
C OUTPUT OPTION IS AVAILABLE ONLY ONCE IN THE PROGRAM
C***********************************************************************
C
C---------------1 MEANS OUTPUT? PRINTS. 0 MEANS INPUT? PRINTS
C--------------- IDLG, ICC, INP, IOUT, IO2, IO3, ICODE ARE
C--------------- INPUT THRU COMMON /IOBLK/.
CALL IOB(1)
1 CALL IOB(0)
C
C***********************************************************************
C DETERMINE TYPE OF FORMAT
C ITYPE = 2 MEANS ONLY F-TYPE FORMAT ALLOWED
C***********************************************************************
C
ITYPE=2
C---------------NOTF, ISTD ARE RETURNED. OTHER ARGS. ARE INPUT. 48=NO.
C--------------- OF OBJ. TIME FORMAT WORDS (3 LINES).
CALL GETFOR(IDLG,ICC,NOTF,ISTD,48,ITYPE)
WRITE(IDLG,200)
200 FORMAT(' ENTER HEADER'/)
READ(ICC,201) ID
201 FORMAT(16A5)
C
C READ NUMBER OF VARIABLES.
C
210 WRITE(IDLG,82)
82 FORMAT(' ENTER NUMBER OF VARIABLES, (INCLUDING DEPENDENT): ',
1$)
READ(ICC,80,END=1000)N
80 FORMAT(I)
IF(N.GE.2)GO TO 750
749 WRITE(IDLG,58)N
58 FORMAT('+NO REGRESSION POSSIBLE ON ',I3,' VARIABLES, TRY
1 AGAIN'/)
IF (ICODE.LT.0) CALL EXIT
GO TO 210
C
C GET CORE ALLOCATED.
C
750 MAX=2*N*(N+2)
CALL ALLCOR(MAX,IERR,IBASE,SPACE(1))
IF(IERR.NE.0)GO TO 749
I1=IBASE
I2=I1+N*N
I3=I2+N*N
I4=I3+N
I5=I4+N
I6=I5+N
C
C PASS CALCULATED ADDRESS TO MAIN SUBROUTINE.
C
C---------------SPACE(I3) IS DUPLICATED TO MAKE BETA AND MN EQUIVALENT
C--------------- IN SUBR. MAINL
C--------------- SINCE ONE CANNOT USE EQUIVALENCE STATEMENTS ON
C--------------- ARGUMENTS PASSED TO A SUBR. THE SAME EXPLANATION
C--------------- APPLIES TO SPACE(I6)
CALL MAINL(N,SPACE(I1),SPACE(I2),SPACE(I3),SPACE(I3),
1SPACE(I4),SPACE(I5),SPACE(I6),SPACE(I6))
C
************************************************************************
C END OF ONE DATA SET, SEE IF MORE ARE TO BE PROCESSED
C***********************************************************************
C
GO TO 1
1000 STOP
END
C
C ..................................................................
C
C---------------N IS INPUT. OTHER ARGS. ARE SPACES RESERVED BY DYN. ALLOC.
SUBROUTINE MAINL(N,A,SXX,BETA,MN,X,XM,SX,MM)
C
C
DIMENSION A(1),SXX(1),BETA(1),X(1),XM(1),SX(1)
DIMENSION MM(1),MN(1),NOTF(48),ID(16)
EQUIVALENCE (BETAI,SXXNN,AII),(RESID,RESIDL,PC,SD),(Y,TB,F,VAR)
COMMON /IOBLK/IDLG,ICC,INP,IOUT,IO2,IO3,ICODE,IBNK,NAMI(2)
COMMON /BLOCK1/ISTD,NOTF,ID
C
C ..................................................................
C
55 FORMAT(' ROW',I4)
58 FORMAT('0OUTSIDE ALLOWABLE RANGE.')
63 FORMAT('0'//6X,4HMEAN,6X,'STANDARD DEVIATION',5X,8HVARIANCE)
67 FORMAT(1H ,2G,5X,G)
71 FORMAT(' REGRESSION CONSTANT =',G)
72 FORMAT('0THE B REGRESSION COEFFICIENTS')
80 FORMAT(I)
81 FORMAT(F)
84 FORMAT(' DEPENDENT VARIABLE IS NUMBER ? ',$)
90 FORMAT('0PARTIAL CORRELATION COEFFICIENTS')
92 FORMAT(' TABLE OF RESIDUALS? (YES OR NO)--',$)
94 FORMAT('0VARIABLE NUMBER',I3,' IS CONSTANT. JOB ABANDONED.')
97 FORMAT('0TOO FEW OBSERVATIONS. JOB TERMINATED.')
99 FORMAT('0NOTE:DEPENDENT VARIABLE IS NUMBER',I4,' WHICH APPEARS',
1I3,A2,' IN ANALYSIS'/)
112 FORMAT('+WITH ANY VARIABLES WHICH APPEARED AFTER DEPENDENT
1 SHIFTED ONE LOWER.'//)
503 FORMAT(7X,2HB(,I3,3H) =,G)
601 FORMAT(1X,8F9.5)
610 FORMAT('0MULTIPLE CORRELATION')
611 FORMAT('0F TEST OF B(1)=B(2)=...=B(N-1)=0.')
612 FORMAT(' F=',F12.5,' WITH (',F6.0,1H,F6.0,') DEGREES OF FREEDOM',
1 /,' F-PROB.=',F8.5)
613 FORMAT('0T TEST OF B(I)=0 FOR I=1,2,...,N-1')
618 FORMAT('0INTERCORRELATION')
619 FORMAT('0STANDARD ERROR OF ESTIMATE')
620 FORMAT(3(7X,G))
621 FORMAT(' TB (',I3,') =',F12.6,' WITH ',F6.0,' DEGREES OF FREEDOM',
1 /,' T-PROB.=',F8.5)
730 FORMAT(A3)
812 FORMAT(' ERROR READING OBSERVATION ',I11/)
813 FORMAT('+REENTER IT!'/)
9872 FORMAT(1H1,12X,1HY,16X,'PREDICTED Y',12X,8HRESIDUAL)
C
C ..................................................................
C
IDSK=1
C
C N IS NUMBER OF VARIABLES.
C NRESID: IS CODE FOR RESIDUALS.
C ND IS NUMBER OF THE DEPENDENT VARIABLE(IN ORIGINAL ORDER)
C
C
C READ NUMBER OF DEPENDENT VARIABLE.
C
40 WRITE (IDLG,84)
READ (ICC,80,END=1000) NDEPNT
ND=NDEPNT
IF(ND.LE.0)ND=N
IF(NDEPNT-N) 111,20,2585
2585 WRITE (IDLG,58)
IF (ICODE.GE.0) GO TO 40
CALL EXIT
20 NDEPNT=0
111 WRITE(IDLG,92)
READ(ICC,730) NRESID
IF (NRESID.EQ.'YES') CALL OFILE(IDSK,'TEMP')
IF (ISTD.NE.1) GO TO 100
NOTF(1)='(20F)'
DO 101 I=2,48
101 NOTF(I)=' '
100 IF (IO2.NE.'TTY') GO TO 103
WRITE(IDLG,102)
102 FORMAT(' ENTER DATA'/)
IF (ISTD.EQ.1) WRITE(IDLG,9999)
9999 FORMAT('+(AT MOST 20 NUMBERS PER LINE,SEPARATED BY COMMAS)'/)
GO TO 225
103 WRITE(IDLG,104)
104 FORMAT(' PLEASE WAIT, YOUR DATA IS BEING PROCESSED'/)
225 TOL1=.1
TOL2=1.
NN=N-1
NE=MOD(N,100)
NC=MOD(N,10)
NTH='TH'
IF((NC.EQ.1).AND.(NE.NE.11))NTH='ST'
IF((NC.EQ.2).AND.(NE.NE.12))NTH='ND'
IF((NC.EQ.3).AND.(NE.NE.13))NTH='RD'
C
C ZERO SX AND SXX
C
CALL ZEROH(SX,SXX,N,N)
T=0.
C
C READ DATA.
C
998 READ (INP,NOTF,ERR=850,END=3) (X(I),I=1,N)
2 T=T+1.
C
C MOVE DEPENDENT VARIABLE TO LAST.
C
IF(NDEPNT) 668,668,669
669 H=X(NDEPNT)
DO 671 J=NDEPNT,NN
671 X(J)=X(J+1)
X(N)=H
C
C FORM SUMS AND LOWER CORNER CROSS-PRODUCTS.
C
C---------------X, N INPUT; SX, SXX OUTPUT.
668 CALL XPRODH(X,SX,SXX,N,N)
IF(NRESID.EQ.'YES')WRITE (IDSK) (X(J),J=1,N)
GO TO 998
C
C END OF INPUT SECTION.
C
3 DO 53 I=2,N
I1=I-1
DO 53 J=1,I1
JI=I1*N+J
IJ=(J-1)*N+I
53 SXX(JI)=SXX(IJ)
WRITE(IOUT,105)
105 FORMAT(1H1)
DO 106 I=1,16
IF (ID(I).NE.' ') GO TO 107
106 CONTINUE
GO TO 160
107 WRITE(IOUT,108) ID
108 FORMAT(1X,16A5)
160 KOUNT=T
WRITE(IOUT,109) KOUNT,N
109 FORMAT('-NUMBER OF OBSERVATIONS =',I6/' NUMBER OF VARIABLES =',
1 I5//)
WRITE(IOUT, 99) ND,N, NTH
IF (ND.NE.N) WRITE(IOUT,112)
C
C CHECK FOR SUFFICIENT DATA.
C
IF (T .LT. 2.) GO TO 98
JJ=-N
DO 95 J=1,N
SXXNN=SX(J)
JJ=JJ+N
DO 5 I=1,J
IJ=JJ+I
A(IJ)=T*SXX(IJ)-SX(I)*SXXNN
JI=(I-1)*N+J
5 A(JI)=A(IJ)
JJJ=JJ+J
95 IF(A(JJJ).EQ.0.) GO TO 96
WRITE (IOUT,63)
II=-N
N1=N+1
JI=0
DO 64 I=1,N
II=II+N1
AII=A(II)
C
C COMPUTE CORRELATION MATRIX.
C
DO 600 J=1,N
JI=JI+1
JJ=(J-1)*N+J
600 SXX(JI)=A(JI)/SQRT(AII*A(JJ))
XM(I)=SX(I)/T
VAR=AII/(T*(T-1.))
SD=SQRT(VAR)
64 WRITE (IOUT,67) XM(I),SD,VAR
WRITE (IOUT,618)
N2=N*N
DO 616 I=1,N
WRITE (IOUT,55) I
WRITE(IOUT,601) (SXX(IJ),IJ=I,N2,N)
IJ=I-N
DO 617 J=1,N
IJ=IJ+N
617 A(IJ)=A(IJ)/T
616 X(I)=A(IJ)
C
C CHECK FOR SUFFICIENT DATA.
C
IF(T.LE.FLOAT(N)) GO TO 98
C
C MINVSQ FORMS INVERSE OF SQUARE MATRIX WITHIN ITSELF.
C
C---------------N=NO. OF VAR., NN=N-1 (ST. 225+2)
C---------------A, NN, TOL1, IOUT, N ARE INPUT. DET, IEXP ARE
C--------------- OUTPUT. SPACE PROVIDED FOR MM AND MN BY
C--------------- ST. 55-7, SUBR. MAINL ARG. LIST, AND CALL MAINL ARG.
C--------------- LIST IN MAIN PROG.
CALL MINVSQ (A,NN,TOL1,MM,MN,N,IOUT,2,DET,IEXP)
IF(DET.EQ.0) GO TO 776
C---------------INVERSE OF SXX USED IN ST. 825-1. INPUT SXX IS
C--------------- CORRELATION MATRIX (SEE ST. 64+1 AND 64+5)
CALL MINVSQ (SXX,N,TOL2,MM,MN,N,IOUT,2,DET,IEXP)
IF(DET.EQ.0) GO TO 776
C
C FORM REGRESSION COEFFICIENTS.
C
WRITE (IOUT,72)
BETAC=XM(N)
UNSQ=0.
DO 1001 I=1,NN
BETAI=0.
IJ=I-N
DO 501 J=1,NN
IJ=IJ+N
501 BETAI=BETAI+A(IJ)*X(J)
BETA(I)=BETAI
BETAC=BETAC-BETAI*XM(I)
UNSQ=UNSQ+BETAI*X(I)
1001 WRITE (IOUT,503) I,BETAI
47 WRITE (IOUT,71) BETAC
DEGES=NN
RESID=X(N)-UNSQ
F=N
DEGER=T-F
UMLRS=UNSQ/X(N)
S2=RESID/DEGER
S=SQRT(S2)
WRITE (IOUT,619)
WRITE (IOUT,620) S
IF(N.LT.3) GO TO 710
WRITE (IOUT,610)
UMLTR=SQRT(UMLRS)
WRITE (IOUT,601) UMLTR
710 WRITE (IOUT,611)
C
C COMPUTE F AND T VALUES.
C
F=(UMLRS*DEGER)/((1.-UMLRS)*DEGES)
FPRB=FISHER(INT(DEGES),INT(DEGER),F)
WRITE (IOUT,612) F,DEGES,DEGER,FPRB
WRITE (IOUT,613)
DO 614 I=1,NN
II=(I-1)*N+I
TB=BETA(I)/SQRT(S2*A(II))
TPRB=FISHER(1,INT(DEGER),TB*TB)
614 WRITE (IOUT,621) I,TB,DEGER,TPRB
IF(N.LT.3) GO TO 720
C
C COMPUTE PARTIAL CORRELATIONS.
C
WRITE (IOUT,90)
SXXNN=SXX(N2)
DO 825 I=1,NN
NI=I*N
II=(I-1)*N+I
PC=-SXX(NI)/SQRT(SXXNN*SXX(II))
825 WRITE (IOUT,601) PC
720 IF(NRESID.NE. 'YES') GO TO 776
C
C NOW FORM AND OUTPUT RESIDUALS.
C
ENDFILE IDSK
CALL IFILE(IDSK,'TEMP')
WRITE (IOUT,9872)
DO 9876 I=1,KOUNT
READ (IDSK,END=775) (X(J),J=1,N)
Y=BETAC
DO 9875 J=1,NN
9875 Y=Y+X(J)*BETA(J)
RESIDL=X(N)-Y
9876 WRITE (IOUT,620) X(N),Y,RESIDL
775 CALL RELEASE(IDSK)
CALL RENAMS(IDSK,1,'TEMP.DAT',0,0)
GO TO 776
96 WRITE (IOUT,94) J
GO TO 776
98 WRITE (IOUT,97)
776 WRITE(IDLG,110)
110 FORMAT('-')
RETURN
1000 STOP
850 KERR=T+1.1
WRITE(IDLG,812)KERR
IF (ICODE.LT.0) CALL EXIT
WRITE(IDLG,813)
GO TO 998
END