Trailing-Edge
-
PDP-10 Archives
-
decuslib10-01
-
43,50144/mreg.f4
There are no other files named mreg.f4 in the archive.
00100 DIMENSION X(90,10),XPX(25,11),A(25,11),AA(275)
00110 DIMENSION XSUM(10),XMEAN(10),XPY(10),XOUT(10)
00120 DIMENSION Y(90),LABEL(25),LL(10)
00130 DIMENSION INA(11),INO(7)
00140 EQUIVALENCE (A,AA,XPX)
00150 10 FORMAT(//21X,'MULTIPLE REGRESSION PROGRAM'//)
00160 12 FORMAT(25X,5HYMEAN 9X,5HXMEAN )
00170 C
00180 TYPE 9000
00190 9000 FORMAT(' MULTIPLE REGRESSION PROGRAM')
00200 C
00210 LOOP=0
00220 1 CALL MREG2(X,Y,N,NP,LOOP)
00230 LOOP=LOOP+1
00240 FN = N
00250 350 TYPE 10
00260 TYPE 12
00270 DO 30 I=1,N
00280 YSUM=YSUM+Y(I)
00290 DO 30 K=1,NP
00300 30 XSUM(K)=XSUM(K)+X(I,K)
00310 YMEAN=YSUM/FN
00320 DO 35 I=1,NP
00330 35 XMEAN(I)=XSUM(I)/FN
00340 TYPE 15,YMEAN,(XMEAN(I),I=1,NP)
00350 15 FORMAT(20X,1P1E14.6/(34X,1P1E14.6))
00360 C CONVERT DATA
00370 DO 45 I=1,N
00380 DO 40 J=1,NP
00390 40 X(I,J)=X(I,J)-XMEAN(J)
00400 45 Y(I) = Y(I)-YMEAN
00410 50 TYPE 9050
00420 9050 FORMAT(//1X,20(1H-)//' NUMBER OF VARIATES
00430 + TO CONSIDER? ',$)
00440, ACCEPT 9052,NVR
00450 IF(NVR)1,1,53
00460 53 IF(NVR-NP)54,54,50
00470 54 TYPE 9051
00480 9051 FORMAT(/' WHICH ONES? ',$)
00490 9052 FORMAT(11I)
00500 ACCEPT 9052,(LL(I),I=1,NVR)
00510 J=0
00520 DO 55 I=1,NVR
00530 IF(LL(I))55,55,51
00540 51 IF(LL(I)-NP)52,52,55
00550 52 J=J+1
00560 LL(J)=LL(I)
00570 55 CONTINUE
00580 IF(J)50,50,56
00590 56 NVR=J
00600 NR=NVR
00610 NC=NVR+1
00620 NP1=NC
00630 DO 60 I=1,275
00640 60 AA(I)=0.
00650 C GET X-PRODS + SUMS
00660 DO 72 II=1,NR
00670 I=LL(II)
00680 DO 72 JJ=II,NR
00690 J=LL(JJ)
00700 DO 70 K=1,N
00710 70 XPX(II,JJ)=XPX(II,JJ)+X(K,I)*X(K,J)
00720 72 XPX(JJ,II)=XPX(II,JJ)
00730 C RT SIDE OF NORMAL EQUATIONS
00740 DO 90 II=1,NR
00750 I=LL(II)
00760 DO 80 K=1,N
00770 80 XPX(II,NP1)=XPX(II,NP1)+Y(K)*X(K,I)
00780 90 XPY(II)=XPX(II,NP1)
00790 C SOLVE
00800 DO 21 J1=1,NR
00810 21 LABEL(J1)=J1
00820 DO 291 J1=1,NR
00830 101 TEMP=0.
00840 DO 125 J2=J1,NR
00850 TEM=ABS(A(J2,J1))-TEMP
00860 IF (TEM) 121,122,122
00870 122 TEMP=TEM
00880 IBIG=J2
00890 121 IF (IBIG-J1) 5001,201,125
00900 125 CONTINUE
00910 GO TO 201
00920 131 DO 141 J2=1,NC
00930 TEMP=A(J1,J2)
00940 A(J1,J2)=A(IBIG,J2)
00950 141 A(IBIG,J2)=TEMP
00960 I=LABEL(J1)
00970 LABEL(J1)=LABEL(IBIG)
00980 LABEL(IBIG)=I
00990 201 TEMP=A(J1,J1)
01000 A(J1,J1)=1.0
01010 DO 221 J2=1,NC
01020 221 A(J1,J2)=A(J1,J2)/TEMP
01030 DO 281 J2=1,NR
01040 IF (J2-J1) 231, 281, 231
01050 231 TEMP=A(J2,J1)
01060 A(J2,J1)=0.
01070 DO 241 J3=1,NC
01080 241 A(J2,J3)=A(J2,J3)-TEMP*A(J1,J3)
01090 281 CONTINUE
01100 291 CONTINUE
01110 301 N1=NR-1
01120 DO 391 J1=1,N1
01130 DO 321 J2=J1,NR
01140 IF (LABEL(J2)-J1) 321, 311, 321
01150 311 IF (J2-J1) 5001, 391, 341
01160 321 CONTINUE
01170 341 DO 361 J3=1,NR
01180 TEMP=A(J3,J1)
01190 A(J3,J1)=A(J3,J2)
01200 361 A(J3,J2)=TEMP
01210 LABEL(J2)=LABEL(J1)
01220 391 CONTINUE
01230 9903 SST=0.
01240 SSR=0.
01250 C SET UP ANOVA
01260 DO 100 I=1,N
01270 100 SST=Y(I)*Y(I)+SST
01280 DO 102 I=1,NVR
01290 102 SSR=SSR+XPX(I,NP1)*XPY(I)
01300 SSE=SST-SSR
01310 DFR=NVR
01320 DFT=N-1
01330 IDFT=DFT
01340 DFE=DFT-DFR
01350 IDFE=DFE
01360 VSR=SSR/DFR
01370 VSE=SSE/DFE
01380 FTEST=VSR/VSE
01390 TYPE 9897
01400 9897 FORMAT(' WANT TO SEE PREDICTED VALUES, TYPE YES OR
01410 + NO.--',$)
01420 ACCEPT 9102, IQ
01430 9102 FORMAT(A3)
01440 IF(IQ-'YES')400,7400,400
01450 7400 TYPE 2000
01460 DO 130 I=1,N
01470 OUT=0.
01480 DO 110 KK=1,NVR
01490 K=LL(KK)
01500 110 OUT=OUT+X(I,K)*XPX(KK,NP1)
01510 OUT=OUT+YMEAN
01520 YOUT=Y(I)+YMEAN
01530 130 TYPE 3000,OUT,YOUT
01540 2000 FORMAT(//24X,22HCALCULATED OBSERVED)
01550 3000 FORMAT(20X,1P2E14.6)
01560 460 FORMAT(5X,65(1H.))
01570 461 FORMAT(5X,1H.11X,25H.DEGREE OF FREE. SUM OF ,
01580 +29HSQUARES. VARIANCE ESTIMATE . )
01590 400 TYPE 460
01600 TYPE 461
01610 TYPE 500,NVR,SSR,VSR
01620 500 FORMAT(5X,13H. REGRESSION. I9,5X,1H.,
01630 +1P1E14.6,3H .1P1E16.6,4H .)
01640 501 FORMAT(5X,13H. REMAINDER .I9,5X,1H.,
01650 +1P1E14.6,3H .1P1E16.6,4H .)
01660 TYPE 501,IDFE,SSE,VSE
01670 502 FORMAT(5X,13H. TOTAL .I9,5X,1H.1P1E14.6,3H .19X,1
01680 + H.)
01690 TYPE 502,IDFT,SST
01700 465 FORMAT(7X,'LEAST SQUARE COEFFICIENTS',4X,
01710 +'COEFFICIENT T CONFIDENCE BAND')
01720 TYPE 465
01730 DO 510 I=1,NVR
01740 SER=SQRT(A(I,I)*VSE)
01750 510 TYPE 466,XPX(I,NP1),SER
01760 466 FORMAT(2(13X,1PE14.6,4X))
01770 468 FORMAT(//' F RATIO(',I3,',',I3,
01780 + ' DEGREES OF FREEDOM) = ',1PE14.6)
01790 TYPE 468,NVR,IDFE,FTEST
01800 TYPE 9896
01810 9896 FORMAT(/' VARIANCE-COVARIANCE MATRIX')
01820 DO 469 I=1,NVR
01830 469 TYPE 1000,(XPX(I,J),J=I,NVR)
01840 1000 FORMAT(1P2E14.6)
01850 GO TO 50
01860 5001 CALL EXIT
01870 END
01000 SUBROUTINE MREG2(X,Y,N,NP,LOOP)
01010 DIMENSION X(90,10),Y(90)
01020 1111 FORMAT(2I)
01030 2222 FORMAT(11F)
01040 LERR=-2
01050 IF(LOOP)11,11,101
01060 11 TYPE 9001
01070 9001 FORMAT(/' DO YOU DESIRE USER INSTRUCTIONS,
01080 + TYPE YES OR NO--',$)
01090 ACCEPT 9002,IREPLY
01100 9002 FORMAT(A3)
01110 IF(IREPLY-'YES')101,201,101
01120 201 TYPE 9201
01130 TYPE 9202
01140 TYPE 9203
01150 TYPE 9204
01160 TYPE 9205
01170 TYPE 9206
01180 TYPE 9207
01185 TYPE 9288
01186 9288 FORMAT(' IF DATA IS READ FROM A FILE, THE FIRST'/
01187 + ' LINE MUST CONTAIN TWO INTEGER VARIABLES:'/
01188 + ' NP, THE NUMBER OF DATA POINTS PER LINE; N, THE'/
01189 + ' NUMBER OF OBSERVATIONS.')
01190 TYPE 9208
01200 9201 FORMAT(/' MULTIPLE REGRESSION IS PRIMARILY AN
01210 +ATTEMPT TO CURVE FIT'/' OBSERVED DATA BY THE MODEL'/)
01220 9202 FORMAT(10X,'Y-YMEAN = A(1)(X(1)-X1MEAN)+...
01230 ++A(P)(X(P)-XPMEAN))'/)
01240 9203 FORMAT(' WHERE'/' YMEAN = MEAN OF OBSERVED
01250 + DATA'/' XIMEAN = MEAN OF I-TH INDEPENDENT
01260 + VARIABLE'/' THE A(I) ARE UNKNOWN COEFFICIENTS
01270 + DETERMINED BY THE'/' LEAST SQUARES PROCESS'/)
01280 9204 FORMAT(' THE DATA IS ENTERED IN MATRIX FORM,'//
01290 +9X,'Y1 X11 X12 ... X1M'/
01300 +9X,'Y2 X21 X22 ... X2M'/
01310 +9X,'Y3 X31 X32 ... X3M'/
01320 +9X,'.'/9X,'.'/9X,'.'/
01330 +9X,'YN XN1 XN2 ... XNM'/)
01340 9205 FORMAT(' Y1 IS THE OBSERVED DATA WHILE X11,X12,
01350 +X1M ARE THE CORRESPONDING'/' INDEPENDENT VARIABLES
01360 + FOR THE FIRST POINT. SIMILARLY UP TO THE'/
01370 + ' N-TH POINT.'/)
01380 9206 FORMAT(/' RESTRICTIONS:'/
01390 +' M MUST NOT EXCEED 10'/' N MUST NOT EXCEED 90'/
01400 +' INPUT FROM KEYBOARD OR DATA FILE MUST BE IN
01410 + UNIFORM LINE LENGTHS'/' WHICH MUST NOT EXCEED
01420 + 11 FIELDS EITHER IN THE FILE OR AT THE KEYBOARD.')
01430 9207 FORMAT(' INPUT IS FREE FIELD IN THE SENSE THAT
01440 + COMMAS ARE INTERPRETED AS'/' SEPARATING
01450 + THE FIELDS.'/)
01460 9208 FORMAT(' NOW YOU TRY IT.')
01470 C
01480 101 TYPE 9101
01490 9101 FORMAT(/' IS DATA TO BE READ FROM A FILE OR
01500 + ENTERED FROM THE KEYBOARD?'/' RESPOND WITH
01510 + "FILE", "KEY", OR DONE"--', $)
01520 9102 FORMAT(A4)
01530 ACCEPT 9102,IREPLY
01540 IF(IREPLY-'FILE')102,301,102
01550 102 IF(IREPLY-'KEY')103,401,103
01560 103 CALL EXIT
01570 C
01580 301 TYPE 9301
01590 9301 FORMAT(/' WHAT IS FILE NAME?--',$)
01600 ACCEPT 9302,NAME
01610 9302 FORMAT(A5)
01620 CALL IFILE(1,NAME)
01630 READ(1,1111)NP,N
01640 IF(NP)399,399,303
01650 303 IF(N)399,399,305
01660 305 IF(NP-10)307,307,399
01670 307 IF(N-90)309,309,399
01680 309 TYPE 9309,NP,N
01690 9309 FORMAT(' FILE INDICATES',I3,' X VARIABLES AND',
01700 + I3,' OBSERVATIONS.'/)
01710 DO 330 I=1,N
01720 READ(1,2222)Y(I),(X(I,J),J=1,NP)
01730 IF(Y(I))330,320,330
01740 320 IF(X(I,1))330,322,330
01750 322 DO 324 J = 2,NP
01760 IF(X(I,J))330,324,330
01770 324 CONTINUE
01780 N=I-1
01790 TYPE 9309,NP,N
01800 GO TO 800
01810 330 CONTINUE
01820 GO TO 800
01830 399 TYPE 9309,NP,N
01840 TYPE 9399,NAME
01850 9399 FORMAT(' NO MORE THAN 10 X"S OR MORE THAN 90
01860 + OBSERVATIONS PERMITTED.'/' CORRECT DATA FILE ',A5/)
01870 CALL EXIT
01880 C
01890 401 TYPE 9401
01900 9401 FORMAT(/' DATA WILL BE ACCEPTED FROM KEYBOARD.'/
01910 + ' GIVE NUMBER OF OBSERVATIONS (N) AND NUMBER OF
01920 + VARIABLES (M)--',$)
01930 9402 FORMAT(2I)
01940 ACCEPT 9402,N,NP
01950 IF(NP)499,499,403
01960 403 IF(N)499,499,405
01970 405 IF(NP-10)407,407,499
01980 407 IF(N-90)409,409,499
01990 409 TYPE 9409,NP
02000 9409 FORMAT(/' ON EACH NUMBERED LINE ENTER FIRST THE
02010 + VALUE OF THE DEPENDENT (Y)'/' VARIABLE. THEN
02020 + IN SUCCESSION ENTER THE VALUES OF THE',I3,
02030 + ' INDEPENDENT'/' VARIABLES. SEPARATE EACH VALUE BY
02040 + A COMMA.'/)
02050 9410 FORMAT('+ LINE',I3,2X,$)
02060 9420 FORMAT(11F)
02070 I=0
02080 DO 430 K=1,N
02090 I=I+1
02100 TYPE 9410,I
02110 ACCEPT 9420,Y(I),(X(I,J),J=1,NP)
02120 IF(Y(I))430,420,430
02130 420 IF(X(I,1))430,422,430
02140 422 DO 424 J=2,NP
02150 IF(X(I,J))430,424,430
02160 424 CONTINUE
02170 I=I-1
02180 TYPE 9424
02190 9424 FORMAT(/' IF YOU WANT TO ENTER MORE DATA, TYPE
02200 + "1". IF ALL DATA HAS BEEN'/' ENTERED, TYPE "0". WHAT
02210 + DO YOU WANT?--',$)
02220 ACCEPT 9402,J
02230 IF(J)426,426,430
02240 426 N=I
02250 GO TO 440
02260 430 CONTINUE
02270 440 TYPE 9440,N
02280 444 TYPE 9444
02290 ACCEPT 9002,IREPLY
02300 IF(IREPLY-'YES')800,450,800
02310 9444 FORMAT(/' ARE THERE ANY CORRECTIONS YOU WANT TO
02320 + MAKE? ',$)
02330 9440 FORMAT(/' YOU HAVE ENTERED',I3,' OBSERVATIONS.'/)
02340 499 TYPE 9499,NP,N
02350 9499 FORMAT(/' YOU GAVE M AS',I3,', BUT IT MUST BE
02360 + BETWEEN 1 AND 10.'/' YOU GAVE N AS',I3,', BUT IT
02370 + MUST BE BETWEEN M AND 90.'/)
02380 LERR=LERR+1
02390 IF(LERR)401,401,101
02400 800 RETURN
02410 450 TYPE 9450
02420 9450 FORMAT(' ENTER LINE NUMBER, THEN VARIABLE NUMBER,
02430 + THEN THE CORRECT VALUE.'/' IF Y IS TO BE
02440 + CORRECTED, THEN ENTER "0" FOR'/
02450 + ' VARIABLE NUMBER.'/)
02460 460 TYPE 9460
02470 9460 FORMAT('+ WHICH LINE, NUMBER, VALUE? ',$)
02480 9461 FORMAT(2I,F)
02490 ACCEPT 9461,LINE,NUMBER,VALUE
02500 IF(LINE)444,444,471
02510 471 IF(NUMBER)444,472,472
02520 472 IF(LINE-N)473,473,444
02530 473 IF(NUMBER-NP)474,474,444
02540 474 IF(NUMBER)475,475,476
02550 475 Y(LINE)=VALUE
02560 GO TO 444
02570 476 X(LINE,NUMBER)=VALUE
02580 GO TO 444
02590 END