Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-05 - 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