Google
 

Trailing-Edge - PDP-10 Archives - decus_20tap5_198111 - decus/20-0137/pathan/pathan.for
There are 2 other files named pathan.for in the archive. Click here to see a list.
C	WESTERN MICHIGAN UNIVERSITY
C	PATHAN.FOR (FILE NAME ON LIBRARY DECTAPE)
C	PATHAN, 1.8.1 (CALLING NAME, SUBLST #)
C	PATH ANALYSIS
C	THIS PROGRAM WAS WRITTEN BY SAM ANEMA AND LATER MODIFIED
C	 BY RUSS R. BARR, DAVE SCHULZ
C	APLIB PROGRAMS USED:  IOB, GETFOR
C	LIBRARY DECTAPE PROGRAMS USED:  USAGE.MAC
C	FORWMU PROGRAMS USED:  MINVSQ, DEVICE, DEVCHG, EXISTS, PRINTS
C	INTERNAL SUBR. USED:  LIST, REGRES
C	ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
	DIMENSION IDN(16),FMT(48),ISD(9),INP(80),IW(6)
C	DEFINED BY:(NVARS)
C---------------50 IN DIMENSIONS CORRESPONDS TO LIMITATION OF
C--------------- 50 VARS. IN WRITE UP.
C---------------IVN IS VECTOR OF LABELS FOR VARS.
	DIMENSION IVN(50),NIND(50),NUM(50)
C	DEFINED BY:(NVARS,NVARS)
	DIMENSION RS(50),RESD(50),COR(50,50),BMAT(50,50),SMAT(50,50)
C	DEFINED BY:(NLEVEL)
	DIMENSION NOR(20)
C	DEFINED BY:(NPLEVL)
	DIMENSION ILS(20),SE(20),BETA(20)
C	DEFINED BY:(NPLEVL+2)
	DOUBLE PRECISION FMA(22)
C	DEFINED BY:(NLEVEL,NPLEVL)
	DIMENSION NORD(20,20)
	COMMON /IOBLK/IDLG,INT,INPUT,IREP,IDEV,IDEVA,ICODE,IC,NAMI(2)
	COMMON /IOBLKA/NAMO(2),IPJ,IPG,NCOPYS,IAUX
	COMMON /VARTMP/IDUM(72),ISAVE(5)
	DOUBLE PRECISION FMA1,FMA2,FMA3,FMA4
	DATA FMA1,FMA2,FMA3,FMA4/'(1X,A5,F7.',',2(7X,A1)',
     #',2F8.4    ',')         '/
	DATA ISD/'1','2','3','4','5','6','7','8','9'/
C
C		INITIALIZE
C
	INT=-4
	IDLG=-1
	IREP=2
	INPUT=3
	NVARS=50
C---------------SEE FIG. 1 OF WRITE UP.  NLEVEL CORRESPONDS
C--------------- TO NO. OF COLS. NPLEVEL CORRESPONDS TO NO. OF
C--------------- VARS. IN A COL.
	NLEVEL=20
	NPLEVL=20
	DO 1470 I=1,NVARS
	RS(I)=0.
1470	RESD(I)=0.
	WRITE(IDLG,10)
10	FORMAT('0---WMU-CSR PATH ANALYSIS PROGRAM'//)
C
C		ENTER PRELIMINARY DOPE
C
C	CALL USAGEB('PATHAN')
	CALL IOB(1)
2	CALL IOB(0)
	NV=NVARS
	CALL GETLAB(NV,IVN,NUM)
	WRITE(IDLG,921)
921	FORMAT(' SAMPLE SIZE? ',$)
	READ(INT,922)XN
922	FORMAT(F)
	CALL GETFOR(IDLG,INT,FMT,ISTD,48,2)
	IF(ISTD.EQ.1)FMT(1)='(10F)'
	IF(IDEV.EQ.'TTY')WRITE(IDLG,855)
855	FORMAT(' ENTER CORRELATION MATRIX',/)
	DO 40 I=1,NV
C
C		READ CORRELATION MATRIX GENERATED BY CORL
C
	READ(INPUT,FMT)(COR(I,J),J=1,NV)
40	CONTINUE
C
C		MORE PARAMETERS ETC.
C
510	WRITE(IDLG,51)
51	FORMAT(' ENTER IDENTIFICATION')
	READ(INT,52)IDN
52	FORMAT(16A5)
	DO 502 I=15,1,-1
	IF(IDN(I).NE.' ') GOTO 503
502	CONTINUE
503	IDN(16)=I
	WRITE(IDLG,60)
60	FORMAT(' ENTER ORDERING OF VARIABLES'/)
	I=0
69	I=I+1
C---------------NL IS NO. OF VARS. SUBMITTED BY USER IN A LEVEL.
C--------------- SEE ST. 1470-3, NV=NO. OF VARS.
	CALL GETLST(NV,NL,IVN,ILS)
	IF(NL.EQ.0) GOTO 200
C---------------SEE COMMENT FOR ST. 1470-4 ABOVE.
	IF(I.GT.NLEVEL)GO TO 70
	IF(NL.GT.NPLEVL)GO TO 72
C---------------STORE ITH LEVEL (COL.) SEQ. VAR. ID INTO THE ITH
C--------------- ROW OF NORD(I,J) AND STORE THE NO. OF VARS. IN
C--------------- THE ITH COLUMN INTO NOR(I).
	DO 90 J=1,NL
90	NORD(I,J)=ILS(J)
	NOR(I)=NL
	GO TO 69
70	WRITE(IDLG,71)NLEVEL
71	FORMAT(' ?MAXIMUM OF ',I3,' LEVELS EXCEEDED',/)
	IF(ICODE.NE.0) CALL EXIT
	GO TO 2
C---------------COME HERE FROM ST. 90-2.
72	WRITE(IDLG,73)NPLEVL
73	FORMAT(' ?MAXIMUM OF ',I3,' VARIABLES PER LEVEL EXCEEDED',/)
	IF(ICODE.NE.0) CALL EXIT
	GO TO 2
C---------------COME HERE FROM ST. 69+2.  NO.=NO. OF LEVELS (COLS)
200	NO=I-1
110	WRITE(IDLG,401)
401	FORMAT(' MANUAL OR AUTOMATIC? '$)
	READ(INT,402)IHOW
402	FORMAT(A1)
	DO 804 I=1,NLEVEL
	DO 804 J=1,NPLEVL
	BMAT(I,J)=0
804	SMAT(I,J)=0
	IF(IHOW.EQ.'A')GO TO 501
	IF(IHOW.NE.'M')GO TO 110
C
C		MANUAL OPERATION
C
	DO 403 I=1,NO
	DO 404 J=1,NOR(I)
420	WRITE(IDLG,405)IVN(NORD(I,J))
405	FORMAT(' WHICH ARE THE VARIABLES CAUSING VARIABLE ',A5/)
398	CALL GETLST(NV,NL,IVN,ILS)
	IF(NL.EQ.0) GOTO 404
	DO 407 K=1,NL
407	NIND(K)=ILS(K)
	NI=NL
	ND=NORD(I,J)
C
C		PERFORM REGRESSION
C
C---------------COR, NIND, NI, ND, XN, ARE INPUT.  BETA, RSQ,
C--------------- SE ARE RETURNED BY SUBR. REGRES.  XN IS READ
C--------------- IN ST. 921+1,  NIND, NI, ND, ARE IN
C--------------- 3 STATEMENTS JUST ABOVE.
	CALL REGRES(COR,NIND,NI,ND,BETA,RSQ,SE,XN)
	R=SQRT(1.0-RSQ)
	RS(ND)=RSQ
	RESD(ND)=R
	WRITE(IDLG,270)IVN(ND),RSQ,R,(IVN(NIND(J3)),BETA(J3),SE(J3),
     #J3=1,NI)
270	FORMAT('0DEP = ',A5/' RSQ = ',F7.4,/' RESIDUAL PATH = ',F7.4/
	1'0IND      BETA     S.E.'/
	1(1X,A5,1XF8.4,1X,F8.4))
	IF(XN-NI.LT.2)WRITE(IDLG,887)
887	FORMAT(' **S.E. NOT CALCULATED BECAUSE NUMBER OF SAMPLES',
     #'  MINUS',/,' NUMBER OF INDEPENDENT VARIABLES IS LESS',
     #' THAN 2.')
422	WRITE(IDLG,409)
409	FORMAT(' OK? ',$)
	READ(INT,402)IYN
	IF(IYN.EQ.'N')GO TO 420
	IF(IYN.NE.'Y')GO TO 422
C---------------SAVE BETA AND SE FOR SUMMARY MATRIX PRINTOUT
C--------------- SEE ST. 880.
	DO 801	J3=1,NI
	BMAT(ND,NIND(J3))=BETA(J3)
801	SMAT(ND,NIND(J3))=SE(J3)
	DO 600 I3=1,NI
600	CONTINUE
404	CONTINUE
403	CONTINUE
C
77	CONTINUE
	GO TO 881
C		AUTOMATIC OPERATION
C
501	WRITE(IDLG,310)
310	FORMAT(' CRITICAL BETA VALUE? ',$)
	READ(INT,320)CB
320	FORMAT(F)
C---------------NO=NO. OF LEVELS (COLS).  NOR(I) CONTAINS NO.
C--------------- OF VARS. IN COL I
	DO 210 I=1,NO
	DO 220 J=1,NOR(I)
C
C		SELECT DEPENDENT VARIABLE
C
	ND=NORD(I,J)
	NI=0
C
C		FIND ALL POSSIBLE INDEPENDENT VARIABLES
C
	DO 230 I1=1,I
	DO 240 J1=1,NOR(I1)
	IF(I.EQ.I1.AND.J.EQ.J1)GO TO 240
	NI=NI+1
	NIND(NI)=NORD(I1,J1)
240	CONTINUE
230	CONTINUE
260	IF(NI.EQ.0)GO TO 220
C
C		PERFORM REGRESSION
C
	CALL REGRES(COR,NIND,NI,ND,BETA,RSQ,SE,XN)
C
C		FIND SMALLEST BETA WEIGHT
C
	SMA=1.E30
	DO 250 I1=1,NI
	IF(ABS(BETA(I1)).GE.SMA)GO TO 250
	SMA=ABS(BETA(I1))
	K=I1
250	CONTINUE
C
C			IF THE SMALLEST BETA WEIGHT IS LARGER THAN THE
C			CRITICAL BETA VALUE DO NOT DELETE ANY MORE 
C			VARIABLES.
C
	IF(SMA.GE.CB)GO TO 276
	IF(K.EQ.NI)GO TO 275
C
C		DELETE VARIABLE WITH THE SMALLEST BETA WEIGHT
C
	DO 271 I1=K,NI-1
271	NIND(I1)=NIND(I1+1)
275	NI=NI-1
	GO TO 260
C
276	R=SQRT(1.0-RSQ)
	RS(ND)=RSQ
	RESD(ND)=R
	WRITE(IDLG,270)IVN(ND),RSQ,R,(IVN(NIND(J3)),BETA(J3),SE(J3),
     #J3=1,NI)
	IF(XN-NI.LT.2)WRITE(IDLG,887)
C---------------STORE BETA AND SE FOR SUMMARY MATRIX
	DO 803	J3=1,NI
	BMAT(ND,NIND(J3))=BETA(J3)
803	SMAT(ND,NIND(J3))=SE(J3)
220	CONTINUE
210	CONTINUE
881	NWIDE=7
	IF(IDEVA.EQ.'TTY')NWIDE=3
	WRITE(IREP,880)(IDN(I),I=1,IDN(16))
880	FORMAT(///,' ***SUMMARY MATRIX***',/1X,16A5)
	NVA=1
	BLK=' '
	FMA(1)=FMA1
	FMA(2)='4,2X,F7.4'
882	IF(NVA.GT.NV)GO TO 789
	NVB=MIN0(NVA+NWIDE-1,NV)
	WRITE(IREP,883)(BLK,IVN(J4),J4=NVA,NVB)
883	FORMAT(///'     INP:',17X,7(A1,3X,A5,7X))
	WRITE(IREP,884)(BLK,J4=NVA,NVB)
884	FORMAT(6X,'  RSQ     RESID  ',7(A1,2X,'BETA',4X,'S.E.',1X))
	WRITE(IREP,888)
888	FORMAT(' DEP')
	DO 885 J3=1,NV
	J5=2
	DO 886	J4=NVA,NVB
	J5=J5+1
	FMA(J5)=FMA2
886	IF(BMAT(J3,J4).NE.0)FMA(J5)=FMA3
	FMA(J5+1)=FMA4
	WRITE(IREP,FMA)IVN(J3),RS(J3),RESD(J3),(BMAT(J3,J4),SMAT(J3,J4),
     #J4=NVA,NVB)
885	CONTINUE
	NVA=NVA+NWIDE
	GO TO 882
789	WRITE(IDLG,788)
788	FORMAT(/)
	GO TO 2
	END
C---------------SEE REFS. 1 AND 2 IN WRITE UP.  ALSO SEE
C--------------- MAIN PROG. ST. 270-4 FOR ARGS. LIST.
	SUBROUTINE REGRES(COR,NIND,NI,ND,BETA,RSQ,SE,XN)
	DIMENSION COR(50,50),NIND(1),BETA(1),A(50,50),XTRA(75),SE(1)
	DO 10 I=1,NI
	DO 20 J=1,NI
	A(I,J)=COR(NIND(I),NIND(J))
20	CONTINUE
10	CONTINUE
C		INVERT MATRIX
C
C---------------MINVSQ IS A MACRO SUBR. IN FORWMU.  A, NI ARE
C--------------- INPUT.  BETA, XTRA, DET, IEXP ARE RETURNED.  A IS
C--------------- RETURNED AS INVERSE OF INPUT A.
	CALL MINVSQ(A,NI,1.,BETA,XTRA,50,5,1,DET,IEXP)
	RSQ=0.0
	DO 30 I=1,NI
	BETA(I)=0.
	SE(I)=0
	DO 40 J=1,NI
C
C		CALCULATE BETAS
C
	BETA(I)=BETA(I)+COR(ND,NIND(J))*A(I,J)
40	CONTINUE
C
C		CALCULATE R-SQUARE
C
	RSQ=RSQ+BETA(I)*COR(ND,NIND(I))
30	CONTINUE
	IF((XN-NI-1).EQ.0)RETURN
	DO 100 I=1,NI
	SE(I)=SQRT((1.0-RSQ)*A(I,I)/(XN-NI-1))
100	CONTINUE
	RETURN
	END