DIMENSION X(90,10),XPX(25,11),A(25,11),AA(275) DIMENSION XSUM(10),XMEAN(10),XPY(10),XOUT(10) DIMENSION Y(90),LABEL(25),LL(10) DIMENSION INA(11),INO(7) EQUIVALENCE (A,AA,XPX) 10 FORMAT(//21X,'MULTIPLE REGRESSION PROGRAM'//) 12 FORMAT(25X,5HYMEAN 9X,5HXMEAN ) C TYPE 9000 9000 FORMAT(' MULTIPLE REGRESSION PROGRAM') C LOOP=0 1 CALL MREG2(X,Y,N,NP,LOOP) LOOP=LOOP+1 FN = N 350 TYPE 10 TYPE 12 DO 30 I=1,N YSUM=YSUM+Y(I) DO 30 K=1,NP 30 XSUM(K)=XSUM(K)+X(I,K) YMEAN=YSUM/FN DO 35 I=1,NP 35 XMEAN(I)=XSUM(I)/FN TYPE 15,YMEAN,(XMEAN(I),I=1,NP) 15 FORMAT(20X,1P1E14.6/(34X,1P1E14.6)) C CONVERT DATA DO 45 I=1,N DO 40 J=1,NP 40 X(I,J)=X(I,J)-XMEAN(J) 45 Y(I) = Y(I)-YMEAN 50 TYPE 9050 9050 FORMAT(//1X,20(1H-)//' NUMBER OF VARIATES + TO CONSIDER? ',$) ACCEPT 9052,NVR IF(NVR)1,1,53 53 IF(NVR-NP)54,54,50 54 TYPE 9051 9051 FORMAT(/' WHICH ONES? ',$) 9052 FORMAT(11I) ACCEPT 9052,(LL(I),I=1,NVR) J=0 DO 55 I=1,NVR IF(LL(I))55,55,51 51 IF(LL(I)-NP)52,52,55 52 J=J+1 LL(J)=LL(I) 55 CONTINUE IF(J)50,50,56 56 NVR=J NR=NVR NC=NVR+1 NP1=NC DO 60 I=1,275 60 AA(I)=0. C GET X-PRODS + SUMS DO 72 II=1,NR I=LL(II) DO 72 JJ=II,NR J=LL(JJ) DO 70 K=1,N 70 XPX(II,JJ)=XPX(II,JJ)+X(K,I)*X(K,J) 72 XPX(JJ,II)=XPX(II,JJ) C RT SIDE OF NORMAL EQUATIONS DO 90 II=1,NR I=LL(II) DO 80 K=1,N 80 XPX(II,NP1)=XPX(II,NP1)+Y(K)*X(K,I) 90 XPY(II)=XPX(II,NP1) C SOLVE DO 21 J1=1,NR 21 LABEL(J1)=J1 DO 291 J1=1,NR 101 TEMP=0. DO 125 J2=J1,NR TEM=ABS(A(J2,J1))-TEMP IF (TEM) 121,122,122 122 TEMP=TEM IBIG=J2 121 IF (IBIG-J1) 5001,201,125 125 CONTINUE GO TO 201 131 DO 141 J2=1,NC TEMP=A(J1,J2) A(J1,J2)=A(IBIG,J2) 141 A(IBIG,J2)=TEMP I=LABEL(J1) LABEL(J1)=LABEL(IBIG) LABEL(IBIG)=I 201 TEMP=A(J1,J1) A(J1,J1)=1.0 DO 221 J2=1,NC 221 A(J1,J2)=A(J1,J2)/TEMP DO 281 J2=1,NR IF (J2-J1) 231, 281, 231 231 TEMP=A(J2,J1) A(J2,J1)=0. DO 241 J3=1,NC 241 A(J2,J3)=A(J2,J3)-TEMP*A(J1,J3) 281 CONTINUE 291 CONTINUE 301 N1=NR-1 DO 391 J1=1,N1 DO 321 J2=J1,NR IF (LABEL(J2)-J1) 321, 311, 321 311 IF (J2-J1) 5001, 391, 341 321 CONTINUE 341 DO 361 J3=1,NR TEMP=A(J3,J1) A(J3,J1)=A(J3,J2) 361 A(J3,J2)=TEMP LABEL(J2)=LABEL(J1) 391 CONTINUE 9903 SST=0. SSR=0. C SET UP ANOVA DO 100 I=1,N 100 SST=Y(I)*Y(I)+SST DO 102 I=1,NVR 102 SSR=SSR+XPX(I,NP1)*XPY(I) SSE=SST-SSR DFR=NVR DFT=N-1 IDFT=DFT DFE=DFT-DFR IDFE=DFE VSR=SSR/DFR VSE=SSE/DFE FTEST=VSR/VSE TYPE 9897 9897 FORMAT(' WANT TO SEE PREDICTED VALUES, TYPE YES OR + NO.--',$) ACCEPT 9102, IQ 9102 FORMAT(A3) IF(IQ-'YES')400,7400,400 7400 TYPE 2000 DO 130 I=1,N OUT=0. DO 110 KK=1,NVR K=LL(KK) 110 OUT=OUT+X(I,K)*XPX(KK,NP1) OUT=OUT+YMEAN YOUT=Y(I)+YMEAN 130 TYPE 3000,OUT,YOUT 2000 FORMAT(//24X,22HCALCULATED OBSERVED) 3000 FORMAT(20X,1P2E14.6) 460 FORMAT(5X,65(1H.)) 461 FORMAT(5X,1H.11X,25H.DEGREE OF FREE. SUM OF , +29HSQUARES. VARIANCE ESTIMATE . ) 400 TYPE 460 TYPE 461 TYPE 500,NVR,SSR,VSR 500 FORMAT(5X,13H. REGRESSION. I9,5X,1H., +1P1E14.6,3H .1P1E16.6,4H .) 501 FORMAT(5X,13H. REMAINDER .I9,5X,1H., +1P1E14.6,3H .1P1E16.6,4H .) TYPE 501,IDFE,SSE,VSE 502 FORMAT(5X,13H. TOTAL .I9,5X,1H.1P1E14.6,3H .19X,1 + H.) TYPE 502,IDFT,SST 465 FORMAT(7X,'LEAST SQUARE COEFFICIENTS',4X, +'COEFFICIENT T CONFIDENCE BAND') TYPE 465 DO 510 I=1,NVR SER=SQRT(A(I,I)*VSE) 510 TYPE 466,XPX(I,NP1),SER 466 FORMAT(2(13X,1PE14.6,4X)) 468 FORMAT(//' F RATIO(',I3,',',I3, + ' DEGREES OF FREEDOM) = ',1PE14.6) TYPE 468,NVR,IDFE,FTEST TYPE 9896 9896 FORMAT(/' VARIANCE-COVARIANCE MATRIX') DO 469 I=1,NVR 469 TYPE 1000,(XPX(I,J),J=I,NVR) 1000 FORMAT(1P2E14.6) GO TO 50 5001 CALL EXIT END SUBROUTINE MREG2(X,Y,N,NP,LOOP) DIMENSION X(90,10),Y(90) 1111 FORMAT(2I) 2222 FORMAT(11F) LERR=-2 IF(LOOP)11,11,101 11 TYPE 9001 9001 FORMAT(/' DO YOU DESIRE USER INSTRUCTIONS, + TYPE YES OR NO--',$) ACCEPT 9002,IREPLY 9002 FORMAT(A3) IF(IREPLY-'YES')101,201,101 201 TYPE 9201 TYPE 9202 TYPE 9203 TYPE 9204 TYPE 9205 TYPE 9206 TYPE 9207 TYPE 9288 9288 FORMAT(' IF DATA IS READ FROM A FILE, THE FIRST'/ + ' LINE MUST CONTAIN TWO INTEGER VARIABLES:'/ + ' NP, THE NUMBER OF DATA POINTS PER LINE; N, THE'/ + ' NUMBER OF OBSERVATIONS.') TYPE 9208 9201 FORMAT(/' MULTIPLE REGRESSION IS PRIMARILY AN +ATTEMPT TO CURVE FIT'/' OBSERVED DATA BY THE MODEL'/) 9202 FORMAT(10X,'Y-YMEAN = A(1)(X(1)-X1MEAN)+... ++A(P)(X(P)-XPMEAN))'/) 9203 FORMAT(' WHERE'/' YMEAN = MEAN OF OBSERVED + DATA'/' XIMEAN = MEAN OF I-TH INDEPENDENT + VARIABLE'/' THE A(I) ARE UNKNOWN COEFFICIENTS + DETERMINED BY THE'/' LEAST SQUARES PROCESS'/) 9204 FORMAT(' THE DATA IS ENTERED IN MATRIX FORM,'// +9X,'Y1 X11 X12 ... X1M'/ +9X,'Y2 X21 X22 ... X2M'/ +9X,'Y3 X31 X32 ... X3M'/ +9X,'.'/9X,'.'/9X,'.'/ +9X,'YN XN1 XN2 ... XNM'/) 9205 FORMAT(' Y1 IS THE OBSERVED DATA WHILE X11,X12, +X1M ARE THE CORRESPONDING'/' INDEPENDENT VARIABLES + FOR THE FIRST POINT. SIMILARLY UP TO THE'/ + ' N-TH POINT.'/) 9206 FORMAT(/' RESTRICTIONS:'/ +' M MUST NOT EXCEED 10'/' N MUST NOT EXCEED 90'/ +' INPUT FROM KEYBOARD OR DATA FILE MUST BE IN + UNIFORM LINE LENGTHS'/' WHICH MUST NOT EXCEED + 11 FIELDS EITHER IN THE FILE OR AT THE KEYBOARD.') 9207 FORMAT(' INPUT IS FREE FIELD IN THE SENSE THAT + COMMAS ARE INTERPRETED AS'/' SEPARATING + THE FIELDS.'/) 9208 FORMAT(' NOW YOU TRY IT.') C 101 TYPE 9101 9101 FORMAT(/' IS DATA TO BE READ FROM A FILE OR + ENTERED FROM THE KEYBOARD?'/' RESPOND WITH + "FILE", "KEY", OR DONE"--',$) 9102 FORMAT(A4) ACCEPT 9102,IREPLY IF(IREPLY-'FILE')102,301,102 102 IF(IREPLY-'KEY')103,401,103 103 CALL EXIT C 301 TYPE 9301 9301 FORMAT(/' WHAT IS FILE NAME?--',$) ACCEPT 9302,NAME 9302 FORMAT(A5) CALL IFILE(1,NAME) READ(1,1111)NP,N IF(NP)399,399,303 303 IF(N)399,399,305 305 IF(NP-10)307,307,399 307 IF(N-90)309,309,399 309 TYPE 9309,NP,N 9309 FORMAT(' FILE INDICATES',I3,' X VARIABLES AND', + I3,' OBSERVATIONS.'/) DO 330 I=1,N READ(1,2222)Y(I),(X(I,J),J=1,NP) IF(Y(I))330,320,330 320 IF(X(I,1))330,322,330 322 DO 324 J = 2,NP IF(X(I,J))330,324,330 324 CONTINUE N=I-1 TYPE 9309,NP,N GO TO 800 330 CONTINUE GO TO 800 399 TYPE 9309,NP,N TYPE 9399,NAME 9399 FORMAT(' NO MORE THAN 10 X"S OR MORE THAN 90 + OBSERVATIONS PERMITTED.'/' CORRECT DATA FILE ',A5/) CALL EXIT C 401 TYPE 9401 9401 FORMAT(/' DATA WILL BE ACCEPTED FROM KEYBOARD.'/ + ' GIVE NUMBER OF OBSERVATIONS (N) AND NUMBER OF + VARIABLES (M)--',$) 9402 FORMAT(2I) ACCEPT 9402,N,NP IF(NP)499,499,403 403 IF(N)499,499,405 405 IF(NP-10)407,407,499 407 IF(N-90)409,409,499 409 TYPE 9409,NP 9409 FORMAT(/' ON EACH NUMBERED LINE ENTER FIRST THE + VALUE OF THE DEPENDENT (Y)'/' VARIABLE. THEN + IN SUCCESSION ENTER THE VALUES OF THE',I3, + ' INDEPENDENT'/' VARIABLES. SEPARATE EACH VALUE BY + A COMMA.'/) 9410 FORMAT('+ LINE',I3,2X,$) 9420 FORMAT(11F) I=0 DO 430 K=1,N I=I+1 TYPE 9410,I ACCEPT 9420,Y(I),(X(I,J),J=1,NP) IF(Y(I))430,420,430 420 IF(X(I,1))430,422,430 422 DO 424 J=2,NP IF(X(I,J))430,424,430 424 CONTINUE I=I-1 TYPE 9424 9424 FORMAT(/' IF YOU WANT TO ENTER MORE DATA, TYPE + "1". IF ALL DATA HAS BEEN'/' ENTERED, TYPE "0". WHAT + DO YOU WANT?--',$) ACCEPT 9402,J IF(J)426,426,430 426 N=I GO TO 440 430 CONTINUE 440 TYPE 9440,N 444 TYPE 9444 ACCEPT 9002,IREPLY IF(IREPLY-'YES')800,450,800 9444 FORMAT(/' ARE THERE ANY CORRECTIONS YOU WANT TO + MAKE? ',$) 9440 FORMAT(/' YOU HAVE ENTERED',I3,' OBSERVATIONS.'/) 499 TYPE 9499,NP,N 9499 FORMAT(/' YOU GAVE M AS',I3,', BUT IT MUST BE + BETWEEN 1 AND 10.'/' YOU GAVE N AS',I3,', BUT IT + MUST BE BETWEEN M AND 90.'/) LERR=LERR+1 IF(LERR)401,401,101 800 RETURN 450 TYPE 9450 9450 FORMAT(' ENTER LINE NUMBER, THEN VARIABLE NUMBER, + THEN THE CORRECT VALUE.'/' IF Y IS TO BE + CORRECTED, THEN ENTER "0" FOR'/ + ' VARIABLE NUMBER.'/) 460 TYPE 9460 9460 FORMAT('+ WHICH LINE, NUMBER, VALUE? ',$) 9461 FORMAT(2I,F) ACCEPT 9461,LINE,NUMBER,VALUE IF(LINE)444,444,471 471 IF(NUMBER)444,472,472 472 IF(LINE-N)473,473,444 473 IF(NUMBER-NP)474,474,444 474 IF(NUMBER)475,475,476 475 Y(LINE)=VALUE GO TO 444 476 X(LINE,NUMBER)=VALUE GO TO 444 END