Trailing-Edge
-
PDP-10 Archives
-
decus_20tap1_198111
-
decus/20-0025/mreg.for
There is 1 other file named mreg.for in the archive. Click here to see a list.
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