*BARTLT SUBROUTINE BARTLT (VAR,DF,N,CH2,CH2CDF,VARP,DFP,IFLAG) C C----------------------------------------------------------------------- C BARTLT WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: PERFORMING BARTLETT'S TEST FOR HOMOGENEITY OF VARIANCES ON C THREE OR MORE VARIANCES (THE F TEST SHOULD BE USED IN THE CASE C OF TWO VARIANCES). IF THE INPUT PARAMETERS ARE NOT VALID AN C ERROR FLAG IS SET AND NOTHING FURTHER IS COMPUTED, OTHERWISE C THE FOLLOWING ARE COMPUTED: C C 1) THE CHI-SQUARED STATISTIC (CH2), C 2) THE CUMULATIVE DISTRIBUTION FUNCTION OF THE CHI-SQUARED C DISTRIBUTION EVALUATED AT CH2 (CH2CDF), AND C 3) THE POOLED VARIANCE (VARP) AND ITS CORRESPONDING C DEGREES OF FREEDOM (DFP) C C THE VALUES IN 3) MAY BE USEFUL ONLY IF THE VARIANCES ARE C DETERMINED TO BE EQUAL. THE VALUE OF CH2CDF IS GOOD TO SIX C DECIMAL PLACES. C C SUBPROGRAMS CALLED: CDFGAM (GAMMA CUMULATIVE DISTRIBUTION FUNCTION) C C CURRENT VERSION COMPLETED FEBRUARY 3, 1987 C C REFERENCES: C C 1) SNEDECOR, GEORGE W. AND COCHRAN, WILLIAM G., 'STATISTICAL C METHODS', 6TH EDITION, IOWA STATE UNIVERSITY PRESS, PP. 296-298. C C 2) BROWNLEE, K.A., 'STATISTICAL THEORY AND METHODOLOGY IN SCIENCE C AND ENGINEERING', JOHN WILEY & SONS, 1960, PP. 225-227. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * VAR = VECTOR (LENGTH N) OF VARIANCES (REAL) C C * DF = VECTOR (LENGTH N) OF DEGREES OF FREEDOM CORRESPONDING C TO THE VARIANCES (REAL) C C * N = NUMBER OF VARIANCES [>2] (INTEGER) C C CH2 = THE CHI-SQUARED STATISTIC ASSOCIATED WITH BARTLETT'S TEST C (REAL) C C CH2CDF = THE CUMULATIVE DISTRIBUTION FUNCTION OF THE CHI-SQUARED C DISTRIBUTION WITH N-1 DEGREES OF FREEDOM EVALUATED AT CH2 C (REAL) C C VARP = THE POOLED VARIANCE DETERMINED FROM THE N VARIANCES (REAL) C C DFP = THE DEGREES OF FREEDOM ASSOCIATED WITH THE POOLED C VARIANCE (REAL) C C IFLAG = THE ERROR FLAG ON OUTPUT (INTEGER) INTERPRETATION: C 0 -> NO ERRORS DETECTED C 1,2 -> ERROR FLAGS FROM SUBROUTINE CDFGAM C 3 -> N<3 C 4 -> AT LEAST ONE DF(I) IS <= 0.0 C 5 -> AT LEAST ONE VARIANCE(I) IS < 0.0 C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION VAR(*),DF(*) C C--- CHECK VALIDITY OF INPUT PARAMETERS C IF (N.LT.3) THEN IFLAG = 3 RETURN ENDIF DO 10 I = 1, N IF (DF(I).LE.0.0) THEN IFLAG = 4 RETURN ENDIF IF (VAR(I).LT.0.0) THEN IFLAG = 5 RETURN ENDIF 10 CONTINUE C C--- COMPUTE NEEDED SUMMATIONS C A = 0.0 VARP = 0.0 C = 0.0 CH2 = 0.0 DO 20 I = 1, N A = A+DF(I) VARP = VARP+DF(I)*VAR(I) C = C+1.0/DF(I) CH2 = CH2+DF(I)*ALOG(VAR(I)) 20 CONTINUE C C--- COMPUTE THE POOLED VARIANCE AND ITS DEGREES OF FREEDOM C VARP = VARP/A DFP = A C C--- COMPUTE THE CHI-SQUARED STATISTIC C CH2 = A*ALOG(VARP)-CH2 A = 1.0+(C-1.0/A)/(3.0*REAL(N-1)) CH2 = CH2/A C C--- COMPUTE THE C.D.F. AT CH2 C X = 0.5*CH2 ALPHA = 0.5*REAL(N-1) EPS = 0.000001 CALL CDFGAM (X,ALPHA,EPS,IFLAG,CH2CDF) RETURN END *CDFBET SUBROUTINE CDFBET (X,P,Q,EPS,IFLAG,CDFX) C C----------------------------------------------------------------------- C CDFBET WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING THE CUMULATIVE DISTRIBUTION FUNCTION OF THE BETA C DISTRIBUTION (ALSO KNOWN AS THE INCOMPLETE BETA RATIO) TO A C SPECIFIED ACCURACY (TRUNCATION ERROR IN THE INFINITE SERIES). C THE ALGORITHM, DESCRIBED IN REFERENCE 2, IS A MODIFICATION OF C THE ALGORITHM OF REFERENCE 1. THREE FEATURES HAVE BEEN ADDED: C C 1) A PRECISE METHOD OF MEETING THE TRUNCATION ACCURACY, C 2) A CONSTANT W USED IN DETERMINING FOR WHICH X VALUES THE C RELATION I(X,P,Q) = 1 - I(1-X,Q,P) IS TO BE USED, AND C 3) A CONSTANT UFLO >= THE UNDERFLOW LIMIT ON THE COMPUTER. C C SUBPROGRAMS CALLED: DGAMLN (LOG OF GAMMA FUNCTION) C C CURRENT VERSION COMPLETED OCTOBER 24, 1986 C C REFERENCES: C C 1) MAJUMDER, K.L. AND BHATTACHARJEE, G.P., 'THE INCOMPLETE BETA C INTEGRAL', ALGORITHM AS 63, APPLIED STATISTICS, VOL. 22, NO. 3, C 1973, PP. 409-411. C C 2) REEVE, CHARLES P., 'AN ALGORITHM FOR COMPUTING THE BETA C.D.F. C TO A SPECIFIED ACCURACY', STATISTICAL ENGINEERING DIVISION C NOTE 86-3, OCTOBER 1986. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * X = VALUE AT WHICH THE C.D.F. IS TO BE COMPUTED (REAL) C C * P = FIRST PARAMETER OF THE BETA FUNCTION (>0) (REAL) C C * Q = SECOND PARAMETER OF THE BETA FUNCTION (>0) (REAL) C C * EPS = THE DESIRED ABSOLUTE ACCURACY OF THE C.D.F. (>0) (REAL) C C IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION: C 0 -> NO ERRORS DETECTED C 1 -> EITHER P OR Q OR EPS IS <= UFLO C 2 -> NUMBER OF TERMS EVALUATED IN THE INFINITE SERIES C EXCEEDS JMAX C C CDFX = THE C.D.F. EVALUATED AT X (REAL) C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C LOGICAL LL DOUBLE PRECISION DP,DQ,DGAMLN CCCCC DATA JMAX,W,UFLO / 5000,20.0,1E-100 / DATA JMAX,W,UFLO / 5000,20.0,1E-30 / CDFX = 0.0 C C--- CHECK FOR VALIDITY OF ARGUMENTS P, Q, AND EPS C IF (P.LE.UFLO.OR.Q.LE.UFLO.OR.EPS.LE.UFLO) THEN IFLAG = 1 RETURN ENDIF IFLAG = 0 C C--- CHECK FOR SPECIAL CASES OF X C IF (X.LE.0.0) RETURN IF (X.GE.1.0) THEN CDFX = 1.0 ELSE C C--- SWITCH ARGUMENTS IF NECESSARY C LL = P+W.GE.(P+Q+2.0*W)*X IF (LL) THEN XY = X YX = 1.0-XY PQ = P QP = Q ELSE YX = X XY = 1.0-YX QP = P PQ = Q ENDIF C C--- EVALUATE THE BETA P.D.F. AND CHECK FOR UNDERFLOW C DP = DBLE(PQ-1.0)*DLOG(DBLE(XY))-DGAMLN(PQ) DQ = DBLE(QP-1.0)*DLOG(DBLE(YX))-DGAMLN(QP) PDFL = SNGL(DGAMLN(PQ+QP)+DP+DQ) IF (PDFL.LT.ALOG(UFLO)) THEN ELSE U = EXP(PDFL)*XY/PQ R = XY/YX 10 IF (QP.LE.1.0) GO TO 20 C C--- INCREMENT PQ AND DECREMENT QP C IF (U.LE.EPS*(1.0-(PQ+QP)*XY/(PQ+1.0))) GO TO 40 CDFX = CDFX+U PQ = PQ+1.0 QP = QP-1.0 U = QP*R*U/PQ GO TO 10 20 V = YX*U YXEPS = YX*EPS C C--- INCREMENT PQ C DO 30 J = 0, JMAX IF (V.LE.YXEPS) GO TO 40 CDFX = CDFX+V PQ = PQ+1.0 V = (PQ+QP-1.0)*XY*V/PQ 30 CONTINUE IFLAG = 2 ENDIF 40 IF (.NOT.LL) CDFX = 1.0-CDFX ENDIF RETURN END *CDFDNF SUBROUTINE CDFDNF (X,DF1,DF2,ALAMB1,ALAMB2,EPS,IFLAG,CDFX) C C----------------------------------------------------------------------- C CDFDNF WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING THE CUMULATIVE DISTRIBUTION FUNCTION OF THE DOUBLY C NONCENTRAL F DISTRIBUTION TO A SPECIFIED ACCURACY (TRUNCATION C ERROR IN THE INFINITE SERIES REPRESENTATION GIVEN BY EQUATION C 2.2 IN REFERENCE 1 BELOW). THE BETA C.D.F. ROUTINE IS CALLED C AT MOST TWO TIMES. FURTHER VALUES OF THE BETA C.D.F. ARE C OBTAINED FROM RECURRENCE RELATIONS GIVEN IN REFERENCE 2. C REFERENCE 3 GIVES A DETAILED DESCRIPTION OF THE ALGORITHM C HEREIN. C C THIS PROGRAM MAY ALSO BE EFFICIENTLY USED TO COMPUTE THE C CUMULATIVE DISTRIBUTION FUNCTIONS OF THE SINGLY NONCENTRAL C AND CENTRAL F DISTRIBUTIONS BY SETTING THE APPROPRIATE C NONCENTRALITY PARAMETERS EQUAL TO ZERO. C C CHECKS ARE MADE TO ASSURE THAT ALL PASSED PARAMETERS ARE C WITHIN VALID RANGES AS GIVEN BELOW. NO UPPER LIMIT IS SET C FOR THE NONCENTRALITY PARAMETERS, BUT VALUES UP TO ABOUT C 10,000 CAN BE HANDLED WITH THE CURRENT DIMENSION LIMITS. THE C COMPUTED VALUE CDFX IS VALID ONLY IF IFLAG=0 ON RETURN. C C NOTE: IN EQUATION 2.2 OF REFERENCE 1 THE AUTHOR HAS MISTAKENLY C REVERSED THE ARGUMENTS OF THE INCOMPLETE BETA FUNCTION. C THEY SHOULD READ [(M/2)+R,(N/2+S)] WHERE M AND N ARE THE C DEGREES OF FREEDOM ASSOCIATED WITH THE NUMERATOR AND C DENOMINATOR RESPECTIVELY OF THE F STATISTIC. TO FURTHER C CONFUSE THE ISSUE, THE AUTHOR HAS REVERSED THE USAGE OF C M AND N IN SECTION 1 OF THE PAPER. C C NOTE: IN SUBROUTINE EDGEF THE DOUBLE PRECISION CONSTANT DEUFLO IS C THE EXPONENTIAL UNDERFLOW LIMIT WHOSE CURRENT VALUE IS SET C AT -69D0. ON A COMPUTER WHERE DEXP(-69D0) CAUSES UNDERFLOW C THIS LIMIT SHOULD BE CHANGED. C C SUBPROGRAMS CALLED: CDFBET (BETA C.D.F.) C DGAMLN (DOUBLE PRECISION LOG OF GAMMA FUNCTION) C POISSF, EDGEF (ATTACHED) C C CURRENT VERSION COMPLETED SEPTEMBER 29, 1988 C C REFERENCES: C C 1. BULGREN, W.G., 'ON REPRESENTATIONS OF THE DOUBLY NONCENTRAL F C DISTRIBUTION', JOURNAL OF THE AMERICAN STATISTICAL ASSOCIATION, C MARCH 1971, VOLUME 66, NO. 333, PP. 184-186. C C 2. ABRAMOWITZ, MILTON, AND STEGUN, IRENE A., 'HANDBOOK OF C MATHEMATICAL FUNCTIONS', NATIONAL BUREAU OF STANDARDS APPLIED C MATHEMATICS SERIES 55, NOVEMBER 1970, P. 944. C C 3. REEVE, CHARLES P., 'AN ALGORITHM FOR COMPUTING THE DOUBLY C NONCENTRAL F C.D.F. TO A SPECIFIED ACCURACY', STATISTICAL C ENGINEERING DIVISION NOTE 86-4, NOVEMBER 1986. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * X = VALUE (>=0) AT WHICH THE C.D.F. IS TO BE COMPUTED (REAL) C C * DF1 = DEGREES OF FREEDOM (>0) IN THE NUMERATOR (REAL) C C * DF2 = DEGREES OF FREEDOM (>0) IN THE DENOMINATOR (REAL) C C * ALAMB1 = THE NONCENTRALITY PARAMETER (>=0) FOR THE NUMERATOR C (REAL) [EQUAL TO ZERO FOR THE CENTRAL F DISTRIBUTION] C C * ALAMB2 = THE NONCENTRALITY PARAMETER (>=0) FOR THE DENOMINATOR C (REAL) [EQUAL TO ZERO FOR THE SINGLY NONCENTRAL F AND C CENTRAL F DISTRIBUTIONS] C C * EPS = THE DESIRED ABSOLUTE ACCURACY OF THE C.D.F. (REAL) C [1 >= EPS >= 10**(-10)] C C IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION: C 0 -> NO ERRORS DETECTED C 1,2 -> ERROR FLAGS FROM SUBROUTINE CDFBET C 3 -> EITHER ALAMB1 OR ALAMB2 IS < 0 C 4 -> EITHER DF1 OR DF2 IS <= 0 C 5 -> EPS IS OUTSIDE THE RANGE [10**(-10),1] C 6 -> VECTOR DIMENSIONS ARE TOO SMALL - INCREASE NX C C CDFX = THE DOUBLY NONCENTRAL F C.D.F. EVALUATED AT X (REAL) C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C PARAMETER (NX=1000) DIMENSION BFI(NX),BFJ(NX),POI(NX),POJ(NX) CDFX = 0.0 C C--- CHECK VALIDITY OF ARGUMENTS C IF (ALAMB1.LT.0.0.OR.ALAMB2.LT.0.0) THEN IFLAG = 3 RETURN ENDIF IF (DF1.LE.0.0.OR.DF2.LE.0.0) THEN IFLAG = 4 RETURN ENDIF IF (EPS.GT.1.0.OR.EPS.LT.1.0E-10) THEN IFLAG = 5 RETURN ENDIF IFLAG = 0 C C--- SET ERROR CRITERION FOR THE BETA C.D.F. (PECULIAR TO CDFBET) C EPS3 = 0.001*EPS C FA = 0.5*ALAMB1 GA = 0.5*ALAMB2 FB = 0.5*DF1 GB = 0.5*DF2 YY = DF2/(DF2+DF1*X) IF (YY.GE.1.0) RETURN XX = 1.0-YY IF (XX.GE.1.0) THEN CDFX = 1.0 RETURN ENDIF C C--- COMPUTE POISSON PROBABILITIES IN VECTORS POI AND POJ C CALL POISSF (FA,EPS,IMIN,NI,POI,NX,IFLAG) IF (IFLAG.NE.0) RETURN FC = FB+REAL(IMIN) CALL POISSF (GA,EPS,JMIN,NJ,POJ,NX,IFLAG) IF (IFLAG.NE.0) RETURN GC = GB+REAL(JMIN) C C--- COMPUTE BETA C.D.F. BY RECURRENCE WHEN I=IMIN AND J=JMIN TO JMAX C CALL EDGEF (NJ,GC,FC,YY,XX,BFJ,CDFX,POJ,POI,EPS3,IFLAG,1) IF (NI.LE.1.OR.IFLAG.NE.0) RETURN C C--- COMPUTE BETA C.D.F. BY RECURRENCE WHEN J=JMIN AND I=IMIN TO IMAX C BFI(1) = BFJ(1) CALL EDGEF (NI,FC,GC,XX,YY,BFI,CDFX,POI,POJ,EPS3,IFLAG,2) IF (NJ.LE.1.OR.IFLAG.NE.0) RETURN C C--- COMPUTE BETA C.D.F. BY RECURRENCE WHEN I>IMIN AND J>JMIN C DO 20 I = 2, NI BFJ(1) = BFI(I) DO 10 J = 2, NJ BFJ(J) = XX*BFJ(J)+YY*BFJ(J-1) CDFX = CDFX+POI(I)*POJ(J)*BFJ(J) 10 CONTINUE 20 CONTINUE RETURN END C SUBROUTINE POISSF (ALAMB,EPS,L,NSPAN,V,NV,IFLAG) C C--- COMPUTE THE POISSON(ALAMB) PROBABILITIES OVER THE RANGE [L,K] C--- WHERE THE TOTAL TAIL PROBABILITY IS LESS THAN EPS/2, SUM THE C--- PROBABILITIES IN DOUBLE PRECISION, AND SHIFT THEM TO THE C--- BEGINNING OF VECTOR V. C DIMENSION V(*) DOUBLE PRECISION DAL,DK,DLIMIT,DSUM,DGAMLN DLIMIT = 1.0D0-0.5D0*DBLE(EPS) K = INT(ALAMB) L = K+1 IF (ALAMB.EQ.0.0) THEN PL = 1.0 ELSE DAL = DBLE(ALAMB) DK = DBLE(K) PL = SNGL(DEXP(DK*DLOG(DAL)-DAL-DGAMLN(REAL(K+1)))) ENDIF PK = ALAMB*PL/REAL(L) NK = NV/2 NL = NK+1 DSUM = 0.0 10 IF (PL.LT.PK) THEN NK = NK+1 IF (NK.GT.NV) THEN IFLAG = 6 RETURN ENDIF V(NK) = PK DSUM = DSUM+DBLE(PK) K = K+1 IF (DSUM.GE.DLIMIT) GO TO 20 PK = ALAMB*PK/REAL(K+1) ELSE NL = NL-1 V(NL) = PL DSUM = DSUM+DBLE(PL) L = L-1 IF (DSUM.GE.DLIMIT) GO TO 20 PL = REAL(L)*PL/ALAMB ENDIF GO TO 10 20 INC = NL-1 DO 30 I = NL, NK V(I-INC) = V(I) 30 CONTINUE NSPAN = NK-INC RETURN END C SUBROUTINE EDGEF (NK,FC,GC,XX,YY,BFK,CDFX,POI,POJ,EPS3,IFLAG,L) C C--- COMPUTE THE BETA C.D.F.'S BY A RECURRENCE RELATION ALONG THE EDGES C--- I = IMIN AND J = JMIN OF A GRID. THE CORRESPONDING COMPONENTS OF C--- THE F" C.D.F. ARE INCLUDED IN THE SUMMATION. TERMS WHICH MIGHT C--- CAUSE UNDERFLOW ARE SET TO ZERO. C DIMENSION BFK(*),POI(*),POJ(*) DOUBLE PRECISION DARG,DEUFLO,DGAMLN DATA DEUFLO / -69.0D0 / FD = FC-1.0 K = MAX0(L,MIN0(NK,INT((GC-1.0)*XX/YY-FD))) FK = FD+REAL(K) CALL CDFBET (XX,FK,GC,EPS3,IFLAG,BFK(K)) IF (IFLAG.NE.0) RETURN IF (L.EQ.1) BFK(K) = 1.0-BFK(K) IF (NK.EQ.1) GO TO 40 DARG = DBLE(FK)*DLOG(DBLE(XX))+DBLE(GC)*DLOG(DBLE(YY))- * DLOG(DBLE(FK))+DGAMLN(FK+GC)-DGAMLN(FK)-DGAMLN(GC) IF (DARG.LT.DEUFLO) THEN DK = 0.0 ELSE DK = SNGL(DEXP(DARG))*(-1.0)**L ENDIF IF (K.GE.NK) GO TO 20 BFK(K+1) = BFK(K)-DK DI = DK KFLAG = 1 DO 10 I = K+1, NK-1 IF (KFLAG.EQ.1) THEN DI = DI*(FD+GC+REAL(I-1))*XX/(FD+REAL(I)) IF (DK+DI.EQ.DK) THEN KFLAG = 0 DI = 0.0 ENDIF ENDIF BFK(I+1) = BFK(I)-DI 10 CONTINUE 20 DI = DK KFLAG = 1 DO 30 I = K-1, L, -1 IF (KFLAG.EQ.1) THEN DI = DI*(FC+REAL(I))/((FD+GC+REAL(I))*XX) IF (DK+DI.EQ.DK) THEN KFLAG = 0 DI = 0.0 ENDIF ENDIF BFK(I) = BFK(I+1)+DI 30 CONTINUE 40 DO 50 I = L, NK CDFX = CDFX+POI(I)*POJ(1)*BFK(I) 50 CONTINUE RETURN END *CDFDNT SUBROUTINE CDFDNT (X,DF,DELTA,ALAMB,EPS,IFLAG,CDFX) C C----------------------------------------------------------------------- C CDFDNT WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING THE CUMULATIVE DISTRIBUTION FUNCTION OF THE DOUBLY C NONCENTRAL T DISTRIBUTION TO A SPECIFIED ACCURACY (TRUNCATION C ERROR IN THE INFINITE SERIES REPRESENTATION GIVEN BY EQUATION C 4 IN REFERENCE 1 BELOW). WHEN X<0 THE C.D.F. IS COMPUTED C FROM CDF(X,DF,DELTA,ALAMB) = 1 - CDF(-X,DF,-DELTA,ALAMB). C THE BETA C.D.F. ROUTINE IS CALLED AT MOST FOUR TIMES. FURTHER C VALUES OF THE BETA C.D.F. ARE OBTAINED FROM RECURRENCE C RELATIONS GIVEN IN REFERENCE 2. REFERENCE 3 GIVES A DETAILED C DESCRIPTION OF THE ALGORITHM HEREIN. C C THIS PROGRAM MAY ALSO BE EFFICIENTLY USED TO COMPUTE THE C CUMULATIVE DISTRIBUTION FUNCTIONS OF THE SINGLY NONCENTRAL C AND CENTRAL T DISTRIBUTIONS BY SETTING THE APPROPRIATE C NONCENTRALITY PARAMETERS EQUAL TO ZERO. C C CHECKS ARE MADE TO ASSURE THAT ALL PASSED PARAMETERS ARE C WITHIN VALID RANGES AS GIVEN BELOW. NO UPPER LIMIT IS SET C FOR THE NONCENTRALITY PARAMETERS, BUT VALUES UP TO ABOUT 100 C FOR DELTA AND 10,000 FOR LAMBDA CAN BE HANDLED WITH THE C CURRENT DIMENSION LIMITS. THE COMPUTED VALUE CDFX IS VALID C ONLY IF IFLAG=0 ON RETURN. C C NOTE: IN SUBROUTINE EDGET THE DOUBLE PRECISION CONSTANT DEUFLO IS C THE EXPONENTIAL UNDERFLOW LIMIT WHOSE CURRENT VALUE IS SET C AT -69D0. ON A COMPUTER WHERE DEXP(-69D0) CAUSES UNDERFLOW C THIS LIMIT SHOULD BE CHANGED. C C SUBPROGRAMS CALLED: CDFBET (BETA C.D.F.) C DGAMLN (DOUBLE PRECISION LOG OF GAMMA FUNCTION) C POISST, EDGET, GRID (ATTACHED) C C CURRENT VERSION COMPLETED SEPTEMBER 29, 1988 C C REFERENCES: C C 1. KRISHNAN, MARAKATHA, 'SERIES REPRESENTATIONS OF THE DOUBLY C NONCENTRAL T DISTRIBUTION', JOURNAL OF THE AMERICAN STATISTICAL C ASSOCIATION, SEPTEMBER 1968, VOLUME 63, NO. 323, PP. 1004-1012. C C 2. ABRAMOWITZ, MILTON, AND STEGUN, IRENE A., 'HANDBOOK OF C MATHEMATICAL FUNCTIONS', NATIONAL BUREAU OF STANDARDS APPLIED C MATHEMATICS SERIES 55, NOVEMBER 1970, P. 944. C C 3. REEVE, CHARLES P., 'AN ALGORITHM FOR COMPUTING THE DOUBLY C NONCENTRAL T C.D.F. TO A SPECIFIED ACCURACY', STATISTICAL C ENGINEERING DIVISION NOTE 86-5, DECEMBER 1986. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * X = VALUE AT WHICH THE C.D.F. IS TO BE COMPUTED (REAL) C C * DF = DEGREES OF FREEDOM (>0) IN THE DENOMINATOR (REAL) C C * DELTA = THE NONCENTRALITY PARAMETER FOR THE NUMERATOR (REAL) C [EQUAL TO ZERO FOR THE CENTRAL T DISTRIBUTION] C C * ALAMB = THE NONCENTRALITY PARAMETER (>=0) FOR THE DENOMINATOR C (REAL) [EQUAL TO ZERO FOR THE SINGLY NONCENTRAL T AND C CENTRAL T DISTRIBUTIONS] C C * EPS = THE DESIRED ABSOLUTE ACCURACY OF THE C.D.F. (REAL) C [1 >= EPS >= 10**(-10)] C C IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION: C 0 -> NO ERRORS DETECTED C 1,2 -> ERROR FLAGS FROM SUBROUTINE CDFBET C 3 -> ALAMB IS < 0 C 4 -> DF IS <= 0 C 5 -> EPS IS OUTSIDE THE RANGE [10**(-10),1] C 6 -> VECTOR DIMENSIONS ARE TOO SMALL - INCREASE NX C C CDFX = THE DOUBLY NONCENTRAL T C.D.F. EVALUATED AT X (REAL) C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C PARAMETER (NX=1000) DIMENSION BFI(NX),BFJ(NX),POI(NX),POJ(NX) DOUBLE PRECISION DARG,DFA LOGICAL LL CDFX = 0.0 C C--- CHECK VALIDITY OF ARGUMENTS C IF (ALAMB.LT.0.0) THEN IFLAG = 3 RETURN ENDIF IF (DF.LE.0.0) THEN IFLAG = 4 RETURN ENDIF IF (EPS.GT.1.0.OR.EPS.LT.1.0E-10) THEN IFLAG = 5 RETURN ENDIF IFLAG = 0 C C--- SET ERROR CRITERION FOR THE BETA C.D.F. (PECULIAR TO CDFBET) C EPS3 = 0.001*EPS C DELSQ = DELTA**2 FA = 0.5*DELSQ print *,'delta,delsq,fa=',delta,delsq,fa GA = 0.5*ALAMB GB = 0.5*DF YY = DF/(DF+X*X) XX = 1.0-YY C C--- IF X<0 SET LL=.TRUE., REVERSE SIGN OF DELTA, AND USE THE C--- IDENTITY DESCRIBED UP FRONT FOR COMPUTING THE C.D.F. C LL = X.LT.0.0 IF (XX.GE.1.0) THEN CDFX = 1.0 GO TO 50 ENDIF SDELTA = DELTA IF (LL) SDELTA = -DELTA C C--- COMPUTE POISSON PROBABILITIES IN VECTOR POI C print *,'before poisst, fa=',fa CALL POISST (FA,EPS,IMIN,NI,POI,NX,IFLAG) print *,'after poisst, fa=',fa IF (IFLAG.NE.0) RETURN IF (YY.GE.1.0) GO TO 10 FC = 0.5+REAL(IMIN) C C--- COMPUTE POISSON PROBABILITIES IN VECTOR POJ C CALL POISST (GA,EPS,JMIN,NJ,POJ,NX,IFLAG) IF (IFLAG.NE.0) RETURN GC = GB+REAL(JMIN) C C--- SUM THE TERMS CORRESPONDING TO 'EVEN' VALUES OF INDEX I C CALL GRID (NI,NJ,FC,GC,BFI,BFJ,POI,POJ,XX,YY,EPS3,CDFX,IFLAG) IF (IFLAG.NE.0) RETURN 10 IF (DELTA.EQ.0.0) THEN NI = 0 SUM = 0.0 IF (YY.GE.1.0) GO TO 40 ELSE C C--- COMPUTE 'POISSON-LIKE' PROBABILITIES IN VECTOR POI C print *,'before k, fa=',fa K = INT(FA) IF (IMIN.GT.0) THEN IMIN = IMIN-1 NI = NI+1 ENDIF DFA = DBLE(FA) DARG = (DBLE(K)+0.5D0)*DLOG(DFA)-DFA-DGAMLN(REAL(K)+1.5) print *,'k,dfa,darg=',k,dfa,darg L = K-IMIN+1 POI(L) = SIGN(SNGL(DEXP(DARG)),SDELTA) SUM = POI(L) DO 20 I = K-1, IMIN, -1 L = L-1 POI(L) = POI(L+1)*(REAL(I)+1.5)/FA SUM = SUM+POI(L) 20 CONTINUE L = K-IMIN+1 DO 30 I = K+1, IMIN+NI-1 L = L+1 POI(L) = POI(L-1)*FA/(REAL(I)+0.5) SUM = SUM+POI(L) 30 CONTINUE IF (YY.GE.1.0) GO TO 40 FC = 1.0+REAL(IMIN) C C--- SUM THE TERMS CORRESPONDING TO 'ODD' VALUES OF INDEX I C CALL GRID (NI,NJ,FC,GC,BFI,BFJ,POI,POJ,XX,YY,EPS3,CDFX,IFLAG) print *,'after call grid, cdfx,iflag=',cdfx,iflag IF (IFLAG.NE.0) RETURN ENDIF C C--- COMPUTE THE NORMAL C.D.F. AT -SDELTA C 40 PHI = 0.5*(1.0-SUM) C C--- COMPUTE THE DOUBLY NONCENTRAL T C.D.F. AT X, USING AN IDENTITY C--- IF X<0 C CDFX = 0.5*CDFX+PHI 50 IF (LL) CDFX = 1.0-CDFX RETURN END C SUBROUTINE POISST (ALAMB,EPS,L,NSPAN,V,NV,IFLAG) C C--- COMPUTE THE POISSON(ALAMB) PROBABILITIES OVER THE RANGE [L,K] C--- WHERE THE TOTAL TAIL PROBABILITY IS LESS THAN EPS/3, SUM THE C--- PROBABILITIES IN DOUBLE PRECISION, AND SHIFT THEM TO THE C--- BEGINNING OF VECTOR V. C DIMENSION V(*) DOUBLE PRECISION DAL,DK,DLIMIT,DSUM,DGAMLN DLIMIT = 1.0D0-2.0D0*DBLE(EPS)/3.0D0 K = INT(ALAMB) L = K+1 IF (ALAMB.EQ.0.0) THEN PL = 1.0 ELSE DAL = DBLE(ALAMB) DK = DBLE(K) PL = SNGL(DEXP(DK*DLOG(DAL)-DAL-DGAMLN(REAL(K+1)))) ENDIF PK = ALAMB*PL/REAL(L) NK = NV/2 NL = NK+1 DSUM = 0.0 10 IF (PL.LT.PK) THEN NK = NK+1 IF (NK.GT.NV) THEN IFLAG = 6 RETURN ENDIF V(NK) = PK DSUM = DSUM+DBLE(PK) K = K+1 IF (DSUM.GE.DLIMIT) GO TO 20 PK = ALAMB*PK/REAL(K+1) ELSE NL = NL-1 V(NL) = PL DSUM = DSUM+DBLE(PL) L = L-1 IF (DSUM.GE.DLIMIT) GO TO 20 PL = REAL(L)*PL/ALAMB ENDIF GO TO 10 20 INC = NL-1 DO 30 I = NL, NK V(I-INC) = V(I) 30 CONTINUE NSPAN = NK-INC RETURN END C SUBROUTINE EDGET (NK,FC,GC,XX,YY,BFK,CDFX,POI,POJ,EPS3,IFLAG,L) C C--- COMPUTE THE BETA C.D.F.'S BY A RECURRENCE RELATION ALONG THE EDGES C--- I = IMIN AND J = JMIN OF A GRID. THE CORRESPONDING COMPONENTS OF C--- THE T" C.D.F. ARE INCLUDED IN THE SUMMATION. TERMS WHICH MIGHT C--- CAUSE UNDERFLOW ARE SET TO ZERO. C DIMENSION BFK(*),POI(*),POJ(*) DOUBLE PRECISION DARG,DEUFLO,DGAMLN DATA DEUFLO / -69.0D0 / FD = FC-1.0 K = MAX0(L,MIN0(NK,INT((GC-1.0)*XX/YY-FD))) FK = FD+REAL(K) CALL CDFBET (XX,FK,GC,EPS3,IFLAG,BFK(K)) IF (IFLAG.NE.0) RETURN IF (L.EQ.1) BFK(K) = 1.0-BFK(K) IF (NK.EQ.1) GO TO 40 DARG = DBLE(FK)*DLOG(DBLE(XX))+DBLE(GC)*DLOG(DBLE(YY))- * DLOG(DBLE(FK))+DGAMLN(FK+GC)-DGAMLN(FK)-DGAMLN(GC) IF (DARG.LT.DEUFLO) THEN DK = 0.0 ELSE DK = SNGL(DEXP(DARG))*(-1.0)**L ENDIF IF (K.GE.NK) GO TO 20 BFK(K+1) = BFK(K)-DK DI = DK KFLAG = 1 DO 10 I = K+1, NK-1 IF (KFLAG.EQ.1) THEN DI = DI*(FD+GC+REAL(I-1))*XX/(FD+REAL(I)) IF (DK+DI.EQ.DK) THEN KFLAG = 0 DI = 0.0 ENDIF ENDIF BFK(I+1) = BFK(I)-DI 10 CONTINUE 20 DI = DK KFLAG = 1 DO 30 I = K-1, L, -1 IF (KFLAG.EQ.1) THEN DI = DI*(FC+REAL(I))/((FD+GC+REAL(I))*XX) IF (DK+DI.EQ.DK) THEN KFLAG = 0 DI = 0.0 ENDIF ENDIF BFK(I) = BFK(I+1)+DI 30 CONTINUE 40 DO 50 I = L, NK CDFX = CDFX+POI(I)*POJ(1)*BFK(I) 50 CONTINUE RETURN END C SUBROUTINE GRID (NI,NJ,FC,GC,BFI,BFJ,POI,POJ,XX,YY,EPS3,CDFX,IFLAG * ) C C--- COMPUTE DOUBLE SUMMATION OF COMPONENTS OF THE T" C.D.F. OVER THE C--- GRID I=IMIN TO IMAX AND J=JMIN TO JMAX C DIMENSION BFI(*),BFJ(*),POI(*),POJ(*) C C--- COMPUTE BETA C.D.F. BY RECURRENCE WHEN I=IMIN, J=JMIN TO JMAX C CALL EDGET (NJ,GC,FC,YY,XX,BFJ,CDFX,POJ,POI,EPS3,IFLAG,1) IF (NI.LE.1.OR.IFLAG.NE.0) RETURN C C--- COMPUTE BETA C.D.F. BY RECURRENCE WHEN J=JMIN, I=IMIN TO IMAX C BFI(1) = BFJ(1) CALL EDGET (NI,FC,GC,XX,YY,BFI,CDFX,POI,POJ,EPS3,IFLAG,2) IF (NJ.LE.1.OR.IFLAG.NE.0) RETURN C C--- COMPUTE BETA C.D.F. BY RECURRENCE WHEN I>IMIN, J>JMIN C DO 20 I = 2, NI BFJ(1) = BFI(I) DO 10 J = 2, NJ BFJ(J) = XX*BFJ(J)+YY*BFJ(J-1) CDFX = CDFX+POI(I)*POJ(J)*BFJ(J) 10 CONTINUE 20 CONTINUE RETURN END *CDFF SUBROUTINE CDFF (X,DF1,DF2,EPS,IFLAG,CDFX) C C----------------------------------------------------------------------- C CDFF WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING THE CUMULATIVE DISTRIBUTION FUNCTION OF THE F C DISTRIBUTION FROM THE CUMULATIVE DISTRIBUTION FUNCTION OF C THE BETA DISTRIBUTION. THE RELATIONSHIP BETWEEN THE TWO C DISTRIBUTIONS IS GIVEN IN THE REFERENCE BELOW. C C SUBPROGRAMS CALLED: CDFBET (BETA C.D.F.) C C CURRENT VERSION COMPLETED AUGUST 15, 1987 C C REFERENCE: REEVE, CHARLES P., 'AN ALGORITHM FOR COMPUTING THE BETA C C.D.F. TO A SPECIFIED ACCURACY', STATISTICAL ENGINEERING C DIVISION NOTE 86-3, OCTOBER 1986. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * X = VALUE AT WHICH THE C.D.F. IS TO BE COMPUTED (REAL) C C * DF1 = FIRST 'DEGREES OF FREEDOM' PARAMETER (>0) (REAL) C C * DF2 = SECOND 'DEGREES OF FREEDOM' PARAMETER (>0) (REAL) C C * EPS = THE DESIRED ABSOLUTE ACCURACY OF THE C.D.F. (>0) (REAL) C C IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION: C 0 -> NO ERRORS DETECTED C 1,2 -> ERROR INDICATORS FROM THE BETA C.D.F. ROUTINE C 3 -> DF1 AND/OR DF2 IS NON-POSITIVE C C CDFX = THE C.D.F. EVALUATED AT X (REAL) C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C CDFX = 0.0 C C--- CHECK FOR NON-POSITIVE DEGREES OF FREEDOM C IF (AMIN1(DF1,DF2).LE.0.0) THEN IFLAG = 3 RETURN C ENDIF H = DF1*X Y = H/(H+DF2) P = 0.5*DF1 Q = 0.5*DF2 CALL CDFBET (Y,P,Q,EPS,IFLAG,CDFX) RETURN C END *CDFGAM SUBROUTINE CDFGAM (X,ALPHA,EPS,IFLAG,CDFX) C C----------------------------------------------------------------------- C CDFGAM WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MD 20899 C C FOR: COMPUTING THE CUMULATIVE DISTRIBUTION FUNCTION OF THE GAMMA C DISTRIBUTION (ALSO KNOWN AS THE INCOMPLETE GAMMA RATIO) TO A C SPECIFIED ACCURACY (TRUNCATION ERROR IN THE INFINITE SERIES). C THE ALGORITHM, DESCRIBED IN REFERENCE 2, IS A MODIFICATION OF C THE ALGORITHM OF REFERENCE 1. THREE FEATURES HAVE BEEN ADDED: C C 1) A PRECISE METHOD OF MEETING THE TRUNCATION ACCURACY, C 2) COMPUTATION OF THE UPPER TAIL AREA BY DECREMENTING ALPHA C WHEN THAT METHOD IS MORE EFFICIENT, AND C 3) A CONSTANT UFLO >= THE UNDERFLOW LIMIT ON THE COMPUTER. C C SUBPROGRAMS CALLED: DGAMLN (LOG OF GAMMA FUNCTION) C C CURRENT VERSION COMPLETED OCTOBER 29, 1986 C C REFERENCES: C C 1) LAU, CHI-LEUNG, 'A SIMPLE SERIES FOR THE INCOMPLETE GAMMA C INTEGRAL', ALGORITHM AS 147, APPLIED STATISTICS, VOL. 29, C NO. 1, 1980, PP. 113-114. C C 2) REEVE, CHARLES P., 'AN ALGORITHM FOR COMPUTING THE GAMMA C.D.F. C TO A SPECIFIED ACCURACY', STATISTICAL ENGINEERING DIVISION C NOTE 86-2, OCTOBER 1986. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * X = VALUE AT WHICH THE C.D.F IS TO BE COMPUTED (REAL) C C * ALPHA = PARAMETER OF THE GAMMA FUNCTION (>0) (REAL) C C * EPS = THE DESIRED ABSOLUTE ACCURACY OF THE C.D.F (>0) (REAL) C C IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION: C 0 -> NO ERRORS DETECTED C 1 -> EITHER ALPHA OR EPS IS <= UFLO C 2 -> NUMBER OF TERMS EVALUATED IN THE INFINITE SERIES C EXCEEDS IMAX. C C CDFX = THE C.D.F. EVALUATED AT X (REAL) C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C LOGICAL LL DOUBLE PRECISION DX,DGAMLN DATA IMAX,UFLO / 5000,1.0E-100 / CDFX = 0.0 C C--- CHECK FOR VALIDITY OF ARGUMENTS ALPHA AND EPS C IF (ALPHA.LE.UFLO.OR.EPS.LE.UFLO) THEN IFLAG = 1 RETURN ENDIF IFLAG = 0 C C--- CHECK FOR SPECIAL CASE OF X C IF (X.LE.0.0) RETURN C C--- EVALUATE THE GAMMA P.D.F. AND CHECK FOR UNDERFLOW C DX = DBLE(X) PDFL = SNGL(DBLE(ALPHA-1.0)*DLOG(DX)-DX-DGAMLN(ALPHA)) IF (PDFL.LT.ALOG(UFLO)) THEN IF (X.GE.ALPHA) CDFX = 1.0 ELSE P = ALPHA U = EXP(PDFL) C C--- DETERMINE WHETHER TO INCREMENT OR DECREMENT ALPHA (A.K.A. P) C LL = .TRUE. IF (X.GE.P) THEN K = INT(P) IF (P.LE.REAL(K)) K = K-1 ETA = P-REAL(K) BL = SNGL(DBLE(ETA-1.0)*DLOG(DX)-DX-DGAMLN(ETA)) LL = BL.GT.ALOG(EPS) ENDIF EPSX=EPS/X IF (LL) THEN C C--- INCREMENT P C DO 10 I = 0, IMAX IF (U.LE.EPSX*(P-X)) RETURN U = X*U/P CDFX = CDFX+U P = P+1.0 10 CONTINUE IFLAG = 2 ELSE C C--- DECREMENT P C DO 20 J = 1, K P = P-1.0 IF (U.LE.EPSX*(X-P)) GO TO 30 CDFX = CDFX+U U = P*U/X 20 CONTINUE 30 CDFX = 1.0-CDFX ENDIF ENDIF RETURN END *CDFNOR SUBROUTINE CDFNOR (Z,EPS,IFLAG,CDFZ) C C----------------------------------------------------------------------- C CDFNOR WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: COMPUTING THE CUMULATIVE DISTRIBUTION FUNCTION OF THE STANDARD C NORMAL DISTRIBUTION TO A SPECIFIED ACCURACY. THE C.D.F. OF C THE GAMMA DISTRIBUTION IS FIRST COMPUTED, THEN TRANSFORMED TO C THE C.D.F. OF THE NORMAL DISTRIBUTION. C C NOTE: TIMING TESTS ON THE CYBER 180/855 COMPUTER AT N.I.S.T. C INDICATE THAT THE AVERAGE CPU TIME FOR COMPUTING ONE C.D.F. C IS ABOUT 0.0004 SECOND. C C SUBPROGRAMS CALLED: CDFGAM (C.D.F. OF THE GAMMA DISTRIBUTION) C C CURRENT VERSION COMPLETED APRIL 7, 1989 C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * Z = THE VALUE FOR WHICH THE NORMAL C.D.F. IS TO BE COMPUTED C [REAL] C C * EPS = THE ABSOLUTE ACCURACY REQUIREMENT FOR THE C.D.F. [REAL] C C IFLAG = ERROR INDICATOR ON OUTPUT [INTEGER] INTERPRETATION: C 0 -> NO ERRORS DETECTED. C 1 -> ERROR FLAG FROM CDFGAM C 2 -> ERROR FLAG FROM CDFGAM C C CDFZ = THE C.D.F. OF THE STANDARD NORMAL DISTRIBUTION EVALUATED C AT X [REAL] C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DEL = 2.0*EPS IF (Z.EQ.0.0) THEN CDFZ = 0.0 ELSE X = 0.5*Z*Z CALL CDFGAM (X,0.5,DEL,IFLAG,CDFX) IF (IFLAG.NE.0) RETURN IF (Z.GT.0.0) THEN CDFZ = 0.5+0.5*CDFX ELSE CDFZ = 0.5-0.5*CDFX ENDIF ENDIF RETURN C END *CDFT SUBROUTINE CDFT (X,DF,EPS,IFLAG,CDFX) C C----------------------------------------------------------------------- C CDFT WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING THE CUMULATIVE DISTRIBUTION FUNCTION OF THE T C DISTRIBUTION FROM THE CUMULATIVE DISTRIBUTION FUNCTION OF C THE BETA DISTRIBUTION. THE RELATIONSHIP BETWEEN THE TWO C DISTRIBUTIONS IS GIVEN IN THE REFERENCE BELOW. C C SUBPROGRAMS CALLED: CDFBET (BETA C.D.F.) C C CURRENT VERSION COMPLETED AUGUST 15, 1987 C C REFERENCE: REEVE, CHARLES P., 'AN ALGORITHM FOR COMPUTING THE BETA C C.D.F. TO A SPECIFIED ACCURACY', STATISTICAL ENGINEERING C DIVISION NOTE 86-3, OCTOBER 1986. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * X = VALUE AT WHICH THE C.D.F. IS TO BE COMPUTED (REAL) C C * DF = 'DEGREES OF FREEDOM' PARAMETER (>0) (REAL) C C * EPS = THE DESIRED ABSOLUTE ACCURACY OF THE C.D.F. (>0) (REAL) C C IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION: C 0 -> NO ERRORS DETECTED C 1,2 -> ERROR INDICATORS FROM THE BETA C.D.F. ROUTINE C 3 -> DF IS NON-POSITIVE C C CDFX = THE C.D.F. EVALUATED AT X (REAL) C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C CDFX = 0.0 C C--- CHECK FOR NON-POSITIVE DEGREES OF FREEDOM C IF (DF.LE.0.0) THEN IFLAG = 3 RETURN C ENDIF H = X*X Y = H/(H+DF) P = 0.5 Q = 0.5*DF CALL CDFBET (Y,P,Q,EPS,IFLAG,CDFY) CDFX = 0.5*(1.0+SIGN(CDFY,X)) RETURN C END *CENSCL SUBROUTINE CENSCL (X,W,N1,N2,IOPT,ITER,P1TAIL,CENTER,SCALE) C C * * * * * * * * * * CPR242 * * * * * * * * * * * C CENSCL WRITTEN BY CHARLES P. REEVE C FOR: COMPUTING STATISTICAL ESTIMATES OF THE CENTER AND SCALE C OF ENTRIES N1 THROUGH N2 (INCLUSIVE) IN A SET OF DATA. C CLASSICAL OR ROBUST ESTIMATES MAY BE SPECIFIED. C EXTERNAL SUBPROGRAMS CALLED: SORT1 C INTERNAL SUBPROGRAMS CALLED: MAD, MEANSD, MEDIAN, VCOPY C CURRENT VERSION COMPLETED MARCH 2, 1983 C * * * * * * * * * * * * * * * * * * * * * * * C C DEFINITION OF PASSED PARAMETERS: C C * X(*) = INPUT: VECTOR (DIMENSION >= N2) CONTAINING DATA C OUTPUT: NO CHANGE C C W(*) = WORKSPACE VECTOR (DIMENSION >= N2) C C * N1 = FIRST ENTRY OF INTEREST IN DATA VECTOR X C C * N2 = LAST ENTRY OF INTEREST IN DATA VECTOR X C C * IOPT = 1 TO COMPUTE MEAN AND STANDARD DEVIATION C 2 TO COMPUTE MEDIAN AND MEDIAN ABSOLUTE DEVIATION C 3 TO COMPUTE TRIMMED MEAN AND STANDARD DEVIATION C 4 TO COMPUTE WINSORIZED MEAN AND STANDARD DEVIATION C 5 TO COMPUTE M(T) AND S(T) (REEVE) C 6 TO COMPUTE M(R) AND S(R) (REEVE) C 7 TO COMPUTE BIWEIGHT MEAN AND STD. DEV. (GROSS) C C * ITER = MAXIMUM NUMBER OF ITERATIONS TO BE USED IN COMPUTATION C OF BIWEIGHT ESTIMATE OF LOCATION. A POSITIVE VALUE WILL C CAUSE A DIAGNOSTIC MESSAGE TO BE PRINTED WHEN 'ITER' C ITERATIONS HAVE OCCURRED. A NEGATIVE VALUE WILL CAUSE C NO DIAGNOSTIC MESSAGE TO BE PRINTED. IN EITHER CASE, THE C LATEST COMPUTED VALUE OF THE BIWEIGHT ESTIMATE IS RETURNED C WHEN 'ABS(ITER)' ITERATIONS HAVE OCCURRED. CONVERGENCE TO C 'X(K)' WILL OCCUR EARLIER ON ITERATION 'K' IF C C ABS( X(K)-X(K-1) ) < 0.0001 * ABS( X(K) ) + .000001 C C C * P1TAIL = PROPORTION OF ORDERED DATA AT EACH EXTREME TO BE C CONSIDERED THE TAIL REGION (USED IF IOPT=3 OR 4) C C CENTER = MEAN IF IOPT=1 C MEDIAN IF IOPT=2 C TRIMMED MEAN IF IOPT=3 C WINSORIZED MEAN IF IOPT=4 C M(T) (REEVE) IF IOPT=5 C M(R) (REEVE) IF IOPT=6 C BIWEIGHT MEAN IF IOPT=7 C C SCALE = STANDARD DEVIATION IF IOPT=1 C MEDIAN ABSOLUTE DEVIATION IF IOPT=2 C TRIMMED STANDARD DEVIATION IF IOPT=3 C WINSORIZED STANDARD DEVIATION IF IOPT=4 C S(T) (REEVE) IF IOPT=5 C S(R) (REEVE) IF IOPT=6 C BIWEIGHT STD. DEV. (GROSS) IF IOPT=7 C C * DENOTES VARIABLES REQUIRING INPUT VALUES C C * * * * * * * * * * * * * * * * * * * * * * * C DIMENSION X(*),W(*) C C--- COMPUTE NUMBER OF DATA POINTS USED C N = N2-N1+1 C C--- CHECK FOR AT LEAST FOUR (4) DATA POINTS -- OTHERWISE, STOP! C IF (N.GE.4) GO TO 10 PRINT *,' *** NUMBER OF DATA POINTS IS < 4' PRINT *,' *** EXECUTION STOPPED IN SUBROUTINE CENSCL' STOP C C--- SET ACCIDENTAL NEGATIVE TAIL PROPORTION TO ZERO C 10 P1TAIL = AMAX1(P1TAIL,0.) C C--- BRANCH TO CLASSICAL OR ROBUST ESTIMATES C GO TO (20,30,40,50,70,70,130), IOPT C C**************************************************************** C*** COMPUTE CLASSICAL ESTIMATES: MEAN AND STANDARD DEVIATION *** C**************************************************************** C 20 CALL MEANSD (X,N1,N2,CENTER,SCALE) RETURN C C********************************************************************** C*** COMPUTE ROBUST ESTIMATES: MEDIAN AND MEDIAN ABSOLUTE DEVIATION *** C********************************************************************** C 30 CALL VCOPY (X,W,N1,N2) CALL MEDIAN (W,N1,N2,CENTER) CALL MAD (W,N1,N2,CENTER,SCALE) RETURN C C********************************************************************* C*** COMPUTE ROBUST ESTIMATES: TRIMMED MEAN AND STANDARD DEVIATION *** C********************************************************************* C 40 CALL VCOPY (X,W,N1,N2) CALL SORT1 (W,N1,N2) NTRIM = MIN0((N-1)/2,NINT(P1TAIL*REAL(N))) NLO = N1+NTRIM NHI = N2-NTRIM CALL MEANSD (W,NLO,NHI,CENTER,SCALE) RETURN C C*************************************************************** C*** COMPUTE ROBUST ESTIMATES: WINSORIZED MEAN AND STD. DEV. *** C*************************************************************** C 50 CALL VCOPY (X,W,N1,N2) CALL SORT1 (W,N1,N2) NTRIM = MIN0((N-1)/2,NINT(P1TAIL*REAL(N))) DO 60 I = 1, NTRIM W(N1+I-1) = W(N1+NTRIM) W(N2-I+1) = W(N2-NTRIM) 60 CONTINUE CALL MEANSD (W,N1,N2,CENTER,SCALE) RETURN C C**************************************************************** C*** COMPUTE ROBUST ESTIMATES: M(T) AND S(T) OR M(R) AND S(R) *** C**************************************************************** C 70 CALL VCOPY (X,W,N1,N2) CALL SORT1 (W,N1,N2) C C--- SET CONSTANT WHICH RELATES SHORTEST INTERVAL TO STANDARD C--- DEVIATION OF THE NORMAL DISTRIBUTION C CNORM = 1.348979 C C--- SET REJECTION LIMITS AS A MULTIPLE OF THE STANDARD DEVIATION C CRHO = 3.0 RMIN = 1E30 L = N/2 NR = (N-1)/2 B = REAL(L-NR+1) A = 4.-B C C--- FIND THE SHORTEST INTERVAL COVERING HALF THE POINTS C DO 80 I = 1, NR J = N1+I-1 K = J+L R = (A*(W(K+1)-W(J))+B*(W(K)-W(J+1)))/4. IF (R.GE.RMIN) GO TO 80 RMIN = R K1 = J 80 CONTINUE K2 = K1+1 K3 = K1+L K4 = K3+1 CENTER = (A*(W(K1)+W(K4))+B*(W(K2)+W(K3)))/8. SCALE = RMIN/CNORM C C--- RETURN M(R) AND S(R) IF IOPT=6 C IF (IOPT.EQ.6) RETURN BOUND = CRHO*SCALE C C--- TRIM OUTLYING OBSERVATIONS AT LOWER END C NLO = N1 90 IF (CENTER-W(NLO).LE.BOUND) GO TO 100 NLO = NLO+1 GO TO 90 C 100 NHI = N2 C C--- TRIM OUTLYING OBSERVATIONS AT UPPER END C 110 IF (W(NHI)-CENTER.LE.BOUND) GO TO 120 NHI = NHI-1 GO TO 110 C 120 CALL MEANSD (W,NLO,NHI,CENTER,SCALE) C C--- RETURN M(T) AND S(T) FOR IOPT=5 C RETURN C C********************************************************************** C*** COMPUTE ROBUST ESTIMATES: BIWEIGHT MEAN AND STANDARD DEVIATION *** C********************************************************************** C--- DEFINE CONVERGENCE CONSTANTS C 130 C1 = 1E-4 C2 = 1E-6 ITERA = IABS(ITER) C C--- COUNT NUMBER OF ITERATIONS FOR BIWEIGHT MEAN COMPUTATION C KOUNT = 0 CALL VCOPY (X,W,N1,N2) CALL MEDIAN (W,N1,N2,CENTER) 140 CALL MAD (W,N1,N2,CENTER,SCALE) S1 = 0. S2 = 0. DO 150 I = N1, N2 U = (W(I)-CENTER)/(6.*SCALE) IF (ABS(U).GE.1.0) GO TO 150 U = (1.-U*U)**2 S1 = S1+U*W(I) S2 = S2+U 150 CONTINUE COLD = CENTER CENTER = S1/S2 IF (ABS(CENTER-COLD).LT.C1*ABS(CENTER)+C2) GO TO 160 KOUNT = KOUNT+1 IF (KOUNT.LE.ITERA) GO TO 140 IF (ITER.LT.0) GO TO 160 N = N2-N1+1 WRITE (6,1010) ITER,COLD,CENTER,N 160 S1 = 0. S2 = 0. DO 170 I = N1, N2 U = (W(I)-CENTER)/(9.*SCALE) IF (ABS(U).GE.1.0) GO TO 170 V = 1.-U*U S1 = S1+(W(I)-CENTER)**2*V**4 S2 = S2+V*(1.-5.*U*U) 170 CONTINUE SCALE = SQRT(REAL(N)*S1/(S2*(S2-1.))) RETURN 1010 FORMAT (1X,'*** ',I3, * ' ITERATION LIMIT ON BIWEIGHT *** OLD CENTER =',G12.5,3X, * 'NEW CENTER =',G12.5,3X,'SAMPLE SIZE =',I6) C END C C SUBROUTINE MEANSD (X,NLO,NHI,AMEAN,SD) C C--- COMPUTE CLASSICAL MEAN AND STANDARD DEVIATION OF DATA BETWEEN C--- POSITIONS X(NLO) AND X(NHI) INCLUSIVE C DIMENSION X(*) AMEAN = 0. DO 10 I = NLO, NHI AMEAN = AMEAN+X(I) 10 CONTINUE AMEAN = AMEAN/REAL(NHI-NLO+1) SD = 0. DO 20 I = NLO, NHI SD = SD+(X(I)-AMEAN)**2 20 CONTINUE SD = SQRT(SD/REAL(NHI-NLO)) RETURN C END C C SUBROUTINE VCOPY (A,B,NLO,NHI) C C--- COPY SEGMENT OF VECTOR A INTO SAME SEGMENT OF VECTOR B C DIMENSION A(*),B(*) DO 10 I = NLO, NHI B(I) = A(I) 10 CONTINUE RETURN C END C C SUBROUTINE MEDIAN (X,NLO,NHI,AMED) C C--- COMPUTE MEDIAN OF DATA BETWEEN POSITIONS X(NLO) AND X(NHI) C--- INCLUSIVE. DATA IN THIS REGION OF VECTOR X ARE ALTERED. C DIMENSION X(*) C C--- SORT REGION OF INTEREST IN VECTOR X C CALL SORT1 (X,NLO,NHI) C C--- COMPUTE MEDIAN C I = (NLO+NHI)/2 J = (NLO+NHI+1)/2 AMED = (X(I)+X(J))/2. RETURN C END C C SUBROUTINE MAD (X,NLO,NHI,C,AMAD) C C--- COMPUTE MEDIAN ABSOLUTE DEVIATION (MAD) OF X FROM C, AN ESTIMATE C--- OF THE CENTER OF DATA, BETWEEN X(NLO) AND X(NHI) INCLUSIVE. VECTOR C--- X IS EXPECTED TO BE SORTED AND IS UNCHANGED BY THIS SUBROUTINE. C--- NOTE: IF THE NUMBER OF ENTRIES OF INTEREST, N, IS EVEN THE MAD IS C--- THE N/2 LARGEST DEVIATION, A SLIGHT OVERESTIMATE. C DIMENSION X(*) MLO = NLO MHI = NHI K = (MHI-MLO)/2 DO 20 I = 1, K IF (X(MHI)+X(MLO).GT.2.0*C) GO TO 10 MLO = MLO+1 GO TO 20 C 10 MHI = MHI-1 20 CONTINUE AMAD = MAX(ABS(X(MHI)-C),ABS(X(MLO)-C)) RETURN C END *CIELIP SUBROUTINE CIELIP (X,Y,N,P,YP,M,A,B,XBAR,S,XP,XL,XU,IOPT) C C----------------------------------------------------------------------- C CIELIP WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, MD C C FOR: INVERSE PREDICTION FOR A STRAIGHT LINE FIT TO DATA. THE C DEPENDENT VARIABLE IS Y (WITH ERROR) AND THE INDEPENDENT C VARIABLE IS X (WITHOUT ERROR). GIVEN N PAIRS (X,Y) A C LEAST SQUARES FIT OF THE MODEL C C Y(I) = ALPHA + BETA * (X - MEAN(X)) + ERROR C C IS MADE. FOR A USER SPECIFIED VALUE Y', (THE MEAN OF M C MEASUREMENTS WITH THE SAME VARIANCE AS THE Y-VALUES) C A PREDICTED X-VALUE, X', IS COMPUTED ALONG WITH ITS C 100(1-P)% CONFIDENCE INTERVAL. THE METHOD IS DUE TO C CHURCHILL EISENHART [SEE REFERENCE 1 BELOW, ESPECIALLY C EQUATIONS (25) AND (26)]. THE CONFIDENCE INTERVAL IS C ASYMMETRIC ABOUT X'. C C NOTE: THE CONFIDENCE INTERVAL DEFINED BY XL AND XU SHOULD BE C INTERPRETED AS FOLLOWS: C C 1) IF XL < XU THEN IT IS THE REGION BETWEEN XL AND XU. C C 2) IF XL > XU THEN IT IS THE UNION OF THE TWO C SEMI-INFINITE RAYS (-INF,XU) AND (XL,+INF). C C 3) IF XL = XU = 0 THEN NO CONFIDENCE INTERVAL EXISTS. C THIS IS CAUSED BY HIGHLY SCATTERED DATA WHICH IS OF C DUBIOUS VALUE FOR INVERSE PREDICTION. C C SUBPROGRAMS CALLED: MDSTI (IMSL) - T PERCENT POINT FUNCTION C C CURRENT VERSION COMPLETED JULY 13, 1984 C C REFERENCES: C C 1. EISENHART, CHURCHILL, 'THE INTERPRETATION OF CERTAIN REGRESSION C METHODS AND THEIR USE IN BIOLOGICAL AND INDUSTRIAL RESEARCH', C THE ANNALS OF MATHEMATICAL STATISTICS, VOLUME X, NO. 2, JUNE C 1939, PP. 162-186. C C 2. REEVE, CHARLES P., 'A COMPARISON OF SOME INVERSE PREDICTION C METHODS FOR STRAIGHT LINE DATA', STATISTICAL ENGINEERING C DIVISION NOTE 84-1, JUNE 1984. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * X(*) = VECTOR (LENGTH N) OF INDEPENDENT VARIABLE C C * Y(*) = VECTOR (LENGTH N) OF DEPENDENT VARIABLE C C * N = NUMBER OF DATA VALUES C C * P = SIGNIFICANCE LEVEL (0.05 FOR 95% C.I., ETC.) C C * YP = Y-VALUE FOR WHICH PREDICTED X-VALUE IS TO BE COMPUTED C C * M = NUMBER OF MEASUREMENTS INVOLVED IN OBTAINING YP. IF C YP IS EXACT, SET M=0 C C A = ESTIMATE OF ALPHA (CONSTANT TERM) C C B = ESTIMATE OF BETA (SLOPE) C C XBAR = MEAN OF THE N X-VALUES C C S = RESIDUAL STANDARD DEVIATION C C XP = PREDICTED X-VALUE CORRESPONDING TO YP C C XL = LOWER BOUND OF CONFIDENCE INTERVAL FOR XP C C XU = UPPER BOUND OF CONFIDENCE INTERVAL FOR XP C C * IOPT = 0 IF THE QUANTITIES A, B, XBAR, AND S ARE TO BE COMPUTED C FROM A LEAST SQUARES FIT USING INPUT VALUES OF X(*), C Y(*), AND N. C = 1 IF INPUT VALUES OF THE QUANTITIES A, B, XBAR, S, AND C N ARE TO BE USED C C * DENOTES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION X(*),Y(*) RN = REAL(N) C C--- BRANCH ACCORDING TO IOPT C IF (IOPT.GT.0) GO TO 40 C C--- COMPUTE LEAST SQUARES ESTIMATES OF COEFFICIENTS C SX = 0.0 SY = 0.0 DO 10 I = 1, N SX = SX+X(I) SY = SY+Y(I) 10 CONTINUE A = SY/RN XBAR = SX/RN SXX = 0.0 SXY = 0.0 DO 20 I = 1, N SXX = SXX+(X(I)-XBAR)**2 SXY = SXY+(X(I)-XBAR)*Y(I) 20 CONTINUE R = 1.0/SXX B = R*SXY C C--- COMPUTE RESIDUAL STANDARD DEVIATION C S = 0.0 DO 30 I = 1, N S = S+(Y(I)-A-B*(X(I)-XBAR))**2 30 CONTINUE DF = REAL(N-2) S = SQRT(S/DF) C C--- COMPUTE T CRITICAL VALUE C CALAN CALL MDSTI (P,DF,T,IER) C C--- COMPUTE INVERSE OF M (UNLESS M=0) C 40 IF (M.EQ.0) THEN AIM = 0.0 ELSE AIM = 1.0/ABS(REAL(M)) ENDIF C C--- COMPUTE PREDICTED X-VALUE FOR YP (SET XP = 0 IF SLOPE = 0) C IF (B.EQ.0.0) THEN XP = 0.0 ELSE XP = XBAR+(YP-A)/B ENDIF C C--- COMPUTE CONFIDENCE INTERVAL FOR PREDICTED X-VALUE C Q = B**2-R*(T*S)**2 RADICL = R*(YP-A)**2+(AIM+1.0/RN)*Q IF (RADICL.LT.0.0) THEN XU = 0.0 XL = 0.0 ELSE TERM1 = XBAR+B*(YP-A)/Q TERM2 = T*S*SQRT(RADICL)/Q XL = TERM1-TERM2 XU = TERM1+TERM2 ENDIF RETURN C END *DGAMLN DOUBLE PRECISION FUNCTION DGAMLN (X) C C----------------------------------------------------------------------- C DGAMLN WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING THE DOUBLE PRECISION LOG OF THE GAMMA FUNCTION WITH C SINGLE PRECISION PARAMETER X>0. THE MAXIMUM TRUNCATION ERROR C IN THE INFINITE SERIES (SEE REFERENCE 1) IS DETERMINED BY THE C CONSTANT XMIN. WHEN X0 (X=0, M>=0, AND M<=N. THAT C VALUE IS RETURNED IN DNCMLN. C C SUBPROGRAMS CALLED: DGAMLN (NATURAL LOGARITHM OF GAMMA FUNCTION) C C CURRENT VERSION COMPLETED MAY 1, 1989 C----------------------------------------------------------------------- C DOUBLE PRECISION D1,D2,D3,DGAMLN IF (N.LT.0) STOP '*** STOP: N<0 IN SUBROUTINE DNCMLN ***' IF (M.LT.0) STOP '*** STOP: M<0 IN SUBROUTINE DNCMLN ***' IF (M.GT.N) STOP '*** STOP: M>N IN SUBROUTINE DNCMLN ***' R1 = REAL(N+1) R2 = REAL(M+1) R3 = REAL(N-M+1) D1 = DGAMLN(R1) D2 = DGAMLN(R2) D3 = DGAMLN(R3) DNCMLN = D1+D2-D3 RETURN C END *DPR1LN DOUBLE PRECISION FUNCTION DPR1LN (N1,M1,N2,M2) C C----------------------------------------------------------------------- C DPR1LN WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: COMPUTING THE DOUBLE PRECISION NATURAL LOGARITHM OF C C N1!/[M1!(N1-M1)!] N2!/[M2!(N2-M2)!] C ----------------------------------- C (N1+N2)!/[(M1+M2)!(N1+N2-M1-M2)!] C C WHERE N1>=0, M1>=0, M1<=N1, N2>=0, M2>=0, AND M2<=N2. THAT C VALUE IS RETURNED IN DPR1LN. C C SUBPROGRAMS CALLED: DNCMLN (NATURAL LOGARITHM OF 'N CHOOSE M') C C CURRENT VERSION COMPLETED MAY 1, 1989 C----------------------------------------------------------------------- C DOUBLE PRECISION D1,D2,D3,DNCMLN IF (N1.LT.0) STOP '*** STOP: N1<0 IN SUBROUTINE DPR1LN ***' IF (M1.LT.0) STOP '*** STOP: M1<0 IN SUBROUTINE DPR1LN ***' IF (M1.GT.N1) STOP '*** STOP: M1>N1 IN SUBROUTINE DPR1LN ***' IF (N2.LT.0) STOP '*** STOP: N2<0 IN SUBROUTINE DPR1LN ***' IF (M2.LT.0) STOP '*** STOP: M2<0 IN SUBROUTINE DPR1LN ***' IF (M2.GT.N2) STOP '*** STOP: M2>N2 IN SUBROUTINE DPR1LN ***' D1 = DNCMLN(N1,M1) D2 = DNCMLN(N2,M2) D3 = DNCMLN(N1+N2,M1+M2) DPR1LN = D1+D2-D3 RETURN C END *DWESD SUBROUTINE DWESD (X,LDX,N,M,U,LDU,V,LDV,E,SD) C C----------------------------------------------------------------------- C DWESD WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING THE EXPECTED VALUE AND VARIANCE OF THE DURBIN-WATSON C STATISTIC D AS APPLIED TO THE RESIDUALS FROM THE LEAST SQUARES C FIT OF THE MODEL Y = XB + ERROR. FOR A GIVEN MATRIX X, EXACT C FORMULAS FOR E(D) AND VAR(D) ARE GIVEN IN THE REFERENCES BELOW. C C SUBPROGRAMS CALLED: MATXXI (STSPAC) C XTAX, XTA2X, TRACES (ATTACHED) C C CURRENT VERSION COMPLETED NOVEMBER 6, 1985 C C REFERENCES: C C 1. DURBIN, J. AND WATSON, G.S., 'TESTING FOR SERIAL CORRELATION IN C LEAST SQUARES REGRESSION. I', BIOMETRIKA 37, 1950, PP. 409-428 C C 2. DURBIN, J. AND WATSON, G.S., 'TESTING FOR SERIAL CORRELATION IN C LEAST SQUARES REGRESSION. II', BIOMETRIKA 38, 1951, PP. 159-178 C C 3. DURBIN, J. AND WATSON, G.S., 'TESTING FOR SERIAL CORRELATION IN C LEAST SQUARES REGRESSION. III', BIOMETRIKA 58, 1971, PP. 1-19 C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * X = MATRIX (SIZE NXM) IN THE MODEL Y = XB + ERROR C C * LDX = LEADING DIMENSION OF X (LDX >=N+2) C C * N = NUMBER OF EQUATIONS AND NUMBER OF ROWS OF X C C * M = NUMBER OF UNKNOWNS AND NUMBER OF COLUMNS OF X C C * U = MATRIX (SIZE MXM) FOR INTERMEDIATE COMPUTATIONS C C * LDU = LEADING DIMENSION OF U (LDU >= M) C C * V = MATRIX (SIZE MXM) FOR INTERMEDIATE COMPUTATIONS C C * LDV = LEADING DIMENSION OF V (LDV >= M) C C E = EXPECTED VALUE OF DURBIN-WATSON STATISTIC FOR X C C SD = STANDARD DEVIATION OF DURBIN-WATSON STATISTIC FOR X C C * DENOTES VARIABLES REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION X(LDX,*),U(LDU,*),V(LDV,*) C C--- CHECK FOR PARAMETERS OUT OF RANGE C IF (M.GE.N.OR.N.LT.5) THEN WRITE (6,1000) STOP C ENDIF C C--- COMPUTE U = X'AX C CALL XTAX (X,LDX,U,LDU,N,M) C C--- COMPUTE V = X'AAX C CALL XTA2X (X,LDX,V,LDV,N,M) C C--- COMPUTE W = INV(X'X) (STORE IN X) C CALL MATXXI (X,LDX,N,M) C C--- COMPUTE TRACE(UW) AND TRACE(UWUW) C CALL TRACES (U,LDU,X,LDX,M,M,TR1,TR2,2) C C--- COMPUTE TRACE(VW) C CALL TRACES (V,LDV,X,LDX,M,M,TR3,DUM,1) C C--- COMPUTE CONSTANTS P AND Q C P = REAL(2*N-2)-TR1 Q = REAL(6*N-8)-2.0*TR3+TR2 C C--- COMPUTE EXPECTED VALUE AND STANDARD DEVIATION OF D C E = P/REAL(N-M) SD = SQRT(2.0*(Q-P*E)/REAL((N-M)*(N-M+2))) RETURN C C 1000 FORMAT (/1X,'*** FATAL ERROR -- MUST HAVE N>4 AND N>M IN DWESD'/) C END C SUBROUTINE XTAX (X,LDX,Z,LDZ,N,M) C C----------------------------------------------------------------------- C XTAX WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING X'AX WHERE X IS NXM AND A IS AN NXN TRIDIAGONAL C MATRIX ASSOCIATED WITH THE DURBIN-WATSON STATISTIC. MATRIX C A HAS THE FORM C C 1 -1 0 0 ... 0 0 C -1 2 -1 0 ... 0 0 C 0 -1 2 -1 ... 0 0 C . . C . . C . . C 0 0 0 0 ... 2 -1 C 0 0 0 0 ... -1 1 . C C THE COMPUTATION TAKES THE FORM Z = (L'X)'(L'X) WHERE L IS C THE CHOLESKY FACTORIZATION OF A (A=LL'). C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED NOVEMBER 6, 1985 C----------------------------------------------------------------------- C DIMENSION X(LDX,*),Z(LDZ,*) DO 30 I = 1, M DO 20 J = 1, I S = 0.0 DO 10 K = 2, N S = S+(X(K-1,I)-X(K,I))*(X(K-1,J)-X(K,J)) 10 CONTINUE Z(I,J) = S Z(J,I) = S 20 CONTINUE 30 CONTINUE RETURN C END C SUBROUTINE XTA2X (X,LDX,Z,LDZ,N,M) C C----------------------------------------------------------------------- C XTA2X WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING X'AAX WHERE X IS NXM AND AA IS AN NXN PENTADIAGONAL C MATRIX ASSOCIATED WITH THE DURBIN-WATSON STATISTIC. MATRIX C AA HAS THE FORM C C 2 -3 1 0 0 ... 0 0 0 C -3 6 -4 1 0 ... 0 0 0 C 1 -4 6 -4 1 ... 0 0 0 C 0 1 -4 6 -4 ... 0 0 0 C . . C . . C . . C 0 0 0 0 0 ... 6 -4 1 C 0 0 0 0 0 ... -4 6 -3 C 0 0 0 0 0 ... 1 -3 2 . C C THE COMPUTATION TAKES THE FORM Z = (AX)'(AX) SINCE A=A'. C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED NOVEMBER 6, 1985 C----------------------------------------------------------------------- C DIMENSION X(LDX,*),Z(LDZ,*) DO 30 I = 1, M DO 20 J = 1, I S = (X(1,I)-X(2,I))*(X(1,J)-X(2,J))+(X(N,I)-X(N-1,I))*(X(N,J * )-X(N-1,J)) DO 10 K = 2, N-1 S = S+(-X(K-1,I)+2.0*X(K,I)-X(K+1,I))*(-X(K-1,J)+2.0*X(K, * J)-X(K+1,J)) 10 CONTINUE Z(I,J) = S Z(J,I) = S 20 CONTINUE 30 CONTINUE RETURN C END C SUBROUTINE TRACES (U,LDU,V,LDV,N,M,TR1,TR2,IND) C C----------------------------------------------------------------------- C TRACES WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING THE TRACES OF THE NXN MATRICES UV AND UVUV WHERE C U IS NXM AND V IS MXN. NOTE THAT NEITHER UV NOR UVUV C ACTUALLY NEEDS TO BE COMPUTED. C C NOTE: IF IND=1 THEN ONLY TRACE(UV) IS COMPUTED C IF IND=2 THEN BOTH TRACE(UV) AND TRACE(UVUV) ARE COMPUTED C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED NOVEMBER 5, 1985 C----------------------------------------------------------------------- C DIMENSION U(LDU,*),V(LDV,*) TR1 = 0.0 DO 20 I = 1, N S = 0.0 DO 10 K = 1, M S = S+U(I,K)*V(K,I) 10 CONTINUE TR1 = TR1+S 20 CONTINUE TR2 = 0.0 IF (IND.NE.2) RETURN DO 50 I = 1, N DO 40 J = 1, N S1 = 0.0 S2 = 0.0 DO 30 K = 1, M S1 = S1+U(I,K)*V(K,J) S2 = S2+U(J,K)*V(K,I) 30 CONTINUE TR2 = TR2+S1*S2 40 CONTINUE 50 CONTINUE RETURN C END C *ELLPTS SUBROUTINE ELLPTS (E,N,X,Y,IFLAG) C C----------------------------------------------------------------------- C ELLPTS WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899. C C FOR: COMPUTING N REGULARLY-SPACED POINTS ON THE PERIMETER OF THE C ELLIPSE DEFINED BY C C [E(1) E(3)] [X-E(4)] C [X-E(4) Y-E(5)] [ ] [ ] = E(6) . C [E(3) E(2)] [Y-E(5)] C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED JANUARY 5, 1987 C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * E = VECTOR (LENGTH 6) OF CONSTANTS DEFINING ELLIPSE (REAL) C C * N = NUMBERS OF POINTS TO BE GENERATED ON PERIMETER OF ELLIPSE C (INTEGER) C C X = VECTOR (LENGTH N) OF X COORDINATES OF POINTS ON PERIMETER C OF ELLIPSE (REAL) C C Y = VECTOR (LENGTH N) OF Y COORDINATES OF POINTS ON PERIMETER C OF ELLIPSE (REAL) C C IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION: C 0 -> NO ERRORS DETECTED C 1 -> N<1 C 2 -> E(1) OR E(2) OR E(6) IS NON-POSITIVE C 3 -> E(1)*E(2)<=E(3)**2 C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION E(*),X(*),Y(*) IF (N.LE.0) THEN IFLAG = 1 ELSEIF (AMIN1(E(1),E(2),E(6)).LE.0.0) THEN IFLAG = 2 ELSEIF (E(1)*E(2).LE.E(3)**2) THEN IFLAG = 3 ELSE U = 2.0*E(3) V = E(1)-E(2) IF (U.EQ.0.0.AND.V.EQ.0.0) THEN THETA = 0.0 ELSE THETA = 0.5*ATAN2(U,V) ENDIF CT = COS(THETA) ST = SIN(THETA) Q = SQRT(V**2+U**2) U = 2.0*E(6) V = E(1)+E(2) A = SQRT(U/(V+Q)) B = SQRT(U/(V-Q)) C = 8.0*ATAN(1.0)/REAL(N) DO 10 I = 1, N PHI = C*REAL(I) CP = COS(PHI) SP = SIN(PHI) X(I) = A*CT*CP-B*ST*SP+E(4) Y(I) = A*ST*CP+B*CT*SP+E(5) 10 CONTINUE IFLAG = 0 ENDIF RETURN END *FACTOR SUBROUTINE FACTOR (N,NFAC,L,IFLAG) C C----------------------------------------------------------------------- C FACTOR WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: FINDING THE PRIME FACTORS OF THE ABSOLUTE VALUE OF N, AN C INTEGER OF 10 DIGITS OR LESS. THE FACTORS ARE RETURNED IN C THE VECTOR NFAC WITH L INDICATING THE NUMBER OF FACTORS. C C ********************************************************** C *** IN THE CALLING PROGRAM THE VECTOR NFAC SHOULD BE *** C *** DIMENSIONED AT LEAST 33 TO ASSURE THAT IT WILL NOT *** C *** BE OVER-INDEXED IN THIS ROUTINE (2**34 > 10**10). *** C ********************************************************** C C UPON THE FIRST CALL TO THIS SUBROUTINE, ALL PRIME INTEGERS IN C THE RANGE [2,100003], IN ASCENDING ORDER, ARE READ FROM A FILE C NAMED 'PRIMES' INTO THE VECTOR IPRM. N IS THEN 'FACTORED BY C DIVISION' USING THESE PRIMES. THE ALGORITHM USED IS A SLIGHT C MODIFICATION OF ALGORITHM 'A' IN THE REFERENCE. THE RETURNED C VALUES FROM THIS ROUTINE ARE VALID ONLY IF THE ERROR FLAG C (IFLAG) IS ZERO ON RETURN. C C SUBPROGRAMS CALLED: -NONE- C C DATA FILES CALLED: PRIMES C C CURRENT VERSION COMPLETED FEBRUARY 4, 1987 C C REFERENCE: KNUTH, DONALD E., 'THE ART OF COMPUTER PROGRAMMING', C VOL. 2/SEMINUMERICAL ALGORITHMS, 2ND EDITION, ADDISON- C WESLEY PUBLISHING CO., 1981, PP. 364-398. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * N = AN INTEGER OF 10 DIGITS OR LESS. THE ABSOLUTE VALUE OF N C WILL BE FACTORED INTO PRIMES. C C NFAC = A VECTOR (LENGTH L) OF PRIME FACTORS OF ABS(N) (INTEGER) C C L = THE NUMBER OF PRIME FACTORS OF N (INTEGER) C C IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION: C 0 -> NO ERRORS DETECTED C 1 -> N HAS MORE THAN 10 DIGITS C 2 -> AN END-OF-FILE ENCOUNTERED IN READING FILE 'PRIMES' C 3 -> A FORMAT ERROR ENCOUNTERED IN READING FILE 'PRIMES' C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C PARAMETER (JX=9593) DIMENSION IPRM(JX),NFAC(*) DATA NCALL / 0 / M = IABS(N) IFLAG = 0 L = 0 IF (M.EQ.0) RETURN IF (M.GT.9999999999) THEN IFLAG = 1 ELSEIF (M.EQ.1) THEN L = 1 NFAC(1) = 1 ELSE IF (NCALL.EQ.0) THEN C C--- READ ALL PRIMES FROM 2 TO 100,003 C OPEN (UNIT=4,FILE='PRIMES') REWIND 4 READ (4,*,ERR=50,END=40) (IPRM(J),J=1,JX) ENDIF C C--- COUNT NUMBER OF CALLS TO THIS SUBROUTINE C NCALL = NCALL+1 C C--- IMPLEMENT ALGORITHM A (FACTORING BY PRIMES) ON PP. 364-5 OF THE C--- REFERENCE C K = 1 10 IF (M.EQ.1) RETURN 20 IQUO = M/IPRM(K) IREM = MOD(M,IPRM(K)) IF (IREM.NE.0) GO TO 30 L = L+1 NFAC(L) = IPRM(K) M = IQUO GO TO 10 30 IF (IQUO.GT.IPRM(K)) THEN K = K+1 GO TO 20 ENDIF L = L+1 NFAC(L) = M ENDIF RETURN C C--- SET ERROR FLAG FOR 'END-OF-FILE' ON UNIT 4 C 40 IFLAG = 2 RETURN C C--- SET ERROR FLAG FOR 'READ ERROR' ON UNIT 4 C 50 IFLAG = 3 RETURN END *FIBMIN SUBROUTINE FIBMIN (FUNC,XLO,XHI,TOL,XMIN,FXMIN) C----------------------------------------------------------------------- C FIBMIN WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING THE X VALUE (XMIN) AT WHICH THE FUNCTION FUNC(X) IS C MINIMAL IN THE INTERVAL [XLO,XHI]. IN THAT INTERVAL FUNC IS C ASSUMED TO BE UNIMODAL. THE TOLERANCE ON XMIN (TOL) MUST BE C SPECIFIED. THE FUNCTION VALUE AT XMIN (FXMIN) IS RETURNED. C THE MINIMUM IS FOUND USING A FIBONACCI SEARCH ALGORITHM. C C NOTE: IF TOL < [XHI-XLO]/[THE JX(TH) FIBONACCI NUMBER], THEN JX C SHOULD BE INCREASED. THE 100(TH) FIBONACCI NUMBER IS ABOUT C 3.5E+20. C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED AUGUST 12, 1988 C----------------------------------------------------------------------- PARAMETER (JX=100) DIMENSION IFIB(JX) DATA IFIB(1),IFIB(2) / 1,1 / XRANGE = XHI-XLO BOUND = XRANGE/TOL J = 2 10 J = J+1 IF (J.GT.JX) STOP ' J>JX IN SUBROUTINE FIBMIN - STOP' IFIB(J) = IFIB(J-2)+IFIB(J-1) IF (REAL(IFIB(J)).LT.BOUND) GO TO 10 DELTA = XRANGE/REAL(IFIB(J)) C1 = XLO C2 = DELTA K = 0 M = IFIB(J-2) N = IFIB(J-1) L = IFIB(J) F1 = FUNC(XLO) F2 = FUNC(XHI) I = J-2 20 I = I-1 IF (I.LT.1) GO TO 30 F1 = FUNC(C1+C2*REAL(M)) F2 = FUNC(C1+C2*REAL(N)) IF (F1.EQ.F2) THEN K = M L = N I = I-2 IF (I.LT.1) GO TO 30 M = K+IFIB(I) N = K+L-M F3 = FUNC(C1+C2*REAL(M)) IF (F3.GE.F1) GO TO 30 F4 = FUNC(C1+C2*REAL(N)) IF (F4.GE.F1) GO TO 30 ELSEIF (F1.GT.F2) THEN K = M M = N N = K+L-M ELSE L = N N = M M = K+L-N ENDIF GO TO 20 30 FBMIN = 0.5*REAL(K+L) XMIN = C1+C2*FBMIN FXMIN = FUNC(XMIN) RETURN C END *IGCD INTEGER FUNCTION IGCD (M,N) C C----------------------------------------------------------------------- C IGCD WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING THE GREATEST COMMON DIVISOR OF THE INTEGERS M AND N. C IF M=N=0 THEN THE G.C.D. IS DEFINED TO BE ZERO. C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED FEBRUARY 3, 1987 C----------------------------------------------------------------------- C I = IABS(N) IGCD = IABS(M) IF (IGCD.NE.0) GO TO 10 IGCD = I RETURN 10 K = MOD(I,IGCD) IF (K.EQ.0) RETURN I = IGCD IGCD = K GO TO 10 END *IGCDM INTEGER FUNCTION IGCDM (M,LDM,NR,NC) C C----------------------------------------------------------------------- C IGCDM WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING THE GREATEST COMMON DIVISOR OF THE NR*NC INTEGERS IN C THE MATRIX M. IF M(1,1)=M(1,2)=...=M(NR,NC)=0 THEN THE G.C.D. C IS DEFINED TO BE ZERO (NR IS THE NUMBER OF ROWS OF M, NC IS C THE NUMBER OF COLUMNS OF M, AND LDM IS THE LEADING DIMENSION C OF M AS DIMENSIONED IN THE CALLING PROGRAM). NO CHECK IS MADE C FOR THE ERROR CONDITION NR>LDM. C C SUBPROGRAMS CALLED: IGCD C C CURRENT VERSION COMPLETED FEBRUARY 3, 1987 C----------------------------------------------------------------------- C DIMENSION M(LDM,*) IGCDM = M(1,1) DO 20 J = 1, NC DO 10 I = 1, NR L = IABS(M(I,J)) IGCDM = IGCD(IGCDM,L) IF (IGCDM.EQ.1) RETURN 10 CONTINUE 20 CONTINUE RETURN END *IGCDV INTEGER FUNCTION IGCDV (M,N) C C----------------------------------------------------------------------- C IGCDV WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING THE GREATEST COMMON DIVISOR OF THE N INTEGERS IN THE C VECTOR M. IF M(1)=M(2)=...=M(N)=0 THEN THE G.C.D. IS DEFINED C TO BE ZERO. C C SUBPROGRAMS CALLED: IGCD C C CURRENT VERSION COMPLETED FEBRUARY 3, 1987 C----------------------------------------------------------------------- C DIMENSION M(*) IGCDV = M(1) DO 10 I = 1, N L = IABS(M(I)) IGCDV = IGCD(IGCDV,L) IF (IGCDV.EQ.1) RETURN 10 CONTINUE RETURN END *LINSYS SUBROUTINE LINSYS (A,LDA,N,B,IPVT,IFLAG) C C----------------------------------------------------------------------- C LINSYS WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: SOLVING THE NXN LINEAR SYSTEM OF EQUATIONS AX=B FOR THE VECTOR C X. GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING IS USED TO C FACTOR A INTO LU. THE TRIANGULAR SYSTEM LY=B IS THEN SOLVED C FOR Y, AND SIMILARLY UX=Y IS SOLVED FOR X. THIS ALGORITHM IS C DESCRIBED ON PAGES 1.11 AND 1.14 OF THE REFERENCE BELOW. C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED JANUARY 6, 1987 C C REFERENCE: DONGARRA, J.J., MOLER, C.B., BUNCH, J.R., AND STEWART, C G.W., "LINPACK USERS' GUIDE", SIAM, PHILADELPHIA, 1979. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * A = MATRIX (SIZE NXN) OF COEFFICIENTS ON INPUT AND LU C FACTORIZATION ON OUTPUT (REAL) C C * LDA = THE LEADING DIMENSION OF MATRIX A (>=N) (INTEGER) C C * N = NUMBER OF EQUATIONS AND UNKNOWNS (INTEGER) C C * B = VECTOR (LENGTH N) CONTAINING THE RIGHT HAND SIDE OF THE C SYSTEM ON INPUT AND THE SOLUTION VECTOR X ON OUTPUT (REAL) C C IPVT = VECTOR (LENGTH N) FOR PIVOTING INFORMATION (INTEGER) C C IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION: C 0 -> NO ERRORS DETECTED C K -> THE PIVOT ELEMENT AT THE KTH STEP OF THE ELIMINATION C (1<=K<=N) IS ZERO. THE LU FACTORIZATION OF A IS C VALID, BUT THE MATRIX U IS SINGULAR AND CANNOT BE USED C FOR SOLVING THE LINEAR SYSTEM. C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION A(LDA,*),B(*),IPVT(*) IFLAG = 0 C C--- COMPUTE LU FACTORIZATION OF A BY GAUSSIAN ELIMINATION WITH C--- PARTIAL PIVOTING. C DO 50 K = 1, N-1 AMAX = 0.0 DO 10 I = K, N IF (ABS(A(I,K)).GE.AMAX) THEN AMAX = ABS(A(I,K)) L = I ENDIF 10 CONTINUE IPVT(K) = L IF (A(L,K).EQ.0.0) THEN IFLAG = K GO TO 50 ENDIF IF (L.NE.K) THEN Q = A(K,K) A(K,K) = A(L,K) A(L,K) = Q ENDIF DO 20 I = K+1, N A(I,K) = -A(I,K)/A(K,K) 20 CONTINUE DO 40 J = K+1, N IF (L.NE.K) THEN Q = A(K,J) A(K,J) = A(L,J) A(L,J) = Q ENDIF DO 30 I = K+1, N A(I,J) = A(I,J)+A(I,K)*A(K,J) 30 CONTINUE 40 CONTINUE 50 CONTINUE IF (A(N,N).EQ.0.0) IFLAG = N IF (IFLAG.NE.0) RETURN C C--- SOLVE LY=B FOR Y BY BACK SUBSTITUTION C DO 70 K = 1, N-1 L = IPVT(K) IF (L.NE.K) THEN Q = B(K) B(K) = B(L) B(L) = Q ENDIF DO 60 I = K+1, N B(I) = B(I)+A(I,K)*B(K) 60 CONTINUE 70 CONTINUE C C--- SOLVE UX=Y FOR X BY BACK SUBSTITUTION C DO 90 K = N, 1, -1 B(K) = B(K)/A(K,K) DO 80 I = 1, K-1 B(I) = B(I)-A(I,K)*B(K) 80 CONTINUE 90 CONTINUE RETURN END *LISYPD SUBROUTINE LISYPD (A,LDA,N,B,IFLAG) C C----------------------------------------------------------------------- C LISYPD WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: SOLVING THE NXN LINEAR SYSTEM OF EQUATIONS AX=B FOR THE VECTOR C X WHERE A IS SYMMETRIC AND POSITIVE DEFINITE. THE CHOLESKY C FACTORIZATION OF A YIELDS A LOWER TRIANGULAR MATRIX L SUCH C THAT A=LL'. THE TRIANGULAR SYSTEM LY=B IS THEN SOLVED FOR C Y, AND SIMILARLY L'X=Y IS SOLVED FOR X. THE SOLUTION VECTOR C X IS RETURNED IN B. THIS ALGORITHM IS DESCRIBED ON PAGES C 3.10 AND 3.11 OF THE REFERENCE BELOW. NOTE THAT ONLY THE C LOWER TRIANGLE OF A IS USED. C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED NOVEMBER 29, 1988 C C REFERENCE: DONGARRA, J.J., MOLER, C.B., BUNCH, J.R., AND STEWART, C G.W., "LINPACK USERS' GUIDE", SIAM, PHILADELPHIA, 1979. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * A = MATRIX (SIZE NXN) OF COEFFICIENTS ON INPUT AND LU C FACTORIZATION ON OUTPUT (REAL) C C * LDA = THE LEADING DIMENSION OF MATRIX A (>=N) (INTEGER) C C * N = NUMBER OF EQUATIONS AND UNKNOWNS (INTEGER) C C * B = VECTOR (LENGTH N) CONTAINING THE RIGHT HAND SIDE OF THE C SYSTEM ON INPUT AND THE SOLUTION VECTOR X ON OUTPUT (REAL) C C IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION: C -3 -> N<1 C 0 -> NO ERRORS DETECTED C K -> THE KTH DIAGONAL ELEMENT OF THE CHOLESKY FACTOR- C IZATION IS <0, THUS THE FACTORIZATION CANNOT BE C COMPLETED. INDICATES A IS NOT POSITIVE DEFINITE. C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION A(LDA,*),B(*) IFLAG = 0 IF (N.LT.1) THEN IFLAG = -3 RETURN C ENDIF C C--- COMPUTE THE CHOLESKY FACTORIZATION A=LL' C DO 40 K = 1, N DO 20 I = 1, K-1 S = A(K,I) DO 10 J = 1, I-1 S = S-A(I,J)*A(K,J) 10 CONTINUE A(K,I) = S/A(I,I) 20 CONTINUE S = A(K,K) DO 30 J = 1, K-1 S = S-A(K,J)**2 30 CONTINUE IF (S.LE.0.0) THEN IFLAG = K RETURN C ENDIF A(K,K) = SQRT(S) 40 CONTINUE C C--- SOLVE LY=B (STORE Y IN B) C DO 60 I = 1, N S = B(I) DO 50 J = 1, I-1 S = S-A(I,J)*B(J) 50 CONTINUE B(I) = S/A(I,I) 60 CONTINUE C C--- SOLVE L'X=Y (STORE X IN B) C DO 80 I = N, 1, -1 S = B(I) DO 70 J = I+1, N S = S-A(J,I)*B(J) 70 CONTINUE B(I) = S/A(I,I) 80 CONTINUE RETURN C END *LSQSVD SUBROUTINE LSQSVD (X,Y,B,WORK,S,E,V,N,M,LDX,MX,K,IFLAG) C C----------------------------------------------------------------------- C LSQSVD WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: SOLVING THE LINEAR SYSTEM Y = XB FOR B GIVEN Y AND X WHERE C Y IS OF LENGTH N AND X IS AN N BY M MATRIX (N>=M). IF X IS C OF FULL RANK (M) THEN B IS THE (UNWEIGHTED) LEAST SQUARES C SOLUTION. IF X IS NOT OF FULL RANK THEN B IS THE MINIMUM C NORM LEAST SQUARES SOLUTION. A LINPACK ROUTINE IS USED TO C PERFORM A SINGULAR VALUE FACTORIZATION OF X. C C NOTE: RNDERR IS A MACHINE DEPENDENT CONSTANT WHICH IS THE MACHINE C ROUNDING ERROR (OR A LITTLE LARGER). IT IS USED TO DETERMINE C WHEN A SINGULAR VALUE IS ZERO. SEE THE DISCUSSION OF THIS C POINT ON PAGE 11.2 OF REFERENCE 1. C C SUBPROGRAMS CALLED: SSVDC (LINPACK) C C CURRENT VERSION COMPLETED SEPTEMBER 28, 1987 C C REFERENCES: C C 1) DONGARRA, J.J., MOLER, C.B., BUNCH, J.R., AND STEWART, G.W., C "LINPACK USERS' GUIDE", SIAM, PHILADELPHIA, 1979, CH. 11. C C 2) LAWSON, CHARLES L. AND HANSON, RICHARD J., "SOLVING LEAST C SQUARES PROBLEMS", PRENTICE-HALL, INC., CH. 7. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * X(LDX,*) = MATRIX (SIZE N BY M). THE SECOND DIMENSION OF X, C DEFINED IN THE CALLING PROGRAM, MUST BE >=M [REAL] C C * Y(*) = VECTOR (LENGTH N) [REAL] C C B(*) = VECTOR (LENGTH M) CONTAINING THE LEAST SQUARES C SOLUTION IF RANK(X)=M, OR THE MINIMUM NORM LEAST C SQUARES SOLUTION IF RANK(X)=N) [INTEGER] C C * MX = LEADING DIMENSION OF MATRIX V (MX>=M) [INTEGER] C C K = NUMBER OF NONZERO SINGULAR VALUES ON RETURN [INTEGER] C C IFLAG = ERROR INDICATOR ON OUTPUT [INTEGER] INTERPRETATION: C 1 -> N>LDX OR M>MX. C 2 -> N INFO<>0 RETURNED FROM SSVDC (SINGULAR VALUES WERE C NOT COMPUTED CORRECTLY, THUS MOORE-PENROSE C PSEUDO-INVERSE NOT COMPUTED). C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION X(LDX,*),Y(*),B(*),WORK(*),S(*),E(*),V(MX,*) DATA RNDERR / 1.0E-14 / IFLAG = 0 IF (N.GT.LDX.OR.M.GT.MX) THEN IFLAG = 1 RETURN C ENDIF IF (N.LT.M) THEN IFLAG = 2 RETURN C ENDIF C C--- COMPUTE LARGEST ELEMENT OF X (IN ABSOLUTE VALUE) C XMAX = 0.0 DO 20 I = 1, N DO 10 J = 1, M XMAX = AMAX1(XMAX,ABS(X(I,J))) 10 CONTINUE 20 CONTINUE C C--- COMPUTE CUTOFF POINT FOR A SINGULAR VALUE BEING ZERO C CUTOFF = 10.0*RNDERR*XMAX C C--- PERFORM SINGULAR VALUE FACTORIZATION USING LINPACK C CALL SSVDC (X,LDX,N,M,S,E,X,LDX,V,MX,WORK,21,INFO) C C--- CHECK WHETHER SINGULAR VALUES HAVE BEEN COMPUTED CORRECTLY C IF (INFO.NE.0) THEN IFLAG = 3 RETURN C ENDIF C C--- DETERMINE NUMBER OF NONZERO SINGULAR VALUES C K = 0 DO 30 J = 1, M IF (ABS(S(J)).GT.CUTOFF) THEN K = K+1 ELSE GO TO 80 C ENDIF 30 CONTINUE C C--- COMPUTE WORK = INV(S)U'Y C 80 DO 50 I = 1, K T = 0.0 DO 40 J = 1, N T = T+X(J,I)*Y(J) 40 CONTINUE WORK(I) = T/S(I) 50 CONTINUE C C--- COMPUTE B = VW C DO 70 I = 1, M T = 0.0 DO 60 J = 1, K T = T+V(I,J)*WORK(J) 60 CONTINUE B(I) = T 70 CONTINUE END *LSTSQR SUBROUTINE LSTSQR (X,LDX,N,M,Y,RES,SRES,HAT,UNK,SDUNK,LU,IND,ID) C C----------------------------------------------------------------------- C LSTSQR WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: ANALYZING THE STATISTICAL MODEL Y = X*B + ERROR BY UNWEIGHTED C LINEAR LEAST SQUARES. THREE SUBROUTINES IN LINPACK ARE CALLED. C COMPUTED QUANTITIES ARE: C C 1) ESTIMATES OF UNKNOWNS (B) C 2) STANDARD ERRORS OF ESTIMATES OF UNKNOWNS C 3) T-VALUES FOR UNKNOWNS C 4) RESIDUAL STANDARD DEVIATION C 5) PREDICTED VALUES OF OBSERVATIONS C 6) RESIDUALS C 7) STANDARDIZED RESIDUALS C 8) DIAGONAL ELEMENTS OF 'HAT MATRIX' C 9) DURBIN-WATSON STATISTIC FOR SEQUENTIAL CORRELATION C OF RESIDUALS C C THESE QUANTITIES ARE AUTOMATICALLY PRINTED AT USER'S OPTION. C ALL BUT 3 AND 5 ARE RETURNED IN THE PASSED ARGUMENTS. AN C OPTION EXISTS TO SUPPRESS THE COMPUTATION OF 7 AND 8 SINCE C THE TIME REQUIREMENTS FOR THAT ARE ON THE ORDER OF N**2. TO C ILLUSTRATE THE TIME REQUIREMENTS, A FIFTH-DEGREE POLYNOMIAL C WAS FIT TO N POINTS WITH AND WITHOUT COMPUTING 7 AND 8. THE C CPU TIMES WERE AS FOLLOWS ON THE CYBER 180/855 AT N.I.S.T.: C C WITHOUT 7,8 WITH 7,8 C N (IND=0) (IND=10) C ---- ----------- -------- C 20 0.005 0.017 C 100 0.013 0.153 C 500 0.050 3.045 C 2500 0.277 83.763 C C THE RESIDUAL STANDARD DEVIATION IS RETURNED IN X(1,1) AND THE C DURBIN-WATSON STATISTIC IS RETURNED IN X(2,1). THE VECTOR OF C OBSERVATIONS, Y, IS NOT CHANGED BY THIS SUBROUTINE. C C NOTE: THE VECTORS WHICH APPEAR IN THE ARGUMENT LIST ARE NECESSARY C FOR WORK SPACE IN THIS SUBROUTINE. THE USER SHOULD NOT C ATTEMPT TO REMOVE ANY OF THEM. C C NOTE: IN THE CALLING PROGRAM VECTORS Y, RES, SRES, AND HAT MUST C BE DIMENSIONED AT LEAST N. VECTORS UNK AND SDUNK MUST BE C DIMENSIONED AT LEAST M. A SAMPLE MAIN PROGRAM MIGHT BE: C C PARAMETER (NX=200,MX=10) C DIMENSION X(NX,MX) C DIMENSION Y(NX),RES(NX),SRES(NX),HAT(NX) C DIMENSION UNK(MX),SDUNK(MX) C CHARACTER*60 ID C ID = '' C LU = 6 C N = C M = C . C . C . C C CALL LSTSQR(X,NX,N,M,Y,RES,SRES,HAT,UNK,SDUNK,LU,IND,ID) C . C . C . C END C C SUBPROGRAMS CALLED: SQRDC, SQRSL, STRSL (LINPACK) C C CURRENT VERSION COMPLETED APRIL 21, 1989 C C REFERENCE: DONGARRA, J.J., MOLER, C.B., BUNCH, J.R., AND STEWART, C G.W., "LINPACK USERS' GUIDE", SIAM, PHILADELPHIA, 1979, C CH. 9. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * X(LDX,*) = MATRIX (SIZE N BY M) IN THE MODEL Y = XB + ERROR. C THE SECOND DIMENSION OF X, DEFINED IN THE CALLING C PROGRAM, MUST BE >=M. THIS MATRIX IS CHANGED C CONSIDERABLY ON OUTPUT. [REAL] C C * LDX = LEADING DIMENSION OF MATRIX X (LDX>=N) [INTEGER] C C * N = NUMBER OF ROWS OF MATRIX X (N<=LDX) [INTEGER] C C * M = NUMBER OF COLUMNS OF MATRIX X (M<=N) [INTEGER] C C * Y(*) = VECTOR (LENGTH N) OF OBSERVATIONS. NOT CHANGED ON C OUTPUT. [REAL] C C RES(*) = VECTOR (LENGTH N) USED FOR WORKSPACE. ON OUTPUT C CONTAINS RESIDUALS [REAL] C C SRES(*) = VECTOR (LENGTH N) USED FOR WORKSPACE. ON OUTPUT C CONTAINS STANDARDIZED RESIDUALS [REAL] C C HAT(*) = VECTOR (LENGTH N) USED FOR WORKSPACE. ON OUTPUT C CONTAINS DIAGONAL ELEMENTS OF THE HAT MATRIX [REAL] C C UNK(*) = VECTOR (LENGTH M) CONTAINING ESTIMATES OF UNKNOWNS C ON OUTPUT [REAL] C C SDUNK(*) = VECTOR (LENGTH M) USED FOR WORKSPACE. ON OUTPUT C CONTAINS STANDARD ERRORS OF ESTIMATES OF UNKNOWNS C [REAL] C C * LU = LOGICAL UNIT OF PRINT FILE WHERE AUTOMATIC PRINTOUT C OCCURS, IF REQUESTED [INTEGER] C C * IND = INTEGER INDICATOR VARIABLE OF THE FORM AB WHERE C A=0 TO SUPPRESS COMPUTING OF HAT MATRIX DIAGONALS C A=1 TO COMPUTE HAT MATRIX DIAGONALS C B=0 FOR NO PRINTOUT TO LOGICAL UNIT LU C B=1 FOR MINIMAL PRINTOUT TO LOGICAL UNIT LU (NO C RESIDUALS, STANDARDIZED RESIDUALS, OR HAT C MATRIX DIAGONAL ELEMENTS PRINTED) C B=2 FOR COMPLETE PRINTOUT TO LOGICAL UNIT LU C C * ID = CHARACTER VARIABLE OF LENGTH 60 FOR PASSING THE C IDENTIFICATION OF THE DATA SET BEING FITTED. C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION X(LDX,*),Y(*),RES(*),SRES(*),HAT(*),UNK(*),SDUNK(*) CHARACTER*1 STAR CHARACTER*60 ID C C--- DETERMINE INDICATORS FOR PRINTOUT AND COMPUTATION OF HAT C--- MATRIX DIAGONAL C IHAT = IND/10 IPRT = IND-10*IHAT C C--- CHECK FOR ILLEGAL INPUT ARGUMENTS C IF (M.LT.1) THEN IFLAG = 1 ELSEIF (N.LE.M) THEN IFLAG = 2 ELSEIF (LDX.LT.N) THEN IFLAG = 3 ELSEIF (IHAT.LT.0.OR.IHAT.GT.1.OR.IPRT.LT.0.OR.IPRT.GT.2) THEN IFLAG = 4 ELSE IF (IPRT.GE.1) WRITE (LU,1000) ID C C--- COMPUTE QR FACTORIZATION OF X C CALL SQRDC (X,LDX,N,M,SDUNK,JDUM,DUM,0) C C--- COMPUTE DIAGONAL OF HAT MATRIX, X*INV(X'X)*X'. IF COMPUTATION OF C--- THIS QUANTITY IS SUPPRESSED, SET EACH ELEMENT TO M/N. C RMN = REAL(M)/REAL(N) HATLIM = 2.0*RMN IF (IHAT.EQ.1) THEN DO 10 I = 1, N SRES(I) = 0.0 10 CONTINUE DO 30 I = 1, N SRES(I) = 1.0 IF (I.GT.1) SRES(I-1) = 0.0 CALL SQRSL (X,LDX,N,M,SDUNK,SRES,DUM,RES,DUM,DUM,DUM, * 01000,INFO) S = 0.0 DO 20 J = 1, M S = S+RES(J)**2 20 CONTINUE HAT(I) = S 30 CONTINUE ELSE DO 40 I = 1, N HAT(I) = RMN 40 CONTINUE ENDIF C C--- COMPUTE ESTIMATES OF UNKNOWNS AND RESIDUALS C CALL SQRSL (X,LDX,N,M,SDUNK,Y,DUM,SRES,UNK,RES,DUM,00110,INFO) C C--- COMPUTE NUMERATOR OF DURBIN-WATSON STATISTIC C DW = 0.0 DO 50 I = 2, N DW = DW+(RES(I)-RES(I-1))**2 50 CONTINUE C C--- COMPUTE RESIDUAL SUM-OF-SQUARES C S = 0.0 DO 60 I = 1, N S = S+RES(I)**2 60 CONTINUE C C--- COMPUTE DURBIN-WATSON STATISTIC, DEGREES OF FREEDOM, AND C--- RESIDUAL STANDARD DEVIATION C DW = DW/S IDF = N-M RSD = SQRT(S/REAL(IDF)) IF (IPRT.GE.1) WRITE (LU,1030) N,M,IDF,RSD IF (IPRT.GE.1) WRITE (LU,1010) C C--- COMPUTE STANDARD ERRORS AND T-VALUES FOR THE UNKNOWNS C DO 90 J = 1, M DO 70 K = 1, M SRES(K) = 0.0 70 CONTINUE SRES(J) = 1.0 CALL STRSL (X,LDX,M,SRES,11,INFO) S = 0.0 DO 80 K = 1, M S = S+SRES(K)**2 80 CONTINUE SDUNK(J) = RSD*SQRT(S) T = UNK(J)/SDUNK(J) IF (IPRT.GE.1) WRITE (LU,1050) J,UNK(J),SDUNK(J),T 90 CONTINUE IF (IPRT.GE.1) WRITE (6,1080) DW IF (IPRT.GE.2) WRITE (6,1020) C C--- COMPUTE STANDARDIZED RESIDUALS AND INDICATE LARGE HAT MATRIX C--- DIAGONAL ELEMENTS, IF ANY C DO 100 I = 1, N IF (HAT(I).GT.HATLIM) THEN STAR = '*' ELSE STAR = ' ' ENDIF PRED = Y(I)-RES(I) IF (HAT(I).GE.1.0) THEN SRES(I) = 0.0 ELSE SRES(I) = RES(I)/(RSD*SQRT(1.0-HAT(I))) ENDIF IF (IPRT.GE.2) WRITE (LU,1040) I,Y(I),PRED,RES(I),SRES(I), * HAT(I),STAR 100 CONTINUE C C--- WRITE MESSAGE ABOUT HAT MATIX DIAGONAL ELEMENTS C IF (IPRT.GE.2) THEN IF (IHAT.EQ.1) THEN WRITE (LU,1060) HATLIM ELSE WRITE (LU,1070) ENDIF ENDIF C C--- RETURN RESIDUAL STANDARD DEVIATION AND DURBIN-WATSON STATISTIC C--- IN MATRIX X (WHICH HAS ALREADY BEEN CHANGED) C X(1,1) = RSD X(2,1) = DW ENDIF RETURN C C 1000 FORMAT ('1',18X,'UNWEIGHTED LINEAR LEAST SQUARES (LSTSQR)'///7X, * 'ID: ',A60) 1010 FORMAT (//15X,'UNKNOWN',8X,'STD. DEV.',8X,'T-VALUE'/) 1020 FORMAT (//50X,'STANDARDIZED',2X,'HAT MATRIX'/11X,'OBSERVED',6X, * 'PREDICTED',5X,'RESIDUAL',5X,'RESIDUAL',5X,'DIAGONAL'/) 1030 FORMAT (//10X,'NUMBER OF OBSERVATIONS = ',I5/14X,'NUMBER OF ', * 'UNKNOWNS = ',I5/8X,'RESIDUAL DEG. OF FREEDOM = ',I5,3X, * 'RESIDUAL STD. DEV. = ',G12.5) 1040 FORMAT (1X,I5,2X,2G15.8,G12.5,F9.3,3X,F9.4,1X,A1) 1050 FORMAT (5X,I5,1X,G16.8,1X,G13.5,2X,G14.3) 1060 FORMAT (//10X,'DIAGONAL ELEMENTS OF HAT MATRIX WHICH EXCEED ', * 'RULE-OF-THUMB'/10X,'LIMIT 2*(# UNKNOWNS)/(# EQUATIONS) = ', * F5.3,' ARE INDICATED BY *'//) 1070 FORMAT (//10X,'DIAGONAL ELEMENTS OF HAT MATRIX WERE NOT COMPU', * 'TED. THEY WERE'/10X,'SET TO THEIR AVERAGE VALUE, (# UNKNOWN', * 'S)/(# EQUATIONS). AS A'/10X,'RESULT, THE STANDARDIZED RESI', * 'DUALS ARE APPROXIMATE.'//) 1080 FORMAT (//10X,'DURBIN-WATSON STATISTIC = ',F6.3) C END *MATCNO SUBROUTINE MATCNO (A,LDA,N,CONDNO,IFLAG) C C----------------------------------------------------------------------- C MATCNO WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: COMPUTING THE CONDITION NUMBER OF THE N BY N MATRIX A. THE C CONDITION NUMBER IS KAPPA = NORM(A)*NORM[INV(A)]. HERE THE C INFINITY-NORM IS USED. THE CONDITION NUMBER MEASURES HOW C "ILL-CONDITIONED" THE MATRIX A IS, AND THUS HOW ACCURATELY THE C SYSTEM A*X = B MAY BE SOLVED. IF A IS SINGULAR THEN KAPPA IS C INFINITE AND AN ERROR FLAG IS RETURNED. SEE THE REFERENCE FOR C A GOOD DISCUSSION OF THE PROBLEM. C C SUBPROGRAMS CALLED: MATINV (STSPAC) C ANORMI (ATTACHED) C C CURRENT VERSION COMPLETED JUNE 28, 1989 C C REFERENCE: GOLUB, G.H. AND VAN LOAN, C.F., "MATRIX COMPUTATIONS", C JOHNS HOPKINS UNIVERSITY PRESS, BALTIMORE, 1985, P. 72. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * A(LDA,*) = MATRIX (SIZE N BY N) WHOSE CONDITION NUMBER IS TO BE C DETERMINED. THE SECOND DIMENSION OF A, DEFINED IN THE C CALLING PROGRAM, MUST BE >=N [REAL] C C * LDA = THE LEADING DIMENSION OF A (>=N) [INTEGER] C C * N = SIZE OF MATRIX A (N<=LDA) [INTEGER] C C CONDNO = THE CONDITION NUMBER OF A (KAPPA) ON RETURN WHEN C IFLAG=0 [REAL] C C IFLAG = ERROR INDICATOR ON OUTPUT. THE OUTPUT VALUE OF IFLAG C FROM SUBROUTINE MATINV IS PASSED DIRECTLY THROUGH THIS C SUBROUTINE. IF IFLAG=0, THEN ALL IS WELL. OTHERWISE, C SEE DOCUMENTATION OF MATINV FOR INTERPRETATION OF THE C THE ERROR FLAG IFLAG [INTEGER] C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION A(LDA,*) FAC1 = ANORMI(A,LDA,N) CALL MATINV (A,LDA,N,IFLAG) IF (IFLAG.NE.0) THEN CONDNO = 0.999999E29 ELSE FAC2 = ANORMI(A,LDA,N) CONDNO = FAC1*FAC2 ENDIF RETURN C END C C======================================================================= C FUNCTION ANORMI (A,LDA,N) C C--- COMPUTE THE INFINITY-NORM OF THE N BY N MATRIX A C DIMENSION A(LDA,*) ANORMI = 0.0 DO 20 I = 1, N T = 0.0 DO 10 J = 1, N T = T+ABS(A(I,J)) 10 CONTINUE ANORMI = AMAX1(ANORMI,T) 20 CONTINUE RETURN C END *MATHDI SUBROUTINE MATHDI (X,LDX,N,M,Y,QTY,HAT,QRAUX) C C----------------------------------------------------------------------- C MATHDI WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: COMPUTING THE DIAGONAL OF THE 'HAT MATRIX' CORRESPONDING TO C THE 'DESIGN MATRIX' X IN THE (UNWEIGHTED) LINEAR LEAST SQUARES C MODEL Y = XB + ERROR. THE HAT MATRIX IS DEFINED TO BE (SEE C REFERENCES 2 AND 3 BELOW) H = X*INV(X'X)*X'. WHEN X HAS N C ROWS AND M COLUMNS (M<=N) THEN H IS N BY N. THE DIAGONAL C ENTRIES OF H MEASURE THE 'LEVERAGE' OR 'INFLUENCE' OF THE C CORRESPONDING OBSERVATIONS. WHEN X IS OF FULL RANK, THE SUM OF C THE DIAGONAL ENTRIES OF H IS M, THEREFORE, THE AVERAGE VALUE C OF THE ENTRIES IS M/N. A RULE-OF-THUMB GIVEN IN THE REFERENCES C IS THAT, FOR A GOOD DESIGN MATRIX, NO DIAGONAL ELEMENT OF H C EXCEEDS 2*M/N. C C THE COMPUTATIONAL PROCEDURE BELOW INCORPORATES THE STANDARD C QR FACTORIZATION OF X AS IMPLEMENTED IN LINPACK. WHEN X IS C NOT OF FULL RANK, THE RANK DEFICIENCY IS EXPRESSED IN THE C FACTOR R. SINCE ONLY THE FACTOR Q IS USED IN COMPUTING THE C HAT MATRIX, THE PROGRAM WILL RETURN RESULTS WHEN X IS NOT OF C FULL RANK. THOSE RESULTS WILL BE ERRONEOUS, HOWEVER. FOR THE C SAKE OF SIMPLICITY, NO CHECK IS MADE ON THE RANK OF X. IT IS C THEREFORE UP TO THE USER TO KNOW THAT X IS OF FULL RANK BEFORE C ACCEPTING THE RESULTS OF THIS ROUTINE. C C SUBPROGRAMS CALLED: SQRDC, SQRSL (LINPACK) C C CURRENT VERSION COMPLETED JANUARY 25, 1989 C C REFERENCES: C C 1) DONGARRA, J.J., MOLER, C.B., BUNCH, J.R., AND STEWART, G.W., C "LINPACK USERS' GUIDE", SIAM, PHILADELPHIA, 1979, CH. 9. C C 2) HOAGLIN, D.C. AND WELSCH, R.E., "THE HAT MATRIX IN REGRESSION C AND ANOVA", THE AMERICAN STATISTICIAN, VOL. 32, NO. 1, FEBRUARY C 1978, PP. 17-22. C C 3) VELLEMAN, P.F. AND WELSCH, R.E., "EFFICIENT COMPUTING OF C REGRESSION DIAGNOSTICS", THE AMERICAN STATISTICIAN, VOL. 35, C NO. 4, NOVEMBER 1981, PP. 234-242. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * X(LDX,*) = MATRIX (SIZE N BY M) IN THE MODEL Y = XB + ERROR. C THE SECOND DIMENSION OF X, DEFINED IN THE CALLING C PROGRAM, MUST BE >=M [REAL] C C * LDX = LEADING DIMENSION OF MATRIX X (LDX>=N) [INTEGER] C C * N = NUMBER OF ROWS OF MATRIX X (N<=LDX) [INTEGER] C C * M = NUMBER OF COLUMNS OF MATRIX X (M<=N) [INTEGER] C C Y(*) = VECTOR (LENGTH N) USED FOR WORKSPACE [REAL] C C QTY(*) = VECTOR (LENGTH N) USED FOR WORKSPACE [REAL] C C HAT(*) = VECTOR (LENGTH N) CONTAINING DIAGONAL ENTRIES OF C THE 'HAT MATRIX' ON OUTPUT [REAL] C C QRAUX(*) = VECTOR (LENGTH M) USED FOR WORKSPACE [REAL] C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION X(LDX,*),Y(*),QTY(*),HAT(*),QRAUX(*) C C--- INITIALIZE UNIT DIRECTIONAL VECTOR C DO 10 I = 1, N Y(I) = 0.0 10 CONTINUE C C--- COMPUTE QR FACTORIZATION OF X (USING LINPACK) C CALL SQRDC (X,LDX,N,M,QRAUX,JDUM,DUM,0) C C--- COMPUTE HAT MATRIX DIAGONAL FOR DESIGN MATRIX X C DO 30 I = 1, N Y(I) = 1.0 IF (I.GT.1) Y(I-1) = 0.0 CALL SQRSL (X,LDX,N,M,QRAUX,Y,DUM,QTY,DUM,DUM,DUM,01000,INFO) S = 0.0 DO 20 J = 1, M S = S+QTY(J)**2 20 CONTINUE HAT(I) = S 30 CONTINUE RETURN C END *MATINV SUBROUTINE MATINV (A,LDA,N,IFLAG) C C----------------------------------------------------------------------- C MATINV WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING THE INVERSE OF A GENERAL N BY N MATRIX IN PLACE, C I.E., THE INVERSE OVERWRITES THE ORIGINAL MATRIX. THE STEPS C OF THE ALGORITHM ARE DESCRIBED BELOW AS THEY OCCUR. ROW C INTERCHANGES ARE DONE AS NEEDED IN ORDER TO INCREASE THE C ACCURACY OF THE INVERSE MATRIX. WITHOUT INTERCHANGES THIS C ALGORITHM WILL FAIL WHEN ANY OF THE LEADING PRINCIPAL C SUBMATRICES ARE SINGULAR OR WHEN THE MATRIX ITSELF IS C SINGULAR. WITH INTERCHANGES THIS ALGORITHM WILL FAIL ONLY C WHEN THE MATRIX ITSELF IS SINGULAR. THE LEADING PRINCIPAL C C [A B C] C SUBMATRICES OF THE MATRIX [D E F] ARE [A] AND [A B] . C [G H I] [D E] C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED JANUARY 15, 1987 C C REFERENCE: STEWART, G.W., 'INTRODUCTION TO MATRIX COMPUTATIONS', C ACADEMIC PRESS, INC., 1973 C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS C C * A = MATRIX (SIZE NXN) TO BE INVERTED (REAL) C C * LDA = LEADING DIMENSION OF MATRIX A [LDA>=N] (INTEGER) C C * N = NUMBER OF ROWS AND COLUMNS OF MATRIX A (INTEGER) C C IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION: C -2 -> TOO MANY ROW INTERCHANGES NEEDED - INCREASE MX C -1 -> N>LDA C 0 -> NO ERRORS DETECTED C K -> MATRIX A FOUND TO BE SINGULAR AT THE KTH STEP OF C THE CROUT REDUCTION (1<=K<=N) C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C PARAMETER (MX=100) DIMENSION A(LDA,*),IEX(MX,2) IFLAG = 0 C C--- CHECK CONSISTENCY OF PASSED PARAMETERS C IF (N.GT.LDA) THEN IFLAG = -1 RETURN ENDIF C C--- COMPUTE A = LU BY THE CROUT REDUCTION WHERE L IS LOWER TRIANGULAR C--- AND U IS UNIT UPPER TRIANGULAR (ALGORITHM 3.4, P. 138 OF THE C--- REFERENCE) C NEX = 0 DO 70 K = 1, N DO 20 I = K, N S = A(I,K) DO 10 L = 1, K-1 S = S-A(I,L)*A(L,K) 10 CONTINUE A(I,K) = S 20 CONTINUE C C--- INTERCHANGE ROWS IF NECESSARY C Q = 0.0 L = 0 DO 30 I = K, N R = ABS(A(I,K)) IF (R.GT.Q) THEN Q = R L = I ENDIF 30 CONTINUE IF (L.EQ.0) THEN IFLAG = K RETURN ENDIF IF (L.NE.K) THEN NEX = NEX+1 IF (NEX.GT.MX) THEN IFLAG = -2 RETURN ENDIF IEX(NEX,1) = K IEX(NEX,2) = L DO 40 J = 1, N Q = A(K,J) A(K,J) = A(L,J) A(L,J) = Q 40 CONTINUE ENDIF C C--- END ROW INTERCHANGE SECTION C DO 60 J = K+1, N S = A(K,J) DO 50 L = 1, K-1 S = S-A(K,L)*A(L,J) 50 CONTINUE A(K,J) = S/A(K,K) 60 CONTINUE 70 CONTINUE C C--- INVERT THE LOWER TRIANGLE L IN PLACE (SIMILAR TO ALGORITHM 1.5, C--- P. 110 OF THE REFERENCE) C DO 100 K = N, 1, -1 A(K,K) = 1.0/A(K,K) DO 90 I = K-1, 1, -1 S = 0.0 DO 80 J = I+1, K S = S+A(J,I)*A(K,J) 80 CONTINUE A(K,I) = -S/A(I,I) 90 CONTINUE 100 CONTINUE C C--- INVERT THE UPPER TRIANGLE U IN PLACE (ALGORITHM 1.5, P. 110 OF C--- THE REFERENCE) C DO 130 K = N, 1, -1 DO 120 I = K-1, 1, -1 S = A(I,K) DO 110 J = I+1, K-1 S = S+A(I,J)*A(J,K) 110 CONTINUE A(I,K) = -S 120 CONTINUE 130 CONTINUE C C--- COMPUTE INV(A) = INV(U)*INV(L) C DO 160 I = 1, N DO 150 J = 1, N IF (J.GT.I) THEN S = 0.0 L = J ELSE S = A(I,J) L = I+1 ENDIF DO 140 K = L, N S = S+A(I,K)*A(K,J) 140 CONTINUE A(I,J) = S 150 CONTINUE 160 CONTINUE C C--- INTERCHANGE COLUMNS OF INV(A) TO REVERSE EFFECT OF ROW C--- INTERCHANGES OF A C DO 180 I = NEX, 1, -1 K = IEX(I,1) L = IEX(I,2) DO 170 J = 1, N Q = A(J,K) A(J,K) = A(J,L) A(J,L) = Q 170 CONTINUE 180 CONTINUE RETURN END *MATIPD SUBROUTINE MATIPD (A,LDA,N,IFLAG) C C----------------------------------------------------------------------- C MATIPD WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING THE INVERSE OF A SYMMETRIC POSITIVE DEFINITE N BY N C MATRIX IN PLACE, I.E., THE INVERSE OVERWRITES THE ORIGINAL C MATRIX. THE STEPS OF THE ALGORITHM ARE DESCRIBED BELOW AS C THEY OCCUR. ONLY THE LOWER TRIANGLE (INCLUDING THE DIAGONAL) C OF THE INPUT MATRIX IS USED IN THE INVERSION. THE COMPLETE C INVERSE MATRIX IS RETURNED ON OUTPUT. IF THE INPUT MATRIX IS C NOT POSITIVE DEFINITE THE INVERSE CANNOT BE COMPUTED BY THE C CHOLESKY FACTORIZATION USED HEREIN, THUS AN ERROR FLAG IS SET. C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED JANUARY 15, 1987 C C REFERENCE: STEWART, G.W., 'INTRODUCTION TO MATRIX COMPUTATIONS', C ACADEMIC PRESS, INC., 1973 C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS C C * A = SYMMETRIC POSITIVE DEFINITE MATRIX (SIZE NXN) TO BE INVERTED C (REAL) C C * LDA = LEADING DIMENSION OF MATRIX A [LDA>=N] (INTEGER) C C * N = NUMBER OF ROWS AND COLUMNS OF MATRIX A (INTEGER) C C IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION: C -1 -> N>LDA C 0 -> NO ERRORS DETECTED C K -> MATRIX A FOUND NOT TO BE POSITIVE DEFINITE AT THE C KTH STEP OF THE CHOLESKY FACTORIZATION (1<=K<=N) C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION A(LDA,*) IFLAG = 0 C C--- CHECK CONSISTENCY OF PASSED PARAMETERS C IF (N.GT.LDA) THEN IFLAG = -1 RETURN ENDIF C C--- COMPUTE THE LOWER TRIANGULAR MATRIX L OF THE CHOLESKY C--- FACTORIZATION A=LL' (ALGORITHM 3.9, P. 142 OF THE REFERENCE) C DO 40 K = 1, N DO 20 I = 1, K-1 S = A(K,I) DO 10 J = 1, I-1 S = S-A(I,J)*A(K,J) 10 CONTINUE A(K,I) = S/A(I,I) 20 CONTINUE S = A(K,K) DO 30 J = 1, K-1 S = S-A(K,J)**2 30 CONTINUE IF (S.LE.0.0) THEN IFLAG = K RETURN ENDIF A(K,K) = SQRT(S) 40 CONTINUE C C--- INVERT THE LOWER TRIANGLE L IN PLACE (SIMILAR TO ALGORITHM 1.5, C--- P. 110 OF THE REFERENCE) C DO 70 K = N, 1, -1 A(K,K) = 1.0/A(K,K) DO 60 I = K-1, 1, -1 S = 0.0 DO 50 J = I+1, K S = S+A(J,I)*A(K,J) 50 CONTINUE A(K,I) = -S/A(I,I) 60 CONTINUE 70 CONTINUE C C--- COMPUTE INV(A)=INV(L)'*INV(L) C DO 100 I = 1, N DO 90 J = I, N S = 0.0 DO 80 K = J, N S = S+A(K,I)*A(K,J) 80 CONTINUE A(I,J) = S A(J,I) = S 90 CONTINUE 100 CONTINUE RETURN END *MATMPI SUBROUTINE MATMPI (X,WORK,S,E,V,N,M,NX,MX,K,IFLAG) C C----------------------------------------------------------------------- C MATMPI WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING THE MOORE-PENROSE PSEUDO-INVERSE OF AN NXM MATRIX C X WHERE N >= M. THE TRANSPOSE OF THE INVERSE IS RETURNED IN C THE ORIGINAL MATRIX. A LINPACK ROUTINE IS USED TO PERFORM C A SINGULAR VALUE DECOMPOSITION OF X FROM WHICH THE INVERSE C X+ IS COMPUTED. IF X IS OF FULL RANK THEN X+ = INV(X'X)*X'. C C NOTE: RNDERR IS A MACHINE DEPENDENT CONSTANT WHICH IS THE MACHINE C ROUNDING ERROR (OR A LITTLE LARGER). IT IS USED TO DETERMINE C WHEN A SINGULAR VALUE IS ZERO (IN CASES WHERE THE USER HAS C REQUESTED AUTOMATIC DETERMINATION OF RANK BY INPUTTING K=0). C SEE DISCUSSION OF THIS POINT ON PAGE 11.2 OF REFERENCE 1. C C SUBPROGRAMS CALLED: SSVDC (LINPACK) C C CURRENT VERSION COMPLETED DECEMBER 14, 1989 C C REFERENCES: C C 1) DONGARRA, J.J., MOLER, C.B., BUNCH, J.R., AND STEWART, G.W., C "LINPACK USERS' GUIDE", SIAM, PHILADELPHIA, 1979, CH. 11. C C 2) LAWSON, CHARLES L. AND HANSON, RICHARD J., "SOLVING LEAST C SQUARES PROBLEMS", PRENTICE-HALL, INC., CH. 7. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * X(NX,*) = MATRIX (SIZE N BY M) WHOSE PSEUDO-INVERSE IS TO BE C COMPUTED. THE SECOND DIMENSION OF X, DEFINED IN THE C CALLING PROGRAM, MUST BE >=M. [REAL] C C WORK(*) = VECTOR (LENGTH N) USED AS WORKSPACE [REAL] C C S(*) = VECTOR (LENGTH M) OF SINGULAR VALUES IN DESCENDING C ORDER ON RETURN, PROVIDED INFO=0 ON RETURN [REAL] C C E(*) = VECTOR (LENGTH M) OF ZEROS ON RETURN, PROVIDED INFO=0 C ON RETURN [REAL] C C V(MX,*) = MATRIX (SIZE M BY M) USED FOR INTERMEDIATE C COMPUTATIONS [REAL] C C * N = NUMBER OF ROWS IN MATRIX X [INTEGER] C C * M = NUMBER OF COLUMNS IN MATRIX X (M<=N) [INTEGER] C C * NX = LEADING DIMENSION OF MATRIX X (NX>=N) [INTEGER] C C * MX = LEADING DIMENSION OF MATRIX V (MX>=M) [INTEGER] C C * K = ON INPUT: K>0 INDICATES KNOWN RANK OF MATRIX X. C K=0 INDICATES RANK SHOULD BE AUTOMATICALLY C DETERMINED BY PROGRAM. C ON OUTPUT: = UNCHANGED IF K>0 ON INPUT. C = COMPUTED RANK OF X IF K=0 ON INPUT. C C IFLAG = ERROR INDICATOR ON OUTPUT [INTEGER] INTERPRETATION: C 1 -> N>NX OR M>MX. C 2 -> N INFO<>0 RETURNED FROM SSVDC (SINGULAR VALUES WERE C NOT COMPUTED CORRECTLY, THUS MOORE-PENROSE C PSEUDO-INVERSE NOT COMPUTED). C 4 -> K<0 ON INPUT. C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION X(NX,*),WORK(*),S(*),E(*),V(MX,*) DATA RNDERR / 1.0E-14 / IFLAG = 0 IF (N.GT.NX.OR.M.GT.MX) THEN IFLAG = 1 RETURN C ENDIF IF (N.LT.M) THEN IFLAG = 2 RETURN C ENDIF IF (K.LT.0) THEN IFLAG = 4 RETURN C ENDIF IF (K.EQ.0) THEN C C--- COMPUTE LARGEST ELEMENT OF X (IN ABSOLUTE VALUE) C XMAX = 0.0 DO 20 I = 1, N DO 10 J = 1, M XMAX = AMAX1(XMAX,ABS(X(I,J))) 10 CONTINUE 20 CONTINUE C C--- COMPUTE CUTOFF POINT FOR A SINGULAR VALUE BEING ZERO C CUTOFF = 10.0*RNDERR*XMAX ENDIF C C--- PERFORM SINGULAR VALUE FACTORIZATION USING LINPACK C CALL SSVDC (X,NX,N,M,S,E,X,NX,V,MX,WORK,21,INFO) C C--- CHECK WHETHER SINGULAR VALUES HAVE BEEN COMPUTED CORRECTLY C IF (INFO.NE.0) THEN IFLAG = 3 RETURN C ENDIF IF (K.EQ.0) THEN C C--- DETERMINE NUMBER OF NONZERO SINGULAR VALUES C K = 0 DO 30 J = 1, M IF (ABS(S(J)).GT.CUTOFF) THEN K = K+1 ELSE GO TO 40 C ENDIF 30 CONTINUE ENDIF C C--- COMPUTE THE MOORE-PENROSE PSEUDO-INVERSE OF X (TRANSPOSED) C 40 DO 60 J = 1, M DO 50 L = 1, K V(J,L) = V(J,L)/S(L) 50 CONTINUE 60 CONTINUE DO 100 I = 1, N DO 80 J = 1, M T = 0.0 DO 70 L = 1, K T = T+V(J,L)*X(I,L) 70 CONTINUE E(J) = T 80 CONTINUE DO 90 J = 1, M X(I,J) = E(J) 90 CONTINUE 100 CONTINUE RETURN C END *MATXXI SUBROUTINE MATXXI (X,LDX,N,M,DET,IFLAG) C C----------------------------------------------------------------------- C MATXXI WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING THE INVERSE OF X'X AND THE DETERMINANT OF X'X WHERE C X IS AN N BY M MATRIX WITH N >= M. THE INVERSE IS RETURNED C IN THE FIRST M ROWS AND M COLUMNS OF THE ORIGINAL MATRIX X. C ENTRIES BELOW ROW M ARE LIKELY TO BE CHANGED. NOTE THAT ROWS C N+1 AND N+2 OF MATRIX X ARE USED FOR SCRATCH AREA, THEREFORE C LDX >= N+2. THE STEPS OF THE ALGORITHM ARE DESCRIBED BELOW C AS THEY OCCUR. C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED APRIL 8, 1988 C C REFERENCE: STEWART, G.W., 'INTRODUCTION TO MATRIX COMPUTATIONS', C ACADEMIC PRESS, INC., 1973 C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * X = MATRIX (SIZE NXM) FOR WHICH INV(X'X) AND DET(X'X) ARE TO C BE COMPUTED (REAL) C C * LDX = LEADING DIMENSION OF X [LDX>=N+2] (INTEGER) C C * N = NUMBER OF ROWS IN MATRIX X (INTEGER) C C * M = NUMBER OF COLUMNS IN MATRIX X [M<=N] (INTEGER) C C DET = DETERMINANT OF (X'X) (REAL) C C IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION: C -2 -> M>N C -1 -> LDX NO ERRORS DETECTED C K -> THE KTH DIAGONAL ELEMENT OF THE TRIANGULAR MATRIX R, C FROM THE QR FACTORIZATION, IS ZERO C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION X(LDX,*) IFLAG = 0 C C--- CHECK CONSISTENCY OF PASSED PARAMETERS C IF (LDX.LT.N+2) THEN IFLAG = -1 RETURN ENDIF IF (N.LT.M) THEN IFLAG = -2 RETURN ENDIF C C--- DEFINE SOME CONSTANTS C N1 = N+1 N2 = N+2 C C--- COMPUTE QR FACTORIZATION OF MATRIX X (ALGORITHM 3.8, P. 236 OF C--- THE REFERENCE) C L = MIN0(N-1,M) DO 70 K = 1, L E = 0.0 DO 10 I = K, N Q = ABS(X(I,K)) IF (Q.GT.E) E = Q 10 CONTINUE IF (E.EQ.0.0) THEN X(N+1,K) = 0.0 ELSE DO 20 I = K, N X(I,K) = X(I,K)/E 20 CONTINUE S = 0.0 DO 30 I = K, N S = S+X(I,K)**2 30 CONTINUE S = SIGN(SQRT(S),X(K,K)) X(K,K) = X(K,K)+S X(N1,K) = S*X(K,K) X(N2,K) = -E*S DO 60 J = K+1, M T = 0.0 DO 40 I = K, N T = T+X(I,K)*X(I,J) 40 CONTINUE T = T/X(N1,K) DO 50 I = K, N X(I,J) = X(I,J)-T*X(I,K) 50 CONTINUE 60 CONTINUE ENDIF 70 CONTINUE IF (N.EQ.M) X(N2,N) = X(N,N) C C--- MOVE DIAGONAL ELEMENTS BACK TO DIAGONAL AND COMPUTE DETERMINANT C DET = 1.0 DO 80 J = 1, M X(J,J) = X(N2,J) DET = DET*X(J,J)**2 IF (X(J,J).EQ.0.0) THEN IFLAG = J RETURN ENDIF 80 CONTINUE C C--- INVERT UPPER TRIANGULAR MATRIX R IN PLACE (ALGORITHM 1.5, P. 110 C--- OF THE REFERENCE) C DO 110 K = M, 1, -1 X(K,K) = 1.0/X(K,K) DO 100 I = K-1, 1, -1 S = 0.0 DO 90 J = I+1, K S = S-X(I,J)*X(J,K) 90 CONTINUE X(I,K) = S/X(I,I) 100 CONTINUE 110 CONTINUE C C--- COMPUTE INV(X'X) = INV(R)*INV(R)' C DO 140 I = 1, M DO 130 J = I, M S = 0.0 DO 120 K = J, M S = S+X(I,K)*X(J,K) 120 CONTINUE X(I,J) = S X(J,I) = S 130 CONTINUE 140 CONTINUE RETURN END *NEXPER SUBROUTINE NEXPER (RED,PINK,BROWN) C----------------------------------------------------------------------- C NEXPER COPIED BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 FROM THE REFERENCE BELOW (COLOR ADDED) C C FOR: COMPUTING THE NEXT PERMUTATION OF THE INTEGERS 1, 2, ..., N. C THE CALLING SEQUENCE IS C C CALL NEXPER (IPERM,N,LL) C C WHERE IPERM IS AN INTEGER VECTOR DIMENSIONED AT LEAST N AND C LL IS A LOGICAL VARIABLE. NONE OF THESE PASSED PARAMETERS C NEEDS TO BE DEFINED ON INPUT. ON OUTPUT IPERM CONTAINS THE C CURRENT PERMUTATION OF THE FIRST N INTEGERS AND LL IS .TRUE. C UNLESS THIS PERMUTATION IS THE LAST PERMUTATION OF THE CYCLE C (IN WHICH CASE LL IS .FALSE.). C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED JANUARY 20, 1987 C C REFERENCE: NIJENHUIS, ALBERT AND WILF, HERBERT S., 'COMBINATORIAL C ALGORITHMS', ACADEMIC PRESS, 1975, PP. 49-59. C----------------------------------------------------------------------- IMPLICIT INTEGER (A-Z) LOGICAL BROWN DIMENSION RED(*) DATA MAROON / 0 / IF (PINK.EQ.MAROON) GO TO 40 10 MAROON = PINK ORANGE = 1 SILVER = 1 PURPLE = 1 DO 20 BLUE = 1, PINK PURPLE = PURPLE*BLUE RED(BLUE) = BLUE 20 CONTINUE 30 BROWN = ORANGE.NE.PURPLE RETURN 40 IF (.NOT.BROWN) GO TO 10 GO TO (50,60), SILVER 50 GOLD = RED(2) RED(2) = RED(1) RED(1) = GOLD SILVER = 2 ORANGE = ORANGE+1 GO TO 30 60 YELLOW = 3 BLACK = ORANGE/2 70 VIOLET = MOD(BLACK,YELLOW) IF (VIOLET.NE.0) GO TO 80 BLACK = BLACK/YELLOW YELLOW = YELLOW+1 GO TO 70 80 BLACK = PINK GREEN = YELLOW-1 DO 90 BLUE = 1, GREEN WHITE = RED(BLUE)-RED(YELLOW) IF (WHITE.LT.0) WHITE = WHITE+PINK IF (WHITE.GE.BLACK) GO TO 90 BLACK = WHITE INDIGO = BLUE 90 CONTINUE GOLD = RED(YELLOW) RED(YELLOW) = RED(INDIGO) RED(INDIGO) = GOLD SILVER = 1 ORANGE = ORANGE+1 RETURN END *PERMAN SUBROUTINE PERMAN (A,LDA,N,IWORK,WORK,APERM) C C----------------------------------------------------------------------- C PERMAN COPIED BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 FROM THE REFERENCE BELOW C C FOR: COMPUTING THE PERMANENT OF THE N BY N MATRIX A. THE PERMANENT C OF A MATRIX IS SIMILAR TO THE DETERMINANT EXCEPT THAT THE C ALTERNATING SIGN CHANGES FOR THE TERMS ARE EXCLUDED. FOR C EXAMPLE, GIVEN THE MATRIX C C [A B C] C [D E F] C [G H I] C C THE DETERMINANT IS AEI+BFG+CDH-CEG-BDI-AFH WHEREAS THE C PERMANENT IS AEI+BFG+CDH+CEG+BDI+AFH . C C NOTE: THE COMPUTING TIME IS PROPORTIONAL TO 2**N, AND ON THE C CYBER 180/855 AT NBS A VALUE OF N=20 REQUIRES ABOUT 30 C SECONDS OF CPU TIME. C C SUBPROGRAMS CALLED: NEXSUB (ATTACHED) C C CURRENT VERSION COMPLETED JANUARY 20, 1987 C C REFERENCE: NIJENHUIS, ALBERT AND WILF, HERBERT S., 'COMBINATORIAL C ALGORITHMS', ACADEMIC PRESS, 1975, CH. 1,19. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * A = MATRIX (SIZE NXN) WHOSE PERMANENT IS TO BE COMPUTED (REAL) C C * LDA = THE LEADING DIMENSION OF MATRIX A [LDA>=N] (INTEGER) C C * N = THE NUMBER OF ROWS AND COLUMNS IN MATRIX A (INTEGER) C C IWORK = VECTOR (LENGTH N) USED AS SCRATCH AREA (INTEGER) C C WORK = VECTOR (LENGTH N) USED AS SCRATCH AREA (REAL) C C APERM = THE PERMANENT OF MATRIX A (REAL) C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- LOGICAL MTC DIMENSION A(LDA,*),IWORK(*),WORK(*) P = 0.0 N1 = N-1 DO 20 I = 1, N SUM = 0.0 DO 10 J = 1, N SUM = SUM+A(I,J) 10 CONTINUE WORK(I) = A(I,N)-SUM/2.0 20 CONTINUE SGN = -1.0 30 SGN = -SGN PROD = SGN CALL NEXSUB (N1,IWORK,MTC,NCARD,J) IF (NCARD.EQ.0) GO TO 50 Z = 2.0*REAL(IWORK(J))-1.0 DO 40 I = 1, N WORK(I) = WORK(I)+Z*A(I,J) 40 CONTINUE 50 DO 60 I = 1, N PROD = PROD*WORK(I) 60 CONTINUE P = P+PROD IF (MTC) GO TO 30 APERM = 2.0*REAL(2*MOD(N,2)-1)*P RETURN END C C======================================================================= C SUBROUTINE NEXSUB (N,IWORK,MTC,NCARD,J) LOGICAL MTC DIMENSION IWORK(*) DATA NLAST / 0 / IF (N.EQ.NLAST) GO TO 30 10 M = 0 MTC = .TRUE. DO 20 I = 1, N IWORK(I) = 0 20 CONTINUE NCARD = 0 NLAST = N RETURN 30 IF (.NOT.MTC) GO TO 10 M = M+1 M1 = M J = 0 40 J = J+1 IF (MOD(M1,2).EQ.1) GO TO 50 M1 = M1/2 GO TO 40 50 L = IWORK(J) IWORK(J) = 1-L NCARD = NCARD+1-2*L MTC = NCARD.NE.1.OR.IWORK(N).EQ.0 RETURN END *PLOTCR SUBROUTINE PLOTCR (X,Y,CH,N,IPR) C C----------------------------------------------------------------------- C PLOTCR MODIFIED BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 FROM THE PROGRAM PLOTC WRITTEN BY JAMES C J. FILLIBEN FOR THE SOFTWARE PACKAGE DATAPAC. C C FOR: CREATING A SIMPLE ONE-PAGE PRINTER PLOT OF Y(I) VS. X(I) WITH C PLOT CHARACTERS CONTROLLED BY CH(I). THE PLOTTING GRID IS C 45 LINES BY 109 COLUMNS WITH THE ENTIRE PLOT TAKING 125 C COLUMNS. THE PLOT CHARACTER USED AT THE Ith PLOTTING POSITION C WILL BE: C C 1 IF CH(I) IS BETWEEN 0.5 AND 1.5 C 2 IF CH(I) IS BETWEEN 1.5 AND 2.5 C : C : C 9 IF CH(I) IS BETWEEN 8.5 AND 9.5 C 0 IF CH(I) IS BETWEEN 9.5 AND 10.5 C A IF CH(I) IS BETWEEN 10.5 AND 11.5 C B IF CH(I) IS BETWEEN 11.5 AND 12.5 C : C : C Y IF CH(I) IS BETWEEN 34.5 AND 35.5 C Z IF CH(I) IS BETWEEN 35.5 AND 36.5 C * IF CH(I) IS OUTSIDE 0.5 TO 36.5 C + TO INDICATE MULTIPLE POINTS. C C THE NUMBER OF CHARACTERS PLOTTED IS N. IF N <= 0 ON INPUT C THEN NO PLOT IS DONE AND NO WARNING IS PRINTED. NO LEADING C OR TRAILING BLANK LINES ARE PRINTED. C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED NOVEMBER 25, 1987 C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C X(*) = VECTOR (LENGTH N) OF X VALUES PLOTTED HORIZONTALLY (REAL) C C Y(*) = VECTOR (LENGTH N) OF Y VALUES PLOTTED VERTICALLY (REAL) C C CH(*) = VECTOR (LENGTH N) OF PLOT CHARACTER CONTROL VALUES (REAL) C C N = NUMBER OF PAIRS (X(I),Y(I)) TO PLOT (INTEGER) C C IPR = LOGICAL UNIT NUMBER OF OUTPUT PRINT FILE (INTEGER) C----------------------------------------------------------------------- C DIMENSION X(*),Y(*),CH(*),YLABEL(11) CHARACTER*1 GR(45,109),PL(37) DATA PL / '1','2','3','4','5','6','7','8','9','0','A','B','C','D', * 'E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T' * ,'U','V','W','X','Y','Z','*'/ C C--- RETURN IF N<=0 C IF (N.LE.0) RETURN C C--- DETERMINE THE VALUES TO BE LISTED ON THE LEFT VERTICAL AXIS C YMIN = Y(1) YMAX = YMIN DO 10 I = 1, N YMIN = AMIN1(YMIN,Y(I)) YMAX = AMAX1(YMAX,Y(I)) 10 CONTINUE DO 20 I = 1, 9 YLABEL(I) = YMAX-(REAL(I-1)/8.0)*(YMAX-YMIN) 20 CONTINUE C C--- DETERMINE THE VALUES TO BE LISTED ON THE BOTTOM HORIZONTAL AXIS C--- DETERMINE XMIN, XMAX, XMID, X25 (=THE 25% POINT), AND C--- X75 (=THE 75% POINT) C XMIN = X(1) XMAX = XMIN DO 30 I = 1, N XMIN = AMIN1(XMIN,X(I)) XMAX = AMAX1(XMAX,X(I)) 30 CONTINUE XMID = 0.5*(XMIN+XMAX) X25 = 0.75*XMIN+0.25*XMAX X75 = 0.25*XMIN+0.75*XMAX C C--- BLANK OUT THE GRAPH C DO 50 I = 1, 45 DO 40 J = 1, 109 GR(I,J) = ' ' 40 CONTINUE 50 CONTINUE C C--- PRODUCE THE VERTICAL AXES C DO 60 I = 3, 43 GR(I,5) = 'I' GR(I,109) = 'I' 60 CONTINUE DO 70 I = 3, 43, 5 GR(I,5) = '-' GR(I,109) = '-' 70 CONTINUE GR(3,1) = '=' GR(3,2) = 'M' GR(3,3) = 'A' GR(3,4) = 'X' GR(23,1) = '=' GR(23,2) = 'M' GR(23,3) = 'I' GR(23,4) = 'D' GR(43,1) = '=' GR(43,2) = 'M' GR(43,3) = 'I' GR(43,4) = 'N' C C--- PRODUCE THE HORIZONTAL AXES C DO 80 J = 7, 107 GR(1,J) = '-' GR(45,J) = '-' 80 CONTINUE DO 90 J = 7, 107, 25 GR(1,J) = 'I' GR(45,J) = 'I' 90 CONTINUE DO 100 J = 20, 107, 25 GR(1,J) = 'I' GR(45,J) = 'I' 100 CONTINUE C C--- DETERMINE THE (X,Y) PLOT POSITIONS C RATIOY = 40.0/(YMAX-YMIN) RATIOX = 100.0/(XMAX-XMIN) DO 110 I = 1, N MX = 7+NINT(RATIOX*(X(I)-XMIN)) MY = 43-NINT(RATIOY*(Y(I)-YMIN)) IARG = NINT(CH(I)) IF (IARG.LT.1.OR.IARG.GT.36) IARG = 37 IF (GR(MY,MX).EQ.' ') THEN GR(MY,MX) = PL(IARG) ELSE GR(MY,MX) = '+' ENDIF 110 CONTINUE C C--- WRITE OUT THE GRAPH C DO 120 I = 1, 45 IP2 = I+2 IFLAG = IP2-(IP2/5)*5 K = IP2/5 IF (IFLAG.NE.0) THEN WRITE (IPR,1000) (GR(I,J),J=1,109) ELSE WRITE (IPR,1010) YLABEL(K),(GR(I,J),J=1,109) ENDIF 120 CONTINUE WRITE (IPR,1020) XMIN,X25,XMID,X75,XMAX RETURN C 1000 FORMAT (1X,15X,109A1) 1010 FORMAT (1X,G15.4,109A1) 1020 FORMAT (1X,7X,G20.4,5X,G20.4,5X,G20.4,5X,G20.4,1X,G20.4) C END *PPFBET SUBROUTINE PPFBET (PR,P,Q,TOL,IFLAG,X) C C----------------------------------------------------------------------- C PPFBET WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: EVALUATING THE INVERSE CUMULATIVE DISTRIBUTION FUNCTION C (PERCENT POINT FUNCTION) OF THE BETA(P,Q) DISTRIBUTION. C FOR A GIVEN PROBABILITY PR. THE PERCENT POINT X IS COMPUTED C TO A SPECIFIED ABSOLUTE ACCURACY TOL WHEN POSSIBLE. THE C METHOD OF BRENT, AS DESCRIBED IN THE REFERENCES BELOW, IS C USED TO COMPUTE THE APPROXIMATE ZERO OF I(X,P,Q)-PR WHERE C I(X,P,Q) IS THE CUMULATIVE DISTRIBUTION FUNCTION OF THE C BETA(P,Q) DISTRIBUTION EVALUATED AT X. THIS METHOD DOES C NOT REQUIRE DERIVATIVES. C C NOTE: THE CONSTANT EPS FOR MACHINE FLOATING POINT PRECISION IS C MACHINE DEPENDENT. C C SUBPROGRAMS CALLED: CDFBET (BETA CUMULATIVE DISTRIBUTION FUNCTION) C C CURRENT VERSION COMPLETED OCTOBER 13, 1987 C C REFERENCES: C C 1) PRESS, WILLIAM H., FLANNERY, BRIAN P., TEUKOLSKY, SAUL A., C AND VETTERLING, WILLIAM T., 'NUMERICAL RECIPES - THE ART OF C SCIENTIFIC COMPUTING', CAMBRIDGE UNIVERSITY PRESS, 1986, C PP. 251-254. C C 2) BRENT, RICHARD P., 'ALGORITHMS FOR MINIMIZATION WITHOUT C DERIVATIVES', PRENTICE-HALL, 1973, CH. 3-4. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * PR = A PROBABILITY VALUE IN THE INTERVAL [0,1] (REAL) C C * P = THE FIRST PARAMETER (>0) OF THE BETA(P,Q) DISTRIBUTION C (REAL) C C * Q = THE SECOND PARAMETER (>0) OF THE BETA(P,Q) DISTRIBUTION C (REAL) C C * TOL = THE REQUIRED ACCURACY (>=1.0E-8) OF THE PERCENT POINT X C (REAL) C C IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION: C 0 -> NO ERRORS DETECTED C 1,2 -> ERROR FLAGS FROM SUBROUTINE CDFBET C 3 -> PR<0 OR PR>1 C 4 -> P<=0 OR Q<=0 C 5 -> TOL<1.0E-8 C 6 -> THE CDF'S AT THE ENDPOINTS HAVE THE SAME SIGN - NO C