	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
