SUBROUTINE POLY(Y,X,W,N,IDEG,IWRITE,B,SDB,S,DF,PRED,RES) C C PURPOSE--THIS SUBROUTINE COMPUTES A LEAST SQUARES C POLYNOMIAL FIT (OF DEGREE = IDEG) OF THE C RESPONSE VARIABLE DATA IN THE SINGLE PRECISION C VECTOR Y AS A FUNCTION OF THE INDEPENDENT C VARIABLE DATA IN THE SINGLE PRECISION C VECTOR X. C INPUT ARGUMENTS--Y = SINGLE PRECISION VECTOR OF C RESPONSE DATA (THAT IS, THE C DEPENDENT VARIABLE). C --X = SINGLE PRECISION VECTOR OF C THE INDEPENDENT VARIABLE. C --W = THE SINGLE PRECISION VECTOR C OF WEIGHTS FOR THE RESPONSE C VARIABLE. C --N = THE INTEGER VALUE OF THE SAMPLE SIZE. C --IDEG = THE INTEGER VALUE OF THE DESIRED C DEGREE OF THE POLYNOMIAL C TO BE FIT. C --IWRITE = THE INTEGER VALUE WHICH IF ZERO WILL C RESULT IN NO PRINTED OUTPUT, AND IF C NON-ZERO (E.G., 1) WILL RESULT IN C SOME LIMITED PRINTED OUTPUT C (COEFFICIENTS, STANDARD DEVIATIONS OF C COEFFICIENTS, RESIDUAL STANDARD DEVIATION). C OUTPUT ARGUMENTS--B = THE SINGLE PRECISION VECTOR OF C ESTIMATED REGRESSION COEFFICIENTS. C --SDB = THE SINGLE PRECISION VECTOR OF C ESTIMATED STANDARD DEVIATIONS OF THE C ESTIMATED REGRESSION COEFFICIENTS. C --S = THE ESTIMATED RESIDUAL STANDARD C DEVIATION. C --DF = THE DEGREES OF FREEDOM C ASSOCIATED WITH THE RESIDUAL C STANDARD DEVIATION = C NUMBER OF OBSERVATIONS MINUS C NUMBER OF PARAMETERS = C N - (IDEG + 1). C --PRED = THE SINGLE PRECISION VECTOR OF C PREDICTED VALUES FROM THE C LEAST SQUARES FIT. C --RES = THE SINGLE PRECISION VECTOR OF C RESIDUALS FROM THE LEAST SQUARES FIT. C SUBROUTINES NEEDED--DECOMP, INVXWX, DOT, FCDF. C WRITTEN BY--JAMES J. FILLIBEN C STATISTICAL ENGINEERING LABORATORY C NATIONAL BUREAU OF STANDARDS C WASHINGTON, D. C. 20234 C PHONE: 301-921-2315 C ORIGINAL VERSION--MARCH 1974 C UPDATED--OCTOBER 1974 C UPDATED--MARCH 1975 C UPDATED--MAY 1975 C UPDATED--JULY 1975 C UPDATED--SEPTEMBER 1975 C UPDATED--NOVEMBER 1975. C UPDATED--FEBRUARY 1976. C UPDATED--JUNE 1976. C UPDATED--OCTOBER 1976. C UPDATED--MAY 1977. C UPDATED--JUNE 1977. C C--------------------------------------------------------------------- C DIMENSION Y(1),X(1),W(1),B(1),SDB(1),PRED(1),RES(1) DIMENSION B2(50) DIMENSION F(3000),WRES(3000),G(50),H(50) DIMENSION Q(10000),R(2500),D(50),IPIVOT(50) COMMON /BLOCK2/ WS(15000) COMMON /BLOCK3/ DUM1(3000),DUM2(3000) EQUIVALENCE (Q(1),WS(1)) EQUIVALENCE (R(1),WS(10001)) EQUIVALENCE (D(1),WS(12501)) EQUIVALENCE (IPIVOT(1),WS(12551)) C IPR=6 AN=N K=IDEG+1 AK=K NK=N*K KK=K*K NMAX=3000 KMAX=50 NKMAX=10000 C C-----START POINT----------------------------------------------------- C C WRITE OUT THE TITLE C IF(IWRITE.EQ.0)GOTO150 WRITE(IPR,999) WRITE(IPR,999) WRITE(IPR,101) WRITE(IPR,102)N WRITE(IPR,103)IDEG C C PRE-SET THE OUTPUT VARIABLES AND VECTORS C TO A LARGE VALUE C IN CASE A PREMATURE EXIT OCCURS DUE TO NUMERICAL INSTABILITY. C 150 VALUE=(10.0**10)+1000.0 S=VALUE IF(K.LE.0)GOTO155 DO160I=1,K B(I)=VALUE SDB(I)=VALUE 160 CONTINUE 155 IF(N.LE.0)GOTO175 DO170I=1,N RES(I)=VALUE 170 CONTINUE C C CHECK THE INPUT ARGUMENTS N AND K C 175 IF(N.LE.0.OR.N.GT.NMAX)WRITE(IPR,105)N,NMAX IF(N.LE.0.OR.N.GT.NMAX)RETURN IF(K.LE.0.OR.K.GT.KMAX)WRITE(IPR,110)K,KMAX IF(K.LE.0.OR.K.GT.KMAX)RETURN IF(K.GT.N)WRITE(IPR,115)K,N IF(K.GT.N)RETURN C C INSPECT THE WEIGHT VECTOR W--IF ALL ELEMENTS ARE IDENTICAL, C THEN RESET ALL ELEMENTS TO 1.0. THIS AVOIDS THE C PROBLEM OF AN UNDEFINED EMPTY WEIGHT VECTOR W WHEN C IN FACT AN EQUAL WEIGHTING SCHEME IS DESIRED. C IWFLAG=0 WHOLD=W(1) DO600I=1,N IF(W(I).EQ.WHOLD)GOTO600 GOTO850 600 CONTINUE IWFLAG=1 850 IF(IWFLAG.EQ.0.AND.IWRITE.NE.0)WRITE(IPR,851) IF(IWFLAG.EQ.1.AND.IWRITE.NE.0)WRITE(IPR,852) C C COMPUTE THE ORIGINAL FORM FOR THE Q MATRIX C WHICH WILL BE IDENTICAL TO THE DATA MATRIX X C OF INDEPENDENT VARIABLES IF THE WEIGHTS C SPECIFIED ARE ALL EQUAL. NOTE THAT THE C DATA MATRIX X IS NEVER COMPUTED AS SUCH. C NOTE THAT THE DATA MATRIX X IS NOT TO BE C CONFUSED WITH THE SINGLE INDEPENDENT C VARIABLE VECTOR X. C THE Q MATRIX WILL BE CHANGED IN THE DECOMP SUBROUTINE. C DO100J=1,K IF(J.EQ.1)GOTO250 IF(J.EQ.2)GOTO350 JM1=J-1 DO200I=1,N IQARG1=(I-1)*K+J IQARG2=(I-1)*K+JM1 Q(IQARG1)=Q(IQARG2)*X(I) 200 CONTINUE GOTO100 250 DO300I=1,N IQARG=(I-1)*K+1 Q(IQARG)=1.0 300 CONTINUE GOTO100 350 DO400I=1,N IQARG=(I-1)*K+2 Q(IQARG)=X(I) 400 CONTINUE 100 CONTINUE IF(IWFLAG.EQ.1)GOTO1440 DO1300I=1,N DO1400J=1,K IQARG=(I-1)*K+J Q(IQARG)=Q(IQARG)*SQRT(W(I)) 1400 CONTINUE 1300 CONTINUE C C COMPUTE ETA AND TOL (FOR THE UNIVAC 1108, ETA = 2**-27) C WHICH WILL BE USED IN THE DECOMP SUBROUTINE C 1440 ETA=1.0 1450 ETA=0.5*ETA ETAP1=ETA+1.0 IF(ETAP1.GT.1.0)GOTO1450 TOL=ETA*AK NM5=N-5 NKM5=NK-5 CCCCC WRITE(IPR,1505)(Y(I),I=1,6) CCCCC WRITE(IPR,1505)(Y(I),I=NM5,N) CCCCC WRITE(IPR,1505)(X(I),I=1,6) CCCCC WRITE(IPR,1505)(X(I),I=NM5,N) CCCCC WRITE(IPR,1505)(W(I),I=1,6) CCCCC WRITE(IPR,1505)(W(I),I=NM5,N) CCCCC WRITE(IPR,1505)(Q(IQARG),IQARG=1,6) CCCCC WRITE(IPR,1505)(Q(IQARG),IQARG=NKM5,NK) CCCCC WRITE(IPR,1505)(Q(IQARG),IQARG=1,NK) C CALL DECOMP(N,K,ETA,TOL,IRANK,INSING) C CCCCC WRITE(IPR,1505)(Q(IQARG),IQARG=1,NK) CCCCC WRITE(IPR,1505)(R(IRARG),IRARG=1,KK) CCCCC WRITE(IPR,1505)(D(J),J=1,K) IF(INSING.EQ.1)GOTO1550 WRITE(IPR,1510)IRANK,K RETURN 1550 KP1=K+1 C C *************************************************************** C C THE PURPOSE OF THIS NEXT SEGMENT (BETWEEN THE STARRED LINES) C IS TO SOLVE FOR THE DESIRED REGRESSION COEFFICIENTS. C ITERATIVE REFINEMENT IS USED. C A SECOND OUTPUT FROM THIS SEGMENT IS THE RESIDUALS FROM THE C FINAL FIT. C A THIRD OUTPUT FROM THIS SEGMENT IS AN INDICATION (IN THE C VARIABLE ICONV) AS TO WHETHER THE ITERATIVE REFINEMENT C CONVERGED OR NOT. C X--USED IN THIS SEGMENT C Q--USED IN THIS SEGMENT C R--USED IN THIS SEGMENT C D--USED IN THIS SEGMENT C IPIVOT--USED IN THIS SEGMENT C C ETA2=ETA*ETA B2(KP1)=-1.0 DO 50010 I=1,N F(I)=Y(I) WRES(I)=0.0 RES(I)=0.0 50010 CONTINUE DO 50020 J=1,K B2(J)=0.0 G(J)=0.0 H(J)=0.0 50020 CONTINUE M=0 AMDB2=0.0 AMDR2=0.0 C C BEGIN THE M-TH ITERATION STEP IN THE ITERATIVE REFINEMENT C 50030 IF (M.LT.2) GO TO 50040 IF (((64.*AMDB2.LT.AMDB1).AND.(AMDB2.GT.ETA2*AMB)).OR.((64.*AMDR2. 1LT.AMDR1).AND.(AMDR2.GT.ETA2*AMR))) GO TO 50040 GO TO 50250 50040 AMDB1=AMDB2 AMDR1=AMDR2 AMDB2=0.0 AMDR2=0.0 IF (M.EQ.0) GO TO 50100 C C BEGIN FORMING NEW RESIDUALS C IF(IWFLAG.EQ.0)GOTO50044 DO50045I=1,N WRES(I)=WRES(I)+F(I) RES(I)=RES(I)+F(I) 50045 CONTINUE GOTO50055 50044 DO 50050 I=1,N WRES(I)=WRES(I)+F(I)*SQRT(W(I)) RES(I)=RES(I)+F(I)/SQRT(W(I)) 50050 CONTINUE 50055 DO 50070 IS=1,IRANK J=IPIVOT(IS) B2(J)=B2(J)+G(IS) DO 50060 L=1,N IF(J.GE.2)GOTO50065 DUM1(L)=1.0 GOTO50060 50065 DUM1(L)=X(L)**(J-1) 50060 CONTINUE C CALL DOT(DUM1,WRES,1,N,0.0,G(IS)) 50070 G(IS)=-G(IS) C DO 50090 I=1,N E=RES(I) DO 50080 L=1,K IF(L.GE.2)GOTO50085 DUM1(L)=1.0 GOTO50080 50085 DUM1(L)=X(I)**(L-1) 50080 CONTINUE DUM1(KP1)=Y(I) C CALL DOT(DUM1,B2,1,KP1,E,F(I)) IF(IWFLAG.EQ.0)F(I)=-F(I)*SQRT(W(I)) IF(IWFLAG.EQ.1)F(I)=-F(I) 50090 CONTINUE C C END FORMING NEW RESIDUALS C 50100 DO 50150 IS=1,IRANK J=IPIVOT(IS) IF (IS.NE.1) GO TO 50110 50110 ISM1=IS-1 DO 50120 L=1,ISM1 IRARG=(L-1)*K+IS 50120 DUM1(L)=R(IRARG) ANEGGI=-G(IS) C CALL DOT(DUM1,H,1,ISM1,ANEGGI,H(IS)) H(IS)=-H(IS) C E=-H(IS) DO 50130 L=1,N IQARG=(L-1)*K+J 50130 DUM1(L)=Q(IQARG) C CALL DOT(DUM1,F,1,N,E,E) E=E/D(IS) C G(IS)=E DO 50140 I=1,N IQARG=(I-1)*K+J 50140 F(I)=F(I)-E*Q(IQARG) 50150 CONTINUE C DO 50210 IS=1,IRANK JS=IRANK+1-IS JSP1=JS+1 DO 50200 L=JSP1,IRANK IRARG=(JS-1)*K+L 50200 DUM1(L)=R(IRARG) ANEGGJ=-G(JS) C CALL DOT(DUM1,G,JSP1,IRANK,ANEGGJ,G(JS)) 50210 G(JS)=-G(JS) C DO 50220 IS=1,IRANK 50220 AMDB2=AMDB2+G(IS)*G(IS) DO 50230 I=1,N 50230 AMDR2=AMDR2+F(I)*F(I) IF (M.NE.0) GO TO 50240 AMB=AMDB2 AMR=AMDR2 50240 CONTINUE CCCCC WRITE(IPR,50505)M CCCCC WRITE(IPR,50506)(B2(I),I=1,K) CCCCC DO5555I=1,N CCCCC WRITE(IPR,50506)Y(I),X(I),RES(I),WRES(I),F(I),G(I),H(I) C5555 CONTINUE 50505 FORMAT(I8) 50506 FORMAT(1H ,8F10.5) C C END THE M-TH ITERATION STEP IN THE ITERATIVE REFINEMENT C M=M+1 GO TO 50030 50250 IF ((AMDR2.GT.4.*ETA2*AMR).AND.(AMDB2.GT.4.*ETA2*AMB)) GO TO 50260 ICONV=1 GO TO 1700 50260 ICONV=0 C C *************************************************************** C 1700 IF(ICONV.EQ.1)GOTO1750 WRITE(IPR,1520) RETURN C C ADJUST THE R MATRIX C WHICH WILL BE USED IN THE INVERT (INVRR) SUBROUTINE. C 1750 DO1800J=1,K IRARG=(J-1)*K+J R(IRARG)=SQRT(D(J)) IF(J.EQ.K)GOTO1800 JP1=J+1 DO1900L=JP1,K IRARG1=(J-1)*K+L IRARG2=(J-1)*K+J R(IRARG1)=R(IRARG1)*R(IRARG2) 1900 CONTINUE 1800 CONTINUE CCCCC WRITE(IPR,1505)(R(IRARG),IRARG=1,KK) C CALL INVXWX(N,K) C CCCCC WRITE(IPR,1505)(R(IRARG),IRARG=1,KK) C C COMPUTE STATISTICAL CALCULATIONS AND THEN WRITE OUT COEFFICIENTS C AND STANDARD DEVIATIONS OF COEFFICIENTS ALONG WITH THE C RESIDUAL STANDARD DEVIATION. C DO3040I=1,K B(I)=B2(I) 3040 CONTINUE DO3070I=1,N IF(RES(I).EQ.0.0)GOTO3070 GOTO3050 3070 CONTINUE WRITE(IPR,3080) 3080 FORMAT(' ',10X,'NOTE THAT AN EXACT FIT HAS BEEN OBTAINED') 3050 SUM=0.0 IF(IWFLAG.EQ.0)GOTO3060 DO3055I=1,N SUM=SUM+RES(I)*RES(I) 3055 CONTINUE RESSS=SUM GOTO3105 3060 DO3100I=1,N SUM=SUM+RES(I)*RES(I)*W(I) 3100 CONTINUE 3105 IF(K.EQ.N.AND.IWRITE.NE.0)WRITE(IPR,3110)K IF(K.EQ.N)GOTO3250 3110 FORMAT(1H ,10X,72HNOTE THAT THE NUMBER OF COEFFICIENTS K 1 = THE SAMPLE SIZE N = ,I8) RESDF=AN-AK DF=RESDF IRESDF=AN-AK+0.5 RESMS=RESSS/RESDF S=SQRT(RESMS) 3250 CONTINUE CCCCC WRITE(IPR,3205) DO3300I=1,N PRED(I)=Y(I)-RES(I) CCCCC WRITE(IPR,3305)Y(I),PRED(I),RES(I) 3300 CONTINUE IF(K.EQ.N)GOTO3450 GOTO3480 3450 IF(IWRITE.NE.0)WRITE(IPR,3455) 3455 FORMAT(1H ,21H J B(J)) DO3460I=1,K IF(IWRITE.NE.0)WRITE(IPR,3405)J,B(J) 3460 CONTINUE RETURN 3480 IF(IWRITE.NE.0)WRITE(IPR,3310)S IF(IWRITE.NE.0)WRITE(IPR,3311)IRESDF IF(IWRITE.NE.0)WRITE(IPR,3312) IF(IWRITE.NE.0)WRITE(IPR,3315) DO3400J=1,K IRARG=(J-1)*K+J SDB(J)=S*SQRT(R(IRARG)) T=B(J)/SDB(J) IF(IWRITE.NE.0)WRITE(IPR,3405)J,B(J),SDB(J),T 3400 CONTINUE C C COMPUTE THE COVARIANCE AND CORRELATION MATRIX OF THE COEFFICIENTS C DO3500I=1,K DO3600J=1,K IRARG=(I-1)*K+J R(IRARG)=R(IRARG)*S*S 3600 CONTINUE 3500 CONTINUE DO2200I=1,K DO2300J=1,K IF(I.EQ.J)GOTO2300 IRARG1=(I-1)*K+J IRARG2=(I-1)*K+I IRARG3=(J-1)*K+J R(IRARG1)=R(IRARG1)/SQRT(R(IRARG2)*R(IRARG3)) 2300 CONTINUE 2200 CONTINUE DO2400J=1,K IRARG=(J-1)*K+J R(IRARG)=1.0 2400 CONTINUE CCCCC WRITE(IPR,2405) DO2500I=1,K IRAMIN=(I-1)*K+1 IRAMAX=I*K CCCCC WRITE(IPR,2505)(R(IRARG),IRARG=IRAMIN,IRAMAX) 2500 CONTINUE C C CHECK FOR REPLICATION. C IF SO, COMPUTE A POOLED STANDARD DEVIATION. C THEN COMPUTE A LACK OF FIT F TEST. C C DETERMINE THE NUMBER OF DISTINCT SUBSETS C NUMSET=0 DO4200I=1,N IF(NUMSET.EQ.0)GOTO4350 DO4300J=1,NUMSET IF(X(I).EQ.DUM1(J))GOTO4200 4300 CONTINUE 4350 NUMSET=NUMSET+1 DUM1(NUMSET)=X(I) 4200 CONTINUE IF(NUMSET.EQ.0)WRITE(IPR,4205) IF(NUMSET.EQ.0)RETURN C C COPY OUT EACH SUBSET INTO THE DUM2 VECTOR C AND ANALYZE IT THEREIN C IREPDF=0 REPSS=0.0 DO4600ISET=1,NUMSET NI=0 DO4700I=1,N IF(X(I).EQ.DUM1(ISET))NI=NI+1 IF(X(I).EQ.DUM1(ISET))DUM2(NI)=Y(I) 4700 CONTINUE ANI=NI SUM=0.0 DO5100I=1,NI SUM=SUM+DUM2(I) 5100 CONTINUE YMEAN=SUM/ANI SUM=0.0 DO5200I=1,NI SUM=SUM+(DUM2(I)-YMEAN)**2 5200 CONTINUE IREPDF=IREPDF+NI-1 REPSS=REPSS+SUM 4600 CONTINUE IF(IREPDF.LE.0)RETURN REPDF=IREPDF REPV=REPSS/REPDF REPSD=SQRT(REPV) IF(IWRITE.NE.0)WRITE(IPR,999) IF(IWRITE.NE.0)WRITE(IPR,4930)REPSD IF(IWRITE.NE.0)WRITE(IPR,4935)IREPDF IF(IWRITE.NE.0)WRITE(IPR,4933)NUMSET IFITDF=IRESDF-IREPDF IF(IFITDF.GE.1)GOTO5250 IF(IWRITE.NE.0)WRITE(IPR,4936) IF(IWRITE.NE.0)WRITE(IPR,4937) IF(IWRITE.NE.0)WRITE(IPR,4938) IF(IWRITE.NE.0)WRITE(IPR,4939) RETURN 5250 FITDF=IFITDF FITSS=RESSS-REPSS FITMS=FITSS/FITDF FSTAT=FITMS/RESMS CALL FCDF(FSTAT,IFITDF,IREPDF,CDF) CDF2=100.0*CDF IF(IWRITE.NE.0)WRITE(IPR,4940)FSTAT,CDF2 IF(IWRITE.NE.0)WRITE(IPR,4945)IFITDF,IREPDF C 999 FORMAT(1H ) 1505 FORMAT(1H ,6E15.8) 1510 FORMAT(1H ,51H*****ERROR--THE MATRIX IS SINGULAR--IT HAS IRANK = , 1I8,51H WHICH IS LESS THAN THE NUMBER OF COEFFICIENTS K = ,I8,6H ** 1***) 1520 FORMAT(1H ,51H*****ERROR--THE ITERATIONS ARE NOT CONVERGING *****) 1605 FORMAT(1H ,I8,5X,E15.8) 2005 FORMAT(1H ,8E15.8) 2405 FORMAT(1H ,34HCORRELATION MATRIX OF COEFFICIENTS) 2505 FORMAT(1H ,8E15.8) 3205 FORMAT(1H ,44H OBSERVED PREDICTED RESIDUALS) 3305 FORMAT(1H ,8E15.8) 3310 FORMAT(1H ,10X,30HRESIDUAL STANDARD DEVIATION = ,E15.8) 3311 FORMAT(1H ,10X,30HRESIDUAL DEGREES OF FREEDOM = ,I8) 3312 FORMAT(1H ,10X,13HCOEFFICIENTS:) 3315 FORMAT(1H ,58H J B(J) SD(B(J)) B(J)/SD 1(B(J))) 3405 FORMAT(1H ,3X,I8,8E15.8) 101 FORMAT(1H ,28HLEAST SQUARES POLYNOMIAL FIT) 102 FORMAT(1H ,10X,16HSAMPLE SIZE N = ,I8) 103 FORMAT(1H ,10X,9HDEGREE = ,I8) 105 FORMAT(1H ,33H*****ERROR--THE SAMPLE SIZE N (= ,I8,40H) IS NON-POS 1ITIVE OR LARGER THAN NMAX = ,I8,6H *****) 110 FORMAT(1H ,52H*****ERROR--THE DESIRED NUMBER OF COEFFICIENTS K (= 1,I8,40H) IS NON-POSITIVE OR LARGER THAN KMAX = ,I8,6H *****) 115 FORMAT(1H ,52H*****ERROR--THE DESIRED NUMBER OF COEFFICIENTS K (= 1,I8,38H) IS LARGER THAN THE SAMPLE SIZE N (= ,I8,7H) *****) 851 FORMAT(1H ,10X,20HUNEQUAL WEIGHTS CASE) 852 FORMAT(1H ,10X,18HEQUAL WEIGHTS CASE) 4205 FORMAT(1H ,38HERROR IN POLY SUBROUTINE--NUMSET = 0) 4930 FORMAT(1H ,44H REPLICATION STANDARD DEVIATION = ,D15.7) 4935 FORMAT(1H ,44H REPLICATION DEGREES OF FREEDOM = ,I8) 4933 FORMAT(1H ,44H NUMBER OF DISTINCT SUBSETS = ,I8) 4936 FORMAT(1H ,51H LACK OF FIT F TEST CANNOT BE DONE BECAUSE) 4937 FORMAT(1H ,44H HAVE ONLY 0 DEGREES OF FREEDOM IN , , 121HNUMERATOR OF F RATIO.) 4938 FORMAT(1H ,49H THIS HAPPENS WHEN NUMBER OF PARAMETERS , 16HFITTED ) 4939 FORMAT(1H ,45H IS IDENTICAL TO NUMBER OF DISTINCT , 18HSUBSETS.) 4940 FORMAT(1H ,32H LACK OF FIT F RATIO = ,F10.4,7H = THE , 1F8.4,14H% POINT OF THE) 4945 FORMAT(1H ,30H F DISTRIBUTION WITH ,I6,5H AND ,I6, 119H DEGREES OF FREEDOM) C RETURN END