*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 VALUE OF X IS DEFINED (THIS SHOULD NEVER OCCUR) C 7 -> MAXIMUM ITERATIONS EXCEEDED - CURRENT VALUE OF X C RETURNED C C X = THE COMPUTED PERCENT POINT (REAL) C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C C--- DEFINE MAXIMUM ITERATIONS AND MACHINE FLOATING POINT PRECISION C DATA ITMAX,EPS / 50,1.0E-12 / C C--- CHECK VALIDITY OF INPUT ARGUMENTS C IFLAG = 0 IF (PR.LT.0.0.OR.PR.GT.1.0) THEN IFLAG = 3 RETURN C ENDIF IF (AMIN1(P,Q).LE.0.0) THEN IFLAG = 4 RETURN C ENDIF IF (TOL.LT.1.0E-8) THEN IFLAG = 5 RETURN C ENDIF A = 0.0 B = 1.0 FA = -PR FB = 1.0-PR C C--- CHECK FOR ZERO BEING BRACKETED C IF (FB*FA.GT.0.0) THEN IFLAG = 6 RETURN C ENDIF FC = FB DO 10 ITER = 1, ITMAX IF (FB*FC.GT.0.0) THEN C C--- RENAME A,B,C AND ADJUST BOUNDING INTERVAL D C C = A FC = FA D = B-A E = D ENDIF IF (ABS(FC).LT.ABS(FB)) THEN A = B B = C C = A FA = FB FB = FC FC = FA ENDIF C C--- CONVERGENCE CHECK C TOL1 = 2.0*EPS*ABS(B)+0.5*TOL XM = 0.5*(C-B) IF (ABS(XM).LE.TOL1.OR.FB.EQ.0.0) THEN X = B RETURN C ENDIF IF (ABS(E).GE.TOL1.AND.ABS(FA).GT.ABS(FB)) THEN C C--- ATTEMPT INVERSE QUADRATIC INTERPOLATION C S = FB/FA IF (A.EQ.C) THEN U = 2.0*XM*S V = 1.0-S ELSE V = FA/FC R = FB/FC U = S*(2.0*XM*V*(V-R)-(B-A)*(R-1.0)) V = (V-1.0)*(R-1.0)*(S-1.0) ENDIF C C--- CHECK WHETHER IN BOUNDS C IF (U.GT.0.0) V = -V U = ABS(U) IF (2.0*U.LT.AMIN1(3.0*XM*V-ABS(TOL1*V),ABS(E*V))) THEN C C--- ACCEPT INTERPLATION C E = D D = U/V ELSE C C--- INTERPOLATION FAILED, USE BISECTION C D = XM E = D ENDIF ELSE C C--- BOUNDS DECREASING TOO SLOWLY, USE BISECTION C D = XM E = D ENDIF C C--- MOVE LAST BEST GUESS TO A C A = B FA = FB C C--- EVALUATE NEW TRIAL ZERO C IF (ABS(D).GT.TOL1) THEN B = B+D ELSE B = B+SIGN(TOL1,XM) ENDIF CALL CDFBET (B,P,Q,EPS,IFLAG,CDF) IF (IFLAG.NE.0) RETURN FB = CDF-PR 10 CONTINUE IFLAG = 7 X = B RETURN C END *QDASIM SUBROUTINE QDASIM (F,A,B,MXPOW2,ABSERR,AREA,IFLAG) C C----------------------------------------------------------------------- C QDASIM WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899. C C FOR: NUMERICAL INTEGRATION USING AN ADAPTIVE SIMPSON'S RULE. THE C NUMBER OF PARTITIONS IS DOUBLED UNTIL SUCCESSIVE INTEGRAL C VALUES AGREE TO ABSERR OR BETTER. C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED JUNE 9, 1986 C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C F = USER-DEFINED REAL FUNCTION WITH ONE ARGUMENT, I.E., F(X) C C A = LEFT ENDPOINT OF REGION OF INTEGRATION (REAL) C C B = RIGHT ENDPOINT OF REGION OF INTEGRATION (REAL) C C MXPOW2 = INTEGER INDICATING THAT THE MAXIMUM NUMBER OF PARTITIONS C TO BE USED IS 2**MXPOW2 C C ABSERR = ABSOLUTE ERROR TOLERANCE ON THE COMPUTED VALUE OF THE C INTEGRAL (REAL) C C AREA = COMPUTED VALUE OF THE INTEGRAL (REAL) C C IFLAG = INTEGER ERROR FLAG ON RETURN. INTERPRETATION: C 0 -> INTEGRAL COMPUTED C 1 -> INTEGRAL NOT COMPUTED SINCE A>B C 2 -> MAXIMUM PARTITION SIZE REACHED WITHOUT CONVERGENCE C----------------------------------------------------------------------- C IFLAG = 0 AREA = 0.0 IF (A.GT.B) THEN IFLAG = 1 ELSE H = 0.5*(B-A) S = 2.0*F(A+H) R = F(A) R = R+F(B) C = H/3.0 AOLD = C*(R+S+S) N = 2 DO 20 I = 1, MXPOW2 N = N+N H = 0.5*H R = R+S S = 0.0 DO 10 J = 1, N, 2 X = A+H*REAL(J) S = S+F(X) 10 CONTINUE S = S+S C = 0.5*C AREA = C*(R+S+S) IF (ABS(AREA-AOLD).LE.ABSERR) RETURN AOLD = AREA 20 CONTINUE IFLAG = 2 ENDIF RETURN C END *QDTANH SUBROUTINE QDTANH (F,A,B,N,AREA,IFLAG) C C----------------------------------------------------------------------- C QDTANH WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899. C C FOR: NUMERICAL INTEGRATION USING THE 'TANH RULE' AS DESCRIBED IN C THE REFERENCE BELOW. THIS RULE IS ESPECIALLY GOOD FOR SMOOTH C FUNCTIONS WHICH ARE INTEGRABLE BUT HAVE SINGULARITIES OR C UNBOUNDED FIRST DERIVATIVES AT THE ENDPOINTS. C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED JUNE 9, 1986 C C REFERENCE: HABER, SEYMOUR, 'THE TANH RULE FOR NUMERICAL C INTEGRATION', SIAM JOURNAL OF NUMERICAL ANALYSIS, C VOL. 14, NO. 4, SEPTEMBER 1977. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C F = USER-DEFINED REAL FUNCTION WITH ONE ARGUMENT, I.E., F(X) C C A = LEFT ENDPOINT OF REGION OF INTEGRATION (REAL) C C B = RIGHT ENDPOINT OF REGION OF INTEGRATION (REAL) C C N = AN INTEGER GREATER THAN ZERO WHICH DETERMINES THE NUMBER C OF ABSCISSAE AT WHICH THE FUNCTION F IS EVALUATED (2*N+1) C C AREA = COMPUTED VALUE OF THE INTEGRAL C C IFLAG = INTEGER ERROR FLAG ON RETURN. INTERPRETATION: C 0 -> INTEGRAL COMPUTED C 1 -> INTEGRAL NOT COMPUTED SINCE N<1 C 2 -> INTEGRAL NOT COMPUTED SINCE A>B C----------------------------------------------------------------------- C IFLAG = 0 AREA = 0.0 IF (N.LT.1) THEN IFLAG = 1 RETURN C ENDIF IF (A.GT.B) THEN IFLAG = 2 ELSE CON = 0.5/REAL(N) H = 4.0*ATAN(1.0)*SQRT(CON)-CON HWIDTH = 0.5*(B-A) X = 0.5*(A+B) AREA = F(X) DO 10 I = 1, N HI = H*REAL(I) Y = HWIDTH*TANH(HI) FSUM = F(X+Y) AREA = AREA+(FSUM+F(X-Y))/COSH(HI)**2 10 CONTINUE AREA = HWIDTH*H*AREA ENDIF RETURN C END *QDTRAP FUNCTION QDTRAP (F,A,B,N,FDA,FDB,IFLAG) C C--- CORRECTED TRAPEZOIDAL RULE FOR INTEGRATION C--- CPR JUNE 5, 1986 REF: CONTE & DEBOOR C H = (B-A)/REAL(N) FA = F(A) QDTRAP = 0.5*(FA+F(B)) DO 10 I = 1, N-1 X = A+H*REAL(I) QDTRAP = QDTRAP+F(X) 10 CONTINUE QDTRAP = H*QDTRAP IF (IFLAG.EQ.1) QDTRAP = QDTRAP+(FDA-FDB)*H*H/12.0 RETURN C END *QSMNMX SUBROUTINE QSMNMX (A,B,C,D,E,F,XA,XB,YA,YB,XMIN,YMIN,ZMIN,XMAX, * YMAX,ZMAX) C C----------------------------------------------------------------------- C QSMNMX WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: FINDING THE MINIMUM (ZMIN) AND MAXIMUM (ZMAX) VALUES OF THE C QUADRATIC SURFACE DEFINED BY THE FUNCTION C C G(X,Y) = A*X**2 + B*X*Y + C*Y**2 + D*X + E*Y + F C C OVER THE RECTANGLE XA<=X<=XB AND YA<=Y<=YB. THE X AND Y C VALUES CORRESPONDING TO ZMIN (XMIN,YMIN) AND ZMAX (XMAX,YMAX) C ARE RETURNED. IF XA>XB OR YA>YB ON INPUT THEN NOTHING IS C DONE. ALL PASSED PARAMETERS ARE OF TYPE REAL. C C SUBPROGRAMS CALLED: MINMAX (ATTACHED) C C CURRENT VERSION COMPLETED JANUARY 15, 1987 C----------------------------------------------------------------------- C G(X,Y) = X*(A*X+B*Y+D)+Y*(C*Y+E)+F C C--- DO NOTHING IF XA>XB OR YA>YB C IF (XA.GT.XB.OR.YA.GT.YB) RETURN C C--- INITIALIZE MINIMUM AND MAXIMUM VALUES OF FUNCTION G C ZMIN = 1.0E30 ZMAX = -ZMIN C C--- EVALUATE FUNCTION AT CORNERS (XA,YA), (XB,YA), (XB,YA), (XB,YB) C Z = G(XA,YA) CALL MINMAX (XA,YA,Z,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX) Z = G(XA,YB) CALL MINMAX (XA,YB,Z,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX) Z = G(XB,YA) CALL MINMAX (XB,YA,Z,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX) Z = G(XB,YB) CALL MINMAX (XB,YB,Z,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX) C C--- CHECK BOUNDARIES Y=YA AND Y=YB FOR CANDIDATE POINTS C IF (A.NE.0.0) THEN X = -0.5*(B*YA+D)/A IF (X.GT.XA.AND.X.LT.XB) THEN Z = G(X,YA) CALL MINMAX (X,YA,Z,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX) ENDIF X = -0.5*(B*YB+D)/A IF (X.GT.XA.AND.X.LT.XB) THEN Z = G(X,YB) CALL MINMAX (X,YB,Z,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX) ENDIF ENDIF C C--- CHECK BOUNDARIES X=XA AND X=XB FOR CANDIDATE POINTS C IF (C.NE.0.0) THEN Y = -0.5*(B*XA+E)/C IF (Y.GT.YA.AND.Y.LT.YB) THEN Z = G(XA,Y) CALL MINMAX (XA,Y,Z,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX) ENDIF Y = -0.5*(B*XB+E)/C IF (Y.GT.YA.AND.Y.LT.YB) THEN Z = G(XB,Y) CALL MINMAX (XB,Y,Z,XMIN,YMIN,ZMIN,XMAX,YMAX,ZMAX) ENDIF ENDIF C C--- CHECK INTERIOR XA 0. IF THIS CONDITION IS NOT MET THEN PROGRAM C EXECUTION IS TERMINATED. C C IF MIN(P,Q) > 1 THEN THE CHENG ALGORITHM BB IS USED. C IF 0 < MIN(P,Q) <= 1 THEN THE CHENG ALGORITHM BC IS USED. C C SUBPROGRAMS CALLED: RDUNI (STSPAC) - UNIFORM(0,1) GENERATOR C C CURRENT VERSION COMPLETED JULY 30, 1984 C C REFERENCE: KENNEDY, WILLIAM J. AND GENTLE, JAMES E., "STATISTICAL C COMPUTING", MARCEL DEKKER, INC., 1980, PP. 216-219 C C REFERENCE: CHENG, R.C.H., "GENERATING BETA VARIATES WITH C NONINTEGRAL SHAPE PARAMETERS", COMMUNICATIONS OF C THE ACM, VOL. 21, NO. 4, APRIL 1978, PP. 317-322. C----------------------------------------------------------------------- C H = AMIN1(P,Q) IF (H.GT.1.0) THEN C C C CHENG ALGORITHM (BB) C A = H B = AMAX1(P,Q) ALPHA = A+B BETA = SQRT((ALPHA-2.)/(2.*A*B-ALPHA)) GAMMA = A+1./BETA C C--- STEP 1 C 10 U1 = RDUNI(0) U2 = RDUNI(0) V = BETA*ALOG(U1/(1.-U1)) W = A*EXP(V) Z = U2*U1**2 R = GAMMA*V-1.3862944 S = A+R-W C C--- STEP 2 C IF (S+2.609438.GE.5.*Z) GO TO 20 C C--- STEP 3 C T = ALOG(Z) IF (S.GE.T) GO TO 20 C C--- STEP 4 C IF (R+ALPHA*ALOG(ALPHA/(B+W)).LT.T) GO TO 10 C C--- STEP 5 C 20 IF (A.EQ.P) THEN RDBETA = W/(B+W) ELSE RDBETA = B/(B+W) ENDIF ELSEIF (H.GT.0.0) THEN C C C CHENG ALGORITHM (BC) C A = AMAX1(P,Q) B = H ALPHA = A+B BETA = 1./B DELTA = 1.+A-B KAPPA1 = DELTA*(0.013888889+0.041666667*B)/(A*BETA-0.77777778) KAPPA2 = 0.25+(0.5+0.25/DELTA)*B C C--- STEP 1 C 30 U1 = RDUNI(0) U2 = RDUNI(0) IF (U1.GE.0.5) GO TO 40 C C--- STEP 2 C Y = U1*U2 Z = U1*Y IF (0.25*U2+Z-Y.GE.KAPPA1) THEN GO TO 30 C ELSE GO TO 50 C ENDIF C C--- STEP 3 C 40 Z = U2*U1**2 IF (Z.LE.0.25) THEN V = BETA*ALOG(U1/(1.-U1)) W = A*EXP(V) GO TO 60 C ENDIF C C--- STEP 4 C IF (Z.GE.KAPPA2) GO TO 30 C C--- STEP 5 C 50 V = BETA*ALOG(U1/(1.-U1)) W = A*EXP(V) IF (ALPHA*(ALOG(ALPHA/(B+W))+V)-1.3862944.LT.ALOG(Z)) GO TO 30 C C--- STEP 6 C 60 IF (A.EQ.P) THEN RDBETA = W/(B+W) ELSE RDBETA = B/(B+W) ENDIF ELSE C C C PARAMETERS OUT-OF-RANGE C C PRINT *,' *** EITHER P OR Q IS <= 0' PRINT *,' *** EXECUTION STOPPED IN FUNCTION RDBETA' STOP C ENDIF RETURN END *RDCHI2 FUNCTION RDCHI2 (DF) C C----------------------------------------------------------------------- C RDCHI2 WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C C FOR: GENERATING A RANDOM DEVIATE FROM THE CHI-SQUARED(DF) C DISTRIBUTION WHERE DF > 0 AND DOES NOT HAVE TO BE AN INTEGER. C A CHI-SQUARED(DF) DEVIATE IS EQUAL TO TWICE A GAMMA(DF/2) C DEVIATE. SEE REFERENCE BELOW. C C SUBPROGRAMS CALLED: RDGAMM (STSPAC) C RDNOR (STSPAC) C C CURRENT VERSION COMPLETED FEBRUARY 28, 1986 C C REFERENCE: KENNEDY, WILLIAM J. AND GENTLE, JAMES E., 'STATISTICAL C COMPUTING', MARCEL DEKKER, INC., 1980, PP. 209-216 C----------------------------------------------------------------------- RDCHI2 = 2.*RDGAMM(DF/2.) RETURN C END *RDF FUNCTION RDF (DF1,DF2) C C----------------------------------------------------------------------- C RDF WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C C FOR: GENERATING A RANDOM DEVIATE FROM THE F(DF1,DF2) DISTRIBUTION C WHERE DF1,DF2 > 0 CAN BE NON-INTEGERS. IF X IS A C BETA(DF1/2,DF2/2) DEVIATE THEN DF2*X/(DF1*(1-X)) IS AN C F(DF1,DF2) DEVIATE. SEE REFERENCE BELOW. C C SUBPROGRAMS CALLED: RDBETA (STSPAC) C C CURRENT VERSION COMPLETED FEBRUARY 28, 1986 C C REFERENCE: KENNEDY, WILLIAM J. AND GENTLE, JAMES E., 'STATISTICAL C COMPUTING', MARCEL DEKKER, INC., 1980, PP. 219-221 C----------------------------------------------------------------------- C P = DF1/2. Q = DF2/2. 10 X = RDBETA(P,Q) IF (X.GE.1.0) THEN PRINT *,' *** A BETA DEVIATE OF 1.0 CANNOT BE TRANSFORMED TO' PRINT *,' *** AN F DEVIATE --- GENERATE ANOTHER BETA DEVIATE' GO TO 10 C ELSE RDF = DF2*X/(DF1*(1.-X)) ENDIF RETURN END *RDGAMM FUNCTION RDGAMM (ALPHA) C C----------------------------------------------------------------------- C RDGAMM WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, MD C C FOR: GENERATING A RANDOM DEVIATE FROM THE GAMMA(ALPHA) DISTRIBUTION. C ONE OF THREE METHODS IS USED DEPENDING ON THE VALUE OF THE C PARAMETER ALPHA (WHICH DOES NOT HAVE TO BE AN INTEGER): C C VALUE OF A METHOD USED C ------------- ------------------ C 0 < ALPHA < 1 AHRENS (GS) [FIRST REFERENCE] C ALPHA = 1 LOG TRANSFORMATION C ALPHA > 1 DIETER-AHRENS (GD) [SECOND REFERENCE] C C IF A <= 0 AN ERROR MESSAGE IS PRINTED AND EXECUTION IS C TERMINATED. C C DESCRIPTIONS OF EACH OF THESE ALGORITHMS CAN BE FOUND IN C THE REFERENCE GIVEN BELOW. C C SUBPROGRAMS CALLED: RDUNI (STSPAC) - UNIFORM(0,1) GENERATOR C RDNOR (STSPAC) - NORMAL(0,1) GENERATOR C C CURRENT VERSION COMPLETED FEBRUARY 28, 1986 C C REFERENCE: KENNEDY, WILLIAM J. AND GENTLE, JAMES E., 'STATISTICAL C COMPUTING', MARCEL DEKKER, INC., 1980, PP. 209-216 C C REFERENCE: AHRENS, J.H. AND DIETER, U., 'GENERATING GAMMA VARIATES C BY A MODIFIED REJECTION TECHNIQUE', COMMUNICATIONS OF C THE ACM, VOL. 25, NO. 1, PP. 47-54, JANUARY 1982. C----------------------------------------------------------------------- C DIMENSION Q(8),A(8),E(6) DATA AP,APP / 0.,0. / DATA Q / +0.041666661,+0.020834040,+0.007970958,+0.001686911,- * 0.000775882,+0.001274709,-0.000506403,+0.000213943 / DATA A / +0.33333332,-0.24999995,+0.20000622,-0.16667748,+ * 0.14236572,-0.12438558,+0.12337954,-0.11275089 / DATA E / 1.,0.50000027,0.1666605,0.04171864,0.00813673,0.00172501 * / IF (ALPHA.GT.1.0) THEN C C C DIETER-AHRENS ALGORITHM (GD) C C--- STEP 1 C IF (ALPHA.NE.AP) THEN AP = ALPHA S2 = ALPHA-0.5 S = SQRT(S2) D = SQRT(32.0)-12.*S ENDIF C C--- STEP 2 C T = RDNOR(0) X = S+0.5*T IF (T.GE.0.0) THEN RDGAMM = X**2 RETURN C ENDIF C C--- STEP 3 C U = RDUNI(0) IF (D*U.LE.T**3) THEN RDGAMM = X**2 RETURN C ENDIF C C--- STEP 4 C IF (ALPHA.NE.APP) THEN APP = ALPHA Q0 = 0. DO 10 K = 1, 8 Q0 = Q0+Q(K)*ALPHA**(-K) 10 CONTINUE IF (ALPHA.LE.3.686) THEN B = 0.463+S+0.178*S2 SIGMA = 1.235 C = 0.195/S-0.079+0.16*S ELSEIF (ALPHA.LE.13.022) THEN B = 1.654+0.0076*S2 SIGMA = 1.68/S+0.275 C = 0.062/S+0.024 ELSE B = 1.77 SIGMA = 0.75 C = 0.1515/S ENDIF ENDIF C C--- STEP 5 C IF (X.LE.0.0) GO TO 30 C C--- STEP 6 C V = T/(S+S) IF (ABS(V).GT.0.25) THEN QQ = Q0-S*T+0.25*T**2+(S2+S2)*ALOG(1.+V) ELSE SUM = 0. DO 20 K = 1, 8 SUM = SUM+A(K)*V**K 20 CONTINUE QQ = Q0+(0.5*T**2)*SUM ENDIF C C--- STEP 7 C IF (ALOG(1.-U).LE.QQ) THEN RDGAMM = X**2 RETURN C ENDIF C C--- STEP 8 C 30 EE = -ALOG(RDUNI(0)) U = 2.*RDUNI(0)-1. T = B+SIGN(EE*SIGMA,U) C C--- STEP 9 C IF (T.LE.-0.71874484) GO TO 30 C C--- STEP 10 C V = T/(S+S) IF (ABS(V).GT.0.25) THEN QQ = Q0-S*T+0.25*T**2+(S2+S2)*ALOG(1.+V) ELSE SUM = 0. DO 40 K = 1, 8 SUM = SUM+A(K)*V**K 40 CONTINUE QQ = Q0+(0.5*T**2)*SUM ENDIF C C--- STEP 11 C IF (QQ.LE.0.0) GO TO 30 ET = EE-0.5*T**2 IF (QQ.LE.0.5) THEN FACTOR = 0. DO 50 K = 1, 6 FACTOR = FACTOR+E(K)*QQ**K 50 CONTINUE EXPROD = FACTOR*EXP(ET) ELSE EXPROD = EXP(QQ+ET)-EXP(ET) ENDIF IF (C*ABS(U).GT.EXPROD) GO TO 30 C C--- STEP 12 C RDGAMM = (S+0.5*T)**2 ELSEIF (ALPHA.EQ.1.0) THEN C C C LOG TRANSFORMATION C C RDGAMM = -ALOG(RDUNI(0)) ELSEIF (ALPHA.GT.0.0) THEN C C C AHRENS ALGORITHM (GS) C C--- STEP 1 C EE = 2.7182818 60 U = RDUNI(0) B = (EE+ALPHA)/EE P = B*U IF (P.GT.1.0) GO TO 70 C C--- STEP 2 C RDGAMM = P**(1.0/ALPHA) V = RDUNI(0) IF (V.GT.EXP(-RDGAMM)) GO TO 60 RETURN C C--- STEP 3 C 70 RDGAMM = -ALOG((B-P)/ALPHA) V = RDUNI(0) IF (V.GT.RDGAMM**(ALPHA-1.0)) GO TO 60 RETURN C ELSE PRINT *,' *** THE GAMMA PARAMETER MUST BE > 0' PRINT *,' *** EXECUTION STOPPED IN FUNCTION RDGAMM' STOP ENDIF RETURN C END *RDMNOR SUBROUTINE RDMNOR (AMU,SIG,LDSIG,N,LTF,ZM,IFLAG) C C----------------------------------------------------------------------- C RDMNOR WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: COMPUTING A VECTOR OF PSEUDO-RANDOM MULTIVARIATE NORMAL C DEVIATES WITH MEAN AMU AND VARIANCE SIG. BEFORE THE FIRST C CALL WITH A GIVEN SIG, THE LOGICAL VARIABLE LTF MUST BE C ASSIGNED THE VALUE .TRUE. SO THAT A CHOLESKY FACTORIZATION C OF SIG IS PERFORMED AS IN THE REFERENCE. THE RESULT IS A C LOWER TRIANGULAR MATRIX L, STORED IN THE LOWER TRIANGLE OF C SIG, SUCH THAT SIG=LL'. FURTHER CALL TO RDMNOR USE ONLY THE C LOWER TRIANGLE OF SIG (L) IN COMPUTING THE DEVIATES UNTIL THE C VALUE OF LTF IS RESET TO .TRUE., EVEN IS SIG IS REDEFINED. C C NOTE: BEFORE THE FIRST CALL TO THIS ROUTINE THE CALL C C Z = RDNOR(ISEED) C C SHOULD BE MADE IN ORDER TO INITIALIZE THE NORMAL RANDOM C NUMBER GENERATOR WHERE ISEED IS A POSITIVE INTEGER. THIS C ALLOWS THE USER TO ESTABLISH A REPEATABLE SEQUENCE OF C DEVIATES. C C SUBPROGRAMS CALLED: RDNOR (PSEUDO-RANDOM NORMAL GENERATOR) C C CURRENT VERSION COMPLETED MAY 15, 1987 C C REFERENCE: STEWART, G.W., 'INTRODUCTION TO MATRIX COMPUTATIONS', C ACADEMIC PRESS, ALGORITHM 3.9, P 142. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * AMU = MEAN VECTOR (LENGTH N) OF THE MULTIVARIATE NORMAL C DEVIATES ZM (REAL) C C * SIG = COVARIANCE MATRIX (SIZE NXN) OF THE MULTIVARIATE NORMAL C DEVIATES ZM (REAL) C C * LDSIG = THE LEADING DIMENSION OF MATRIX SIG (>=N) (INTEGER) C C * N = THE LENGTH OF THE VECTOR OF DEVIATES ZM (INTEGER) C C * LTF = AN INDICATOR VARIABLE FOR PERFORMING A CHOLESKY C FACTORIZATION OF A NEW COVARIANCE MATRIX SIG (LOGICAL) C C ZM = A PSEUDO-RANDOM MULTIVARIATE NORMAL VECTOR (LENGTH N) C WITH MEAN AMU AND VARIANCE SIG (REAL) C C IFLAG = AN ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION: C 0 -> NO ERRORS DETECTED C 1 -> THE MATRIX SIG IS NOT POSITIVE SEMIDEFINITE, THUS C CANNOT BE A COVARIANCE MATRIX - NO DEVIATE GENERATED C C * INDICATES VARIABLES REQUIRING INPUT VALUES C----------------------------------------------------------------------- DIMENSION AMU(*),SIG(LDSIG,*),ZM(*) LOGICAL LTF C C--- IF NEW MATRIX SIG, PERFORM CHOLESKY FACTORIZATION. SET ERROR C--- FLAG IF SIG IS NOT POSITIVE DEFINITE C IF (LTF) THEN DO 40 K = 1, N DO 20 I = 1, K-1 S = 0.0 DO 10 J = 1, I-1 S = S+SIG(I,J)*SIG(K,J) 10 CONTINUE SIG(K,I) = (SIG(K,I)-S)/SIG(I,I) 20 CONTINUE S = 0.0 DO 30 J = 1, K-1 S = S+SIG(K,J)**2 30 CONTINUE Q = SIG(K,K)-S IF (Q.LT.0.0) THEN IFLAG = 1 RETURN ELSE SIG(K,K) = SQRT(Q) ENDIF 40 CONTINUE LTF = .FALSE. ENDIF IFLAG = 0 C C--- COMPUTE N INDEPENDENT N(0,1) PSEUDO-RANDOM DEVIATES IN ZM C DO 50 I = 1, N ZM(I) = RDNOR(0) 50 CONTINUE C C--- COMPUTE THE PSEUDO-RANDOM MULTIVARIATE NORMAL DEVIATES IN ZM C DO 70 I = N, 1, -1 S = 0.0 DO 60 J = 1, I S = S+SIG(I,J)*ZM(J) 60 CONTINUE ZM(I) = AMU(I)+S 70 CONTINUE RETURN END *RDNF FUNCTION RDNF (DF1,DF2,ALAMB1,ALAMB2) C C----------------------------------------------------------------------- C RDNF WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C C FOR: GENERATING A DOUBLY NONCENTRAL F DEVIATE AS THE RATIO OF TWO C NONCENTRAL CHI-SQUARED DEVIATES WITH DEGREES OF FREEDOM DF1 C AND DF2 (BOTH >= 1) AND NONCENTRALITY PARAMETERS ALAMB1 AND C ALAMB2 (BOTH >= 0). C C SUBPROGRAMS CALLED: RDCHI2 (STSPAC) - NONCENTRAL CHI-SQUARED DEVIATE C RDNOR (STSPAC) - NORMAL(0,1) DEVIATE C C CURRENT VERSION COMPLETED FEBRUARY 28, 1986 C----------------------------------------------------------------------- IF (DF1.GE.1.0.AND.DF2.GE.1.0) GO TO 10 PRINT *,' *** BOTH DEGREES OF FREEDOM MUST BE >= 1.0' PRINT *,' *** EXECUTION STOPPED IN FUNCTION RDNF' STOP 10 X1=(RDNOR(0)+SQRT(ALAMB1))**2 IF(DF1.GT.1.0) X1=X1+RDCHI2(DF1-1.0) X2=(RDNOR(0)+SQRT(ALAMB2))**2 IF(DF2.GT.1.0) X2=X2+RDCHI2(DF2-1.0) RDNF = DF2*X1/(DF1*X2) RETURN C END *RDNOR FUNCTION RDNOR (JD) C C----------------------------------------------------------------------- C RDNOR OBTAINED BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MD 20899 FROM THE SOFTWARE LIBRARY CMLIB AT NBS C C FOR: GENERATING NORMAL PSEUDO-RANDOM DEVIATES WITH MEAN ZERO AND C VARIANCE ONE. THE INITIAL CALL TO RDNOR SHOULD BE OF THE C FORM Z=RDNOR(K) WHERE K IS A POSITIVE INTEGER. FURTHER C CALLS SHOULD BE OF THE FORM Z=RDNOR(0). MDIG IS A LOWER C BOUND ON THE NUMBER OF BINARY DIGITS AVAILABLE FOR C REPRESENTING INTEGERS, INCLUDING THE SIGN BIT. THIS VALUE C MUST BE AT LEAST 16. TO ALLOW THE LONGEST SEQUENCE OF C RANDOM DEVIATES WITHOUT CYCLING, SET MDIG TO ITS MAXIMUM C VALUE. C C THE NBS CENTRAL COMPUTERS WERE FOUND TO GENERATE DEVIATES C AT THE FOLLOWING RATES USING THE INDICATED MAXIMUM VALUES C OF MDIG: C C CYBER 180/855 -> 54,893 DEV/SEC (MDIG=49) C CYBER 205/622 -> 89,098 DEV/SEC (MDIG=48) C C SUBPROGRAMS CALLED: RDUNI C C CURRENT VERSION COMPLETED FEBRUARY 9, 1987 C C REFERENCE: MARSAGLIA & TSANG, 'FAST, EASILY IMPLEMENTED METHOD C FOR SAMPLING FROM DECREASING OR SYMMETRIC UNIMODAL C DENSITY FUNCTIONS', SIAM J SISC 1983. C----------------------------------------------------------------------- C REAL V(65),W(65) INTEGER M(17) SAVE I1,J1,M,M1,M2,RMAX DATA AA,B,C,RMAX / 12.37586,.4878992,12.67706,3.0518509E-5 / DATA C1,C2,PC,XN / 0.9689279,1.301198,.1958303E-1,2.776994 / DATA V / 0.3409450,0.4573146,0.5397793,0.6062427,0.6631691, * 0.7136975,0.7596125,0.8020356,0.8417227,0.8792102,0.9148948, * 0.9490791,0.9820005,1.0138492,1.0447810,1.0749254,1.1043917, * 1.1332738,1.1616530,1.1896010,1.2171815,1.2444516,1.2714635, * 1.2982650,1.3249008,1.3514125,1.3778399,1.4042211,1.4305929, * 1.4569915,1.4834526,1.5100121,1.5367061,1.5635712,1.5906454, * 1.6179680,1.6455802,1.6735255,1.7018503,1.7306045,1.7598422, * 1.7896223,1.8200099,1.8510770,1.8829044,1.9155830,1.9492166, * 1.9839239,2.0198430,2.0571356,2.0959930,2.1366450,2.1793713, * 2.2245175,2.2725185,2.3239338,2.3795007,2.4402218,2.5075117, * 2.5834658,2.6713916,2.7769943,2.7769943,2.7769943,2.7769943 / DATA W / 0.10405134E-04,0.13956560E-04,0.16473259E-04,0.18501623E- * 04,0.20238931E-04,0.21780983E-04,0.23182241E-04,0.24476931E-04, * 0.25688121E-04,0.26832186E-04,0.27921226E-04,0.28964480E-04, * 0.29969191E-04,0.30941168E-04,0.31885160E-04,0.32805121E-04, * 0.33704388E-04,0.34585827E-04,0.35451919E-04,0.36304851E-04, * 0.37146564E-04,0.37978808E-04,0.38803170E-04,0.39621114E-04, * 0.40433997E-04,0.41243096E-04,0.42049621E-04,0.42854734E-04, * 0.43659562E-04,0.44465208E-04,0.45272764E-04,0.46083321E-04, * 0.46897980E-04,0.47717864E-04,0.48544128E-04,0.49377973E-04, * 0.50220656E-04,0.51073504E-04,0.51937936E-04,0.52815471E-04, * 0.53707761E-04,0.54616606E-04,0.55543990E-04,0.56492112E-04, * 0.57463436E-04,0.58460740E-04,0.59487185E-04,0.60546402E-04, * 0.61642600E-04,0.62780711E-04,0.63966581E-04,0.65207221E-04, * 0.66511165E-04,0.67888959E-04,0.69353880E-04,0.70922996E-04, * 0.72618816E-04,0.74471933E-04,0.76525519E-04,0.78843526E-04, * 0.81526890E-04,0.84749727E-04,0.84749727E-04,0.84749727E-04, * 0.84749727E-04 / DATA M / 30788,23052,2053,19346,10646,19427,23975,19049,10949, * 19693,29746,26748,2796,23890,29168,31924,16499 / DATA M1,M2,I1,J1 / 32767,256,5,17 / DATA MDIG / 49 / IF (JD.NE.0) GO TO 40 C C--- EXECUTE THIS SEGMENT EACH CALL C 10 I = M(I1)-M(J1) IF (I.LT.0) I = I+M1 M(J1) = I I1 = I1-1 IF (I1.EQ.0) I1 = 17 J1 = J1-1 IF (J1.EQ.0) J1 = 17 J = MOD(I,64)+1 RDNOR = I*W(J+1) IF (((I/M2)/2)*2.EQ.(I/M2)) RDNOR = -RDNOR IF (ABS(RDNOR).LE.V(J)) RETURN C C--- SLOW PART - AA IS A*F(0) C X = (ABS(RDNOR)-V(J))/(V(J+1)-V(J)) Y = RDUNI(0) S = X+Y IF (S.GT.C2) GO TO 30 IF (S.LE.C1) RETURN IF (Y.GT.C-AA*EXP(-0.5*(B-B*X)**2)) GO TO 30 IF (EXP(-0.5*V(J+1)**2)+Y*PC/V(J+1).LE.EXP(-0.5*RDNOR**2)) RETURN C C--- TAIL PART - 3.855849 IS 0.5*XN**2 C 20 S = XN-ALOG(RDUNI(0))/XN IF (3.855849+ALOG(RDUNI(0))-XN*S.GT.-0.5*S**2) GO TO 20 RDNOR = SIGN(S,RDNOR) RETURN 30 RDNOR = SIGN(B-B*X,RDNOR) RETURN C C--- EXECUTE THIS SEGMENT ONLY IF JD<>0 C 40 M1 = 2**(MDIG-2)+(2**(MDIG-2)-1) M2 = 2**(MDIG/2) JSEED = MIN0(IABS(JD),M1) IF (MOD(JSEED,2).EQ.0) JSEED = JSEED-1 K0 = MOD(9069,M2) K1 = 9069/M2 J0 = MOD(JSEED,M2) J1 = JSEED/M2 DO 50 I = 1, 17 JSEED = J0*K0 J1 = MOD(JSEED/M2+J0*K1+J1*K0,M2/2) J0 = MOD(JSEED,M2) 50 M(I) = J0+M2*J1 J1 = 17 I1 = 5 RMAX = 1.0/REAL(M1) C C--- SEED UNIFORM (0,1) GENERATOR (JUST A DUMMY CALL) C RDNOR = RDUNI(JD) DO 60 I = 1, 65 60 W(I) = RMAX*V(I) GO TO 10 END *RDNOR3 FUNCTION RDNOR3 (IX,IY,IZ) C C----------------------------------------------------------------------- C RDNOR3 WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: GENERATING STANDARD NORMAL PSEUDO-RANDOM DEVIATES BY THE C POLAR METHOD AS DESCRIBED IN THE REFERENCE BELOW. UNIFORM C PSEUDO-RANDOM DEVIATES IN THE INTERVAL (0,1) ARE SUPPLIED C BY A WICHMAN-HILL GENERATOR WHICH REQUIRES THREE INTEGER C ARGUMENTS IX, IY, AND IZ. THE NBS CENTRAL COMPUTERS WERE C FOUND TO GENERATE NORMAL DEVIATES AT THE FOLLOWING RATES: C C CYBER 180/855 -> 22,546 DEVIATES/SECOND C CYBER 205/622 -> 42,997 DEVIATES/SECOND . C C BEFORE THE FIRST CALL TO THIS ROUTINE THE INTEGERS IX, IY, AND C IZ SHOULD BE SET TO VALUES WITHIN THE RANGE [1,30000]. AFTER C THAT THEIR VALUES SHOULD NOT BE CHANGED EXCEPT BY THIS ROUTINE. C C SUBPROGRAMS CALLED: RDUNWH (WICHMAN-HILL UNIFORM GENERATOR) C C CURRENT VERSION COMPLETED MAY 15, 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. 117-118. C----------------------------------------------------------------------- C LOGICAL LL SAVE LL DATA LL / .FALSE. / IF (LL) THEN RDNOR3 = V*Q ELSE 10 U = 2.0*RDUNWH(IX,IY,IZ)-1.0 V = 2.0*RDUNWH(IX,IY,IZ)-1.0 S = U*U+V*V IF (S.GE.1.0) GO TO 10 Q = SQRT(-2.0*ALOG(S)/S) RDNOR3 = U*Q ENDIF LL = .NOT.LL RETURN END *RDT FUNCTION RDT (DF) C C----------------------------------------------------------------------- C RDT WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C C FOR: GENERATING A RANDOM DEVIATE FROM THE T(DF) DISTRIBUTION. C ONE OF THREE METHODS IS USED DEPENDING ON THE VALUE OF THE C PARAMETER DF (WHICH DOES NOT HAVE TO BE AN INTEGER): C C VALUE OF DF METHOD USED C ------------- ------------------ C 0 < DF < 1 NORMAL/SQRT(CHI-SQUARED/DF) C DF = 1 TANGENT TRANSFORMATION (CST) C DF > 1 KINDERMAN-MONAHAN-RAMAGE (TIR) C C IF DF <= 0 AN ERROR MESSAGE IS PRINTED AND EXECUTION IS C TERMINATED. C C DESCRIPTIONS OF EACH OF THESE ALGORITHMS CAN BE FOUND IN C THE REFERENCE GIVEN BELOW. C C SUBPROGRAMS CALLED: RDUNI (STSPAC) - UNIFORM(0,1) GENERATOR C RDNOR (STSPAC) - NORMAL(0,1) GENERATOR C RDCHI2 (STSPAC) - CHI-SQUARED GENERATOR C C CURRENT VERSION COMPLETED FEBRUARY 28, 1986 C C REFERENCE: KINDERMAN, A.J., MONAHAN, J.F., AND RAMAGE, J.G., C "COMPUTER METHODS FOR SAMPLING FROM STUDENT'S T C DISTRIBUTION", MATHEMATICS OF COMPUTATION, VOLUME 31, C NUMBER 140, OCTOBER 1977, PP. 1009-1018 C----------------------------------------------------------------------- C F(X,A) = (1.0+X*X/A)**(-(A+1.0)/2.0) IF (DF.GT.1.0) THEN C C C KINDERMAN-MONAHAN-RAMAGE ALGORITHM (TIR) C C--- STEP 1 C 10 U = RDUNI(0) IF (U.GE.0.23079283) GO TO 20 RDT = 4.0*U-0.46158566 C C--- STEP 2 C V = RDUNI(0) IF (V.LE.1.0-0.5*ABS(RDT)) RETURN IF (V.LE.F(RDT,DF)) RETURN GO TO 10 C C--- STEP 3 C 20 IF (U.GE.0.5) GO TO 40 S = 4.0*U-1.46158566 RDT = SIGN(ABS(S)+0.46158566,S) V = RDUNI(0) C C--- STEP 4 C 30 IF (V.LE.1.0-0.5*ABS(RDT)) RETURN IF (V.GE.1.2130613/(1.0+RDT*RDT)) GO TO 10 IF (V.LE.F(RDT,DF)) RETURN GO TO 10 C C--- STEP 5 C 40 IF (U.GE.0.75) GO TO 50 S = 8.0*U-5.0 RDT = 2.0/SIGN(ABS(S)+1.0,S) V = RDUNI(0)/(RDT*RDT) GO TO 30 C C--- STEP 6 C 50 RDT = 2.0/(8.0*U-7.0) V = RDUNI(0) IF (V.LT.RDT*RDT*F(RDT,DF)) RETURN GO TO 10 C ELSEIF (DF.EQ.1.0) THEN C C C SYNTHETIC TANGENT ALGORITHM (CST) C C--- STEP 1 C 60 U = RDUNI(0) V = 2.0*RDUNI(0)-1.0 C C--- STEP 2 C IF (U*U+V*V.GT.1.0) GO TO 60 RDT = V/U RETURN C ELSEIF (DF.GT.0.0) THEN C C C RATIO OF STANDARD NORMAL AND SQUARE ROOT OF C CHI-SQUARED DIVIDED BY ITS DEGREES OF FREEDOM C C D = SQRT(RDCHI2(DF)/DF) RDT = RDNOR(0)/D ELSE PRINT *,' *** DEGREES OF FREEDOM MUST BE > 0' PRINT *,' *** EXECUTION STOPPED IN FUNCTION RDT' STOP C ENDIF RETURN END *RDUNI FUNCTION RDUNI (JD) C C----------------------------------------------------------------------- C RDUNI OBTAINED BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MD 20899 FROM THE SOFTWARE LIBRARY CMLIB AT NBS C C FOR: GENERATING UNIFORM PSEUDO-RANDOM DEVIATES IN THE INTERVAL C [0,1) USING A LAGGED FIBONACCI GENERATOR. THE INITIAL CALL C TO RDUNI SHOULD BE OF THE FORM Z=RDUNI(K) WHERE K IS A C POSITIVE INTEGER. FURTHER CALLS SHOULD BE OF THE FORM C Z=RDUNI(0). MDIG IS A LOWER BOUND ON THE NUMBER OF BINARY C DIGITS AVAILABLE FOR REPRESENTING INTEGERS, INCLUDING THE C SIGN BIT. THIS VALUE MUST BE AT LEAST 16. TO ALLOW THE C LONGEST SEQUENCE OF RANDOM NUMBERS WITHOUT CYCLING, SET C MDIG TO ITS MAXIMUM VALUE. C C THE NBS CENTRAL COMPUTERS WERE FOUND TO GENERATE DEVIATES C AT THE FOLLOWING RATES USING THE INDICATED MAXIMUM VALUES C OF MDIG: C C CYBER 180/855 -> 100,050 DEV/SEC (MDIG=49) C CYBER 205/622 -> 105,471 DEV/SEC (MDIG=48) C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED FEBRUARY 9, 1987 C C REFERENCE: MARSAGLIA, GEORGE, 'COMMENTS ON THE PERFECT UNIFORM C RANDOM NUMBER GENERATOR', UNPUBLISHED NOTES, C WASHINGTON STATE UNIVERSITY. C----------------------------------------------------------------------- C INTEGER M(17) SAVE I,J,M,M1,M2 DATA M / 30788,23052,2053,19346,10646,19427,23975,19049,10949, * 19693,29746,26748,2796,23890,29168,31924,16499 / DATA M1,M2,I,J / 32767,256,5,17 / DATA MDIG / 49 / IF (JD.EQ.0) GO TO 20 C C--- EXECUTE THIS SEGMENT ONLY IF JD<>0 C M1 = 2**(MDIG-2)+(2**(MDIG-2)-1) M2 = 2**(MDIG/2) JSEED = MIN0(IABS(JD),M1) IF (MOD(JSEED,2).EQ.0) JSEED = JSEED-1 K0 = MOD(9069,M2) K1 = 9069/M2 J0 = MOD(JSEED,M2) J1 = JSEED/M2 DO 10 I = 1, 17 JSEED = J0*K0 J1 = MOD(JSEED/M2+J0*K1+J1*K0,M2/2) J0 = MOD(JSEED,M2) M(I) = J0+M2*J1 10 CONTINUE I = 5 J = 17 C C--- EXECUTE THIS SEGMENT EACH CALL C 20 K = M(I)-M(J) IF (K.LT.0) K = K+M1 M(J) = K I = I-1 IF (I.EQ.0) I = 17 J = J-1 IF (J.EQ.0) J = 17 RDUNI = REAL(K)/REAL(M1) RETURN END *RDUNLL FUNCTION RDUNLL (ISEED) C C----------------------------------------------------------------------- C RDUNLL WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: GENERATING UNIFORM PSEUDO-RANDOM DEVIATES IN THE INTERVAL C (0,1) USING A CONGRUENTIAL GENERATOR. THIS GENERATOR WAS C DEVELOPED AT LAWRENCE LIVERMORE LABORATORIES UNDER THE NAME C 'LLRANDOMI' AND IS CURRENTLY USED IN THE IMSL LIBRARY. IT C HAS A CYCLE LENGTH OF 2**31-2 = 2,147,483,646. BEFORE THE C FIRST CALL TO THIS ROUTINE ISEED SHOULD BE ASSIGNED A VALUE C IN THE RANGE [1,2147483646]. THEREAFTER THE VALUE OF ISEED C SHOULD NOT BE CHANGED. C C THE NBS CENTRAL COMPUTERS WERE FOUND TO GENERATE DEVIATES C AT THE FOLLOWING RATES USING THIS ROUTINE: C C CYBER 180/855 -> 125,707 DEV/SEC C CYBER 205/622 -> 260,373 DEV/SEC C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED FEBRUARY 10, 1987 C----------------------------------------------------------------------- C DATA M,RM / 2147483647,2147483647.0 / ISEED = MOD(16807*ISEED,M) RDUNLL = REAL(ISEED)/RM RETURN END *RDUNWH FUNCTION RDUNWH (IX,IY,IZ) C C----------------------------------------------------------------------- C RDUNWH COPIED BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 FROM THE REFERENCE BELOW WITH SLIGHT C MODIFICATIONS C C FOR: GENERATING UNIFORM PSEUDO-RANDOM DEVIATES IN THE INTERVAL C [0,1) BY THE WICHMAN-HILL METHOD. AN EXACT VALUE OF ZERO C IS POSSIBLE BUT NOT LIKELY. THE CYCLE LENGTH EXCEEDS 6.95E12. C TWO METHODS, WHICH GENERATE IDENTICAL DEVIATES, ARE AVAILABLE C FOR USE DEPENDING ON THE LIMITS OF THE INTEGER ARITHMETIC. C THE FIRST METHOD MAY BE USED ON COMPUTERS WITH 16-BIT WORDS. C THE METHOD NOT BEING USED CAN BE 'COMMENTED' OUT. USING THE C SECOND METHOD, THE NBS CENTRAL COMPUTERS WERE FOUND TO C GENERATE DEVIATES AT THE FOLLOWING RATES: C C CYBER 180/855 -> 48,995 DEVIATES/SECOND C CYBER 205/622 -> 126,564 DEVIATES/SECOND . C C BEFORE THE FIRST CALL TO THIS ROUTINE THE INTEGERS IX, IY, AND C IZ SHOULD BE SET TO VALUES WITHIN THE RANGE [1,30000]. AFTER C THAT THEIR VALUES SHOULD NOT BE CHANGED EXCEPT BY THIS ROUTINE. C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED MAY 15, 1987 C C REFERENCE: GRIFFITHS, P. AND HILL, I.D. (EDS.), 'APPLIED STATISTICS C ALGORITHMS', THE ROYAL STATISTICAL SOCIETY/ELLIS HARWOOD C LIMITED, 1985, ALGORITHM AS 183, PP. 238-242. C----------------------------------------------------------------------- C C--- THE NEXT SIX LINES REQUIRE INTEGER ARITHMETIC UP TO 30,323 C C IX = 171*MOD(IX,177)-2*(IX/177) C IY = 172*MOD(IY,176)-35*(IY/176) C IZ = 170*MOD(IZ,178)-63*(IZ/178) C IF (IX.LT.0) IX = IX+30269 C IF (IY.LT.0) IY = IY+30307 C IF (IZ.LT.0) IZ = IZ+30323 C C--- THE NEXT THREE LINES REQUIRE INTEGER ARITHMETIC UP TO 5,212,632 C IX = MOD(171*IX,30269) IY = MOD(172*IY,30307) IZ = MOD(170*IZ,30323) SUM = REAL(IX)/30269.0+REAL(IY)/30307.0+REAL(IZ)/30323.0 RDUNWH = AMOD(SUM,1.0) RETURN END *RDCONS SUBROUTINE RDCONS (P,LDP,N,H,PSI,IFLAG) C C----------------------------------------------------------------------- C RDCONS WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS ON THE C SURFACE OF A CONE WITH HEIGHT H AND HALF-ANGLE PSI. NO C POINTS WILL BE GENERATED UNLESS H>0 AND 00 ON INPUT, THEN POINTS WILL BE C GENERATED ON BOTH THE CURVED SURFACE AND THE BASE. C C POINTS ARE GENERATED ON THE CURVED SURFACE BY TRANSFORMATIONS, C AND ON THE BASE (IF REQUESTED) BY A REJECTION PROCEDURE C WHICH REJECTS 1 - PI/4 = 0.2146 OF THE POINTS. C C SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES C C CURRENT VERSION COMPLETED MAY 26, 1989 C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C P(LDP,*) = MATRIX (SIZE 3 BY N) WHOSE COLUMNS ARE THE N RANDOM C POINTS GENERATED ON THE CONE SURFACE. THE LEADING C DIMENSION OF P IS LDP (>=3) AND THE SECOND DIMENSION, C DEFINED IN THE CALLING PROGRAM, MUST BE AT LEAST N. C [REAL] C C * LDP = LEADING DIMENSION OF MATRIX P (>=3) [INTEGER] C C * N = NUMBER OF POINTS TO BE GENERATED [INTEGER] C C * H = HEIGHT OF THE CONE (>0) [REAL] C C * PSI = HALF-ANGLE OF THE CONE (>0) [REAL] C C * IFLAG = INDICATOR FOR WHERE TO GENERATE POINTS ON INPUT C (0 -> CURVED SURFACE ONLY; OTHERWISE -> CURVED SURFACE C AND BASE) [INTEGER] C C ERROR FLAG ON OUTPUT. INTERPRETATION: C 0 -> NO ERRORS DETECTED. C 1 -> H<=0 OR PSI<=0 OR PSI>=PI/2 C 2 -> N<1. C 3 -> LDP<3. C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION P(LDP,*) TWOPI = 8.0*ATAN(1.0) HALFPI = 0.25*TWOPI IF (H.LE.0.0.OR.PSI.LE.0.0.OR.PSI.GE.HALFPI) THEN IFLAG = 1 ELSEIF (N.LT.1) THEN IFLAG = 2 ELSEIF (LDP.LT.3) THEN IFLAG = 3 ELSE TANPSI = TAN(PSI) UC = 1.0/(1.0+SIN(PSI)) RR = H*TANPSI IF (IFLAG.EQ.0) THEN DO 10 I = 1, N Z = H*SQRT(RDUNI(0)) R = Z*TANPSI PHI = TWOPI*RDUNI(0) P(1,I) = R*COS(PHI) P(2,I) = R*SIN(PHI) P(3,I) = Z 10 CONTINUE ELSE DO 30 I = 1, N U1 = RDUNI(0) IF (U1.LE.UC) THEN Z = H*SQRT(U1/UC) R = Z*TANPSI PHI = TWOPI*RDUNI(0) P(1,I) = R*COS(PHI) P(2,I) = R*SIN(PHI) P(3,I) = Z ELSE 20 U1 = 2.0*RDUNI(0)-1.0 U2 = 2.0*RDUNI(0)-1.0 Q = U1**2+U2**2 IF (Q.GT.1.0) GO TO 20 P(1,I) = RR*U1 P(2,I) = RR*U2 P(3,I) = H ENDIF 30 CONTINUE ENDIF IFLAG = 0 ENDIF RETURN C END *RDCONV SUBROUTINE RDCONV (P,LDP,N,H,PSI,IFLAG) C C----------------------------------------------------------------------- C RDCONV WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS WITHIN THE C CONE WITH HEIGHT H AND HALF-ANGLE PSI. NO POINTS WILL BE C GENERATED UNLESS H>0 AND 0=3) AND THE SECOND DIMENSION, C DEFINED IN THE CALLING PROGRAM, MUST BE AT LEAST N. C [REAL] C C * LDP = LEADING DIMENSION OF MATRIX P (>=3) [INTEGER] C C * N = NUMBER OF POINTS TO BE GENERATED [INTEGER] C C * H = HEIGHT OF THE CONE (>0) [REAL] C C * PSI = HALF-ANGLE OF THE CONE (>0) [REAL] C C IFLAG = ERROR FLAG ON OUTPUT. INTERPRETATION: C 0 -> NO ERRORS DETECTED. C 1 -> H<=0 OR PSI<=0 OR H N<1. C 3 -> LDP<3. C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION P(LDP,*) TWOPI = 8.0*ATAN(1.0) HALFPI = 0.25*TWOPI IF (H.LE.0.0.OR.PSI.LE.0.0.OR.PSI.GT.HALFPI) THEN IFLAG = 1 ELSEIF (N.LT.1) THEN IFLAG = 2 ELSEIF (LDP.LT.3) THEN IFLAG = 3 ELSE IFLAG = 0 CUBERT = 1.0/3.0 TANPSI = TAN(PSI) DO 10 I = 1, N Z = H*RDUNI(0)**CUBERT R = Z*TANPSI*SQRT(RDUNI(0)) PHI = TWOPI*RDUNI(0) P(1,I) = R*COS(PHI) P(2,I) = R*SIN(PHI) P(3,I) = Z 10 CONTINUE ENDIF RETURN C END *RDCYLS SUBROUTINE RDCYLS (P,LDP,N,H,R,IFLAG) C C----------------------------------------------------------------------- C RDCYLS WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS ON THE C SURFACE OF THE CYLINDER WITH HEIGHT H AND RADIUS R. NO C POINTS WILL BE GENERATED UNLESS H>=0 AND R>=0. THE AXIS C OF THE CYLINDER COINCIDES WITH THE POSITIVE Z AXIS, AND THE C ENDS OF THE CYLINDER ARE FORMED BY THE PLANES Z=0 AND Z=H. C ON OUTPUT, EACH OF THE N COLUMNS OF MATRIX P CONTAINS THE X, C Y, AND Z COORDINATES OF THE GENERATED POINTS. THE LEADING C DIMENSION OF P (LDP) MUST BE 3 OR MORE IN THE CALLING PROGRAM. C C AT THE USERS OPTION, POINTS MAY BE GENERATED ON THE CURVED C SURFACE ONLY OR ON BOTH THE CURVED AND FLAT SURFACES. POINTS C ARE GENERATED ON THE CURVED SURFACE BY A TRANSFORMATION, AND C POINTS ARE GENERATED ON THE FLAT SURFACES BY A REJECTION C PROCEDURE WHEREBY THE PROPORTION OF POINTS REJECTED IS 0.2146 C = 1-PI/4. C C NOTE THAT POINTS CAN BE GENERATED ON A STRAIGHT LINE BY C SETTING H>0 AND R=0. POINTS CAN BE GENERATED ON A CIRCLE BY C SETTING H=0 AND R>0. C C SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES C C CURRENT VERSION COMPLETED MAY 30, 1989 C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C P(LDP,*) = MATRIX (SIZE 3 BY N) WHOSE COLUMNS ARE THE N RANDOM C POINTS GENERATED ON THE CYLINDER SURFACE. THE LEADING C DIMENSION OF P IS LDP (>=3) AND THE SECOND DIMENSION, C DEFINED IN THE CALLING PROGRAM, MUST BE AT LEAST N. C [REAL] C C * LDP = LEADING DIMENSION OF MATRIX P (>=3) [INTEGER] C C * N = NUMBER OF POINTS TO BE GENERATED [INTEGER] C C * H = HEIGHT OF THE CYLINDER (>=0) [REAL] C C * R = RADIUS OF THE CYLINDER (>=0) [REAL] C C * IFLAG = INDICATOR FOR WHERE TO GENERATE POINTS ON INPUT C (0 -> CURVED SURFACE ONLY; OTHERWISE -> CURVED SURFACE C AND BASE) [INTEGER] C C ERROR FLAG ON OUTPUT. INTERPRETATION: C 0 -> NO ERRORS DETECTED. C 1 -> H<0. C 2 -> R<0. C 3 -> LDP<3. C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION P(LDP,*) TWOPI = 8.0*ATAN(1.0) IF (H.LT.0.0) THEN IFLAG = 1 ELSEIF (R.LT.0.0) THEN IFLAG = 2 ELSE IF (IFLAG.EQ.0) THEN C C--- GENERATE POINTS ON THE CURVED SURFACE ONLY C DO 10 I = 1, N P(3,I) = H*RDUNI(0) U = TWOPI*RDUNI(0) P(1,I) = R*COS(U) P(2,I) = R*SIN(U) 10 CONTINUE ELSE C C--- GENERATE POINTS ON BOTH THE CURVED AND FLAT SURFACES C UC1 = H/(H+R) UC2 = 0.5*(1.0+UC1) DO 30 I = 1, N U = RDUNI(0) IF (U.LE.UC1) THEN P(3,I) = H*U/UC1 U = TWOPI*RDUNI(0) P(1,I) = R*COS(U) P(2,I) = R*SIN(U) ELSE 20 U1 = 2.0*RDUNI(0)-1.0 U2 = 2.0*RDUNI(0)-1.0 Q = U1*U1+U2*U2 IF (Q.GT.1.0) GO TO 20 P(1,I) = R*U1 P(2,I) = R*U2 IF (U.LE.UC2) THEN P(3,I) = 0.0 ELSE P(3,I) = H ENDIF ENDIF 30 CONTINUE ENDIF IFLAG = 0 ENDIF RETURN C END *RDCYLV SUBROUTINE RDCYLV (P,LDP,N,H,R,IFLAG) C C----------------------------------------------------------------------- C RDCYLV WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS WITHIN THE C CYLINDER WITH HEIGHT H AND RADIUS R. NO POINTS WILL BE C GENERATED UNLESS H>=0 AND R>=0. THE AXIS OF THE CYLINDER C COINCIDES WITH THE POSITIVE Z AXIS, AND THE ENDS OF THE C CYLINDER ARE FORMED BY THE PLANES Z=0 AND Z=H. ON OUTPUT, C EACH OF THE N COLUMNS OF MATRIX P CONTAINS THE X, Y, AND Z C COORDINATES OF THE GENERATED POINTS. THE LEADING DIMENSION C OF P (LDP) MUST BE 3 OR MORE IN THE CALLING PROGRAM. C C POINTS ARE GENERATED BY A REJECTION PROCEDURE WHEREBY THE C PROPORTION OF POINTS REJECTED IS 0.2146 = 1-PI/4. C C NOTE THAT POINTS CAN BE GENERATED (RATHER INEFFICIENTLY) ON C A STRAIGHT LINE BY SETTING H>0 AND R=0. POINTS CAN BE C GENERATED ON A (TWO-DIMENSIONAL) DISK BY SETTING H=0 AND R>0. C C SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES C C CURRENT VERSION COMPLETED MAY 30, 1989 C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C P(LDP,*) = MATRIX (SIZE 3 BY N) WHOSE COLUMNS ARE THE N RANDOM C POINTS GENERATED INSIDE THE CYLINDER. THE LEADING C DIMENSION OF P IS LDP (>=3) AND THE SECOND DIMENSION, C DEFINED IN THE CALLING PROGRAM, MUST BE AT LEAST N. C [REAL] C C * LDP = LEADING DIMENSION OF MATRIX P (>=3) [INTEGER] C C * N = NUMBER OF POINTS TO BE GENERATED [INTEGER] C C * H = HEIGHT OF THE CYLINDER (>=0) [REAL] C C * R = RADIUS OF THE CYLINDER (>=0) [REAL] C C IFLAG = ERROR FLAG ON OUTPUT. INTERPRETATION: C 0 -> NO ERRORS DETECTED. C 1 -> H<0. C 2 -> R<0. C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION P(LDP,*) IF (H.LT.0.0) THEN IFLAG = 1 ELSEIF (R.LT.0.0) THEN IFLAG = 2 ELSE IFLAG = 0 DO 20 I = 1, N P(3,I) = H*RDUNI(0) 10 U = 2.0*RDUNI(0)-1.0 V = 2.0*RDUNI(0)-1.0 Q = U*U+V*V IF (Q.GT.1.0) GO TO 10 P(1,I) = R*U P(2,I) = R*V 20 CONTINUE ENDIF RETURN C END *RDELLS SUBROUTINE RDELLS (P,LDP,N,M,A,IFLAG) C C----------------------------------------------------------------------- C RDELLS WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS ON THE C SURFACE OF THE M-DIMENSIONAL ELLIPSOID DEFINED BY C C M C SUM (X /A )**2 = 1 C J=1 J J C C WHERE EACH A(J)>0. THE CENTER OF THE ELLIPSOID IS AT THE C ORIGIN. NOTE THAT M STANDARD NORMAL DEVIATES AND 1 STANDARD C UNIFORM DEVIATE ARE REQUIRED TO GENERATE EACH POINT. IF THE C A(J) ARE NOT ALL EQUAL, THERE IS A POSITIVE PROBABILITY THAT C THE POINT WILL BE REJECTED. C C SUBPROGRAMS CALLED: RDNOR (STSPAC) - N(0,1) RANDOM DEVIATES C RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES C C CURRENT VERSION COMPLETED JANUARY 8, 1990 C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C P(LDP,*) = MATRIX (SIZE M BY N) WHOSE COLUMNS ARE THE N RANDOM C POINTS GENERATED ON THE SURFACE OF THE ELLIPSOID. THE C LEADING DIMENSION OF P IS LDP (>=M) AND THE SECOND C DIMENSION, DEFINED IN THE CALLING PROGRAM, MUST BE AT C LEAST N. [REAL] C C * LDP = LEADING DIMENSION OF MATRIX P (>=M) [INTEGER] C C * N = NUMBER OF POINTS TO BE GENERATED [INTEGER] C C * M = DIMENSION OF THE ELLIPSOID (>0) [INTEGER] C C * A(*) = VECTOR (LENGTH M) OF ELLIPSOID SEMIAXES (>0) [REAL] C C IFLAG = ERROR FLAG ON OUTPUT. INTERPRETATION: C 0 -> NO ERRORS DETECTED. C 1 -> N<1 OR M<1. C 2 -> M>LDP. C 3 -> AT LEAST ONE A(J)<=0. C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION A(*),P(LDP,*) IF (MIN0(N,M).LE.0) THEN IFLAG = 1 ELSEIF (M.GT.LDP) THEN IFLAG = 2 ELSE AMIN = A(1) DO 10 J = 2, M AMIN = AMIN1(AMIN,A(J)) 10 CONTINUE IF (AMIN.LE.0.0) THEN IFLAG = 3 ELSE IFLAG = 0 DO 50 I = 1, N 20 Q = 0.0 R = 0.0 DO 30 J = 1, M P(J,I) = RDNOR(0) Q = Q+P(J,I)**2 R = R+(P(J,I)/A(J))**2 30 CONTINUE Q = SQRT(Q) R = SQRT(R) U1 = RDUNI(0) IF (AMIN*R.LT.Q*U1) GO TO 20 DO 40 J = 1, M P(J,I) = A(J)*P(J,I)/Q 40 CONTINUE 50 CONTINUE ENDIF ENDIF RETURN C END *RDELLV SUBROUTINE RDELLV (P,LDP,N,M,A,IFLAG) C C----------------------------------------------------------------------- C RDELLV WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS INSIDE THE C M-DIMENSIONAL ELLIPSOID DEFINED BY C C M C SUM (X /A )**2 = 1 C J=1 J J C C WHERE EACH A(J)>0. THE CENTER OF THE ELLIPSOID IS AT THE C ORIGIN. NOTE THAT M STANDARD NORMAL DEVIATES AND ONE C UNIFORM DEVIATE ARE REQUIRED TO GENERATE EACH POINT. C C SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES C RDNOR (STSPAC) - N(0,1) RANDOM DEVIATES C C CURRENT VERSION COMPLETED MARCH 16, 1989 C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C P(LDP,*) = MATRIX (SIZE M BY N) WHOSE COLUMNS ARE THE N RANDOM C POINTS GENERATED INSIDE THE ELLIPSOID. THE C LEADING DIMENSION OF P IS LDP (>=M) AND THE SECOND C DIMENSION, DEFINED IN THE CALLING PROGRAM, MUST BE AT C LEAST N. [REAL] C C * LDP = LEADING DIMENSION OF MATRIX P (>=M) [INTEGER] C C * N = NUMBER OF POINTS TO BE GENERATED [INTEGER] C C * M = DIMENSION OF THE ELLIPSOID (>0) [INTEGER] C C * A(*) = VECTOR (LENGTH M) OF ELLIPSOID SEMIAXES (>0) [REAL] C C IFLAG = ERROR FLAG ON OUTPUT. INTERPRETATION: C 0 -> NO ERRORS DETECTED. C 1 -> N<1 OR M<1. C 2 -> M>LDP. C 3 -> AT LEAST ONE A(J)<=0. C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION A(*),P(LDP,*) IF (MIN0(N,M).LE.0) THEN IFLAG = 1 ELSEIF (M.GT.LDP) THEN IFLAG = 2 ELSE AMIN = A(1) DO 10 J = 2, M AMIN = AMIN1(AMIN,A(J)) 10 CONTINUE IF (AMIN.LE.0.0) THEN IFLAG = 3 ELSE IFLAG = 0 POW = 1.0/REAL(M) DO 40 I = 1, N Q = 0.0 DO 20 J = 1, M P(J,I) = RDNOR(0) Q = Q+P(J,I)**2 20 CONTINUE R = RDUNI(0)**POW FAC = R/SQRT(Q) DO 30 J = 1, M P(J,I) = FAC*A(J)*P(J,I) 30 CONTINUE 40 CONTINUE ENDIF ENDIF RETURN C END *RDHLX SUBROUTINE RDHLX (P,LDP,N,R,D,H,PHI,KAPPA,IFLAG) C C----------------------------------------------------------------------- C RDHLX WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS ON A HELIX C IN THREE-DIMENSIONAL SPACE. THE HELIX MAY BE 'RIGHT-HANDED' C OR 'LEFT-HANDED' AT THE USERS OPTION. THE AXIS OF THE HELIX C COINCIDES WITH THE Z-AXIS BETWEEN Z=0 AND Z=H. THE HELIX C RADIUS IS R, AND THE SEPARATION BETWEEN ADJACENT POINTS WITH C THE SAME (X,Y) COORDINATES IS D. THE INITIAL HELIX ANGLE C (AT Z=0) IS PHI. C C SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES C C CURRENT VERSION COMPLETED JANUARY 9, 1990 C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C P(LDP,*) = MATRIX (SIZE 3 BY N) WHOSE COLUMNS ARE THE N RANDOM C POINTS GENERATED ON THE HELIX. THE LEADING C DIMENSION OF P IS LDP (>=3) AND THE SECOND DIMENSION, C DEFINED IN THE CALLING PROGRAM, MUST BE AT LEAST N. C [REAL] C C * LDP = LEADING DIMENSION OF MATRIX P (>=3) [INTEGER] C C * N = NUMBER OF POINTS TO BE GENERATED (N>0) [INTEGER] C C * R = THE RADIUS OF THE HELIX (R>=0) [REAL] C C * D = THE Z-DISTANCE BETWEEN ADJACENT POINTS WHICH HAVE THE C SAME (X,Y) COORDINATES (D>0) [REAL] C C * H = THE HEIGHT OF THE HELIX (H>0) [REAL] C C * PHI = INITIAL HELIX ANGLE WHEN Z=0 [REAL] C C * KAPPA = 1 FOR A RIGHT-HANDED HELIX, C -1 FOR A LEFT-HANDED HELIX [INTEGER] C C IFLAG = ERROR FLAG ON OUTPUT. INTERPRETATION: C 0 -> NO ERRORS DETECTED. C 1 -> N<1 OR D<=0 OR H<=0. C 2 -> R<0. C 3 -> ABS(KAPPA)<>1. C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION P(LDP,*) IF (N.LT.1.OR.D.LE.0.0.OR.H.LE.0.0) THEN IFLAG = 1 ELSEIF (R.LT.0.0) THEN IFLAG = 2 ELSEIF (IABS(KAPPA).NE.1) THEN IFLAG = 3 ELSE IFLAG = 0 CON = 8.0*ATAN(1.0)*REAL(KAPPA)/D DO 10 I=1,N U1 = RDUNI(0) Z = H*U1 THETA = PHI+CON*Z P(1,I) = R*COS(THETA) P(2,I) = R*SIN(THETA) P(3,I) = Z 10 CONTINUE ENDIF RETURN C END *RDRECS SUBROUTINE RDRECS (P,LDP,N,A,B,C,IFLAG) C C----------------------------------------------------------------------- C RDRECS WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS ON THE C SURFACE OF A 3-DIMENSIONAL RECTANGULAR SOLID WITH SIDES OF C LENGTH A, B, AND C. NO POINTS WILL BE GENERATED UNLESS C MIN(A,B,C)>=0. ONE CORNER OF THE RECTANGULAR SOLID HAS C COORDINATES (0,0,0) AND THE OPPOSITE CORNER HAS COORDINATES C (A,B,C). ON OUTPUT, EACH OF THE N COLUMNS OF MATRIX P C CONTAINS THE X, Y, AND Z COORDINATES OF THE GENERATED POINTS. C THE LEADING DIMENSION OF P (LDP) MUST BE 3 OR MORE IN THE C CALLING PROGRAM. C C NOTE THAT POINTS MAY BE GENERATED IN A PLANE BY SETTING ONE C OF (A,B,C) TO ZERO, AND POINTS MAY BE GENERATED ON A LINE BY C SETTING TWO OF (A,B,C,) TO ZERO. C C SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES C C CURRENT VERSION COMPLETED MAY 31, 1989 C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C P(LDP,*) = MATRIX (SIZE 3 BY N) WHOSE COLUMNS ARE THE N RANDOM C POINTS GENERATED ON THE RECTANGULAR SOLID SURFACE. THE C LEADING DIMENSION OF P IS LDP (>=3) AND THE SECOND C DIMENSION, DEFINED IN THE CALLING PROGRAM, MUST BE AT C LEAST N. [REAL] C C * LDP = LEADING DIMENSION OF MATRIX P (>=3) [INTEGER] C C * N = NUMBER OF POINTS TO BE GENERATED [INTEGER] C C * A = LENGTH OF RECTANGULAR SOLID IN X DIRECTION (>=0) [REAL] C C * B = LENGTH OF RECTANGULAR SOLID IN Y DIRECTION (>=0) [REAL] C C * C = LENGTH OF RECTANGULAR SOLID IN Z DIRECTION (>=0) [REAL] C C IFLAG = ERROR FLAG ON OUTPUT. INTERPRETATION: C 0 -> NO ERRORS DETECTED. C 1 -> MIN(A,B,C)<0. C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION P(LDP,*) IF (AMIN1(A,B,C).LT.0.0) THEN IFLAG = 1 ELSE SA = A*(B+C)+B*C C1 = A*B/SA C2 = A*(B+C)/SA C C--- GENERATE POINTS ON THE SURFACE OF THE RECTANGULAR SOLID C DO 10 I = 1, N U = RDUNI(0) V = RDUNI(0) IF (U.LT.C1) THEN P(1,I) = A*U/C1 IF (V.LE.0.5) THEN P(2,I) = 2.0*B*V P(3,I) = 0.0 ELSE P(2,I) = 2.0*B*(V-0.5) P(3,I) = C ENDIF ELSEIF (U.LT.C2) THEN P(1,I) = A*(U-C1)/(C2-C1) IF (V.LE.0.5) THEN P(3,I) = 2.0*C*V P(2,I) = 0.0 ELSE P(3,I) = 2.0*C*(V-0.5) P(2,I) = B ENDIF ELSE P(2,I) = B*(U-C2)/(1.0-C2) IF (V.LE.0.5) THEN P(3,I) = 2.0*C*V P(1,I) = 0.0 ELSE P(3,I) = 2.0*C*(V-0.5) P(1,I) = A ENDIF ENDIF 10 CONTINUE IFLAG = 0 ENDIF RETURN C END *RDRECV SUBROUTINE RDRECV (P,LDP,N,A,B,C,IFLAG) C C----------------------------------------------------------------------- C RDRECV WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS INSIDE A C 3-DIMENSIONAL RECTANGULAR SOLID WITH SIDES OF LENGTH A, B, C AND C. NO POINTS WILL BE GENERATED UNLESS MIN(A,B,C)>=0. ONE C CORNER OF THE RECTANGULAR SOLID HAS COORDINATES (0,0,0) AND C THE OPPOSITE CORNER HAS COORDINATES (A,B,C). ON OUTPUT, C EACH OF THE N COLUMNS OF MATRIX P CONTAINS THE X, Y, AND Z C COORDINATES OF THE GENERATED POINTS. THE LEADING DIMENSION C OF P (LDP) MUST BE 3 OR MORE IN THE CALLING PROGRAM. C C NOTE THAT POINTS MAY BE GENERATED IN A PLANE BY SETTING ONE C OF (A,B,C) TO ZERO, AND POINTS MAY BE GENERATED ON A LINE BY C SETTING TWO OF (A,B,C,) TO ZERO. C C SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES C C CURRENT VERSION COMPLETED MAY 31, 1989 C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C P(LDP,*) = MATRIX (SIZE 3 BY N) WHOSE COLUMNS ARE THE N RANDOM C POINTS GENERATED INSIDE THE RECTANGULAR SOLID. THE C LEADING DIMENSION OF P IS LDP (>=3) AND THE SECOND C DIMENSION, DEFINED IN THE CALLING PROGRAM, MUST BE AT C LEAST N. [REAL] C C * LDP = LEADING DIMENSION OF MATRIX P (>=3) [INTEGER] C C * N = NUMBER OF POINTS TO BE GENERATED [INTEGER] C C * A = LENGTH OF RECTANGULAR SOLID IN X DIRECTION (>=0) [REAL] C C * B = LENGTH OF RECTANGULAR SOLID IN Y DIRECTION (>=0) [REAL] C C * C = LENGTH OF RECTANGULAR SOLID IN Z DIRECTION (>=0) [REAL] C C IFLAG = ERROR FLAG ON OUTPUT. INTERPRETATION: C 0 -> NO ERRORS DETECTED. C 1 -> MIN(A,B,C)<0. C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION P(LDP,*) IF (AMIN1(A,B,C).LT.0.0) THEN IFLAG = 1 ELSE C C--- GENERATE POINTS INSIDE A RECTANGULAR SOLID C DO 10 I = 1, N P(1,I) = A*RDUNI(0) P(2,I) = B*RDUNI(0) P(3,I) = C*RDUNI(0) 10 CONTINUE IFLAG = 0 ENDIF RETURN C END *RDSPSH SUBROUTINE RDSPSH (P,LDP,L,N,R1,R2,IFLAG) C C----------------------------------------------------------------------- C RDSPSH WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: GENERATING N RANDOM POINTS 'UNIFORMLY' SPACED IN A REGION C BETWEEN THE SURFACES OF TWO SPHERES OF RADIUS R1 AND R2. C THE ONLY REQUIREMENTS ARE THAT N>0, L>0, R1>=0, AND R2>=0. C NOTE THAT SETTING R1=0 GENERATES POINTS WITHIN A SPHERE OF C RADIUS R2. C C THE COMPUTATIONAL PROCEDURE INVOLVES FIRST GENERATING POINTS C UNIFORMLY ON THE SURFACE OF AN L-DIMENSIONAL UNIT SPHERE BY C THE METHOD DESCRIBED IN THE REFERENCE BELOW. THEN THE POINTS C ARE SCALED BY THE RANDOM VARIABLE R SO THAT THEIR LENGTHS ARE C IN THE INTERVAL [R1,R2]. WHEN R1=0 AND R2=1, R IS DISTRIBUTED C AS THE LTH ROOT OF A UNIFORM(0,1) RANDOM VARIABLE. C C SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES C RDNOR (STSPAC) - N(0,1) RANDOM DEVIATES C C CURRENT VERSION COMPLETED FEBRUARY 9, 1989 C C REFERENCE: MULLER, M.E., 'A NOTE ON A METHOD FOR GENERATING POINTS C UNIFORMLY ON N-DIMENSIONAL SPHERES', COMMUNICATIONS OF C THE ACM, 2, PP. 19-20, 1959. C----------------------------------------------------------------------- C DEFINITIONS OF PASSED PARAMETERS: C C P(LDP,*) = MATRIX (SIZE L BY N) OF POINT COORDINATES. ON OUTPUT C THE JTH COLUMN CONTAINS THE L COORDINATES OF THE JTH C POINT (J=1,2,...,N) [REAL] C C * LDP = LEADING DIMENSION OF MATRIX P AS DEFINED IN THE CALLING C PROGRAM (LDP>=L) [INTEGER] C C * L = DIMENSION OF THE POINTS AND SPHERES (L<=LDP) [INTEGER] C C * N = NUMBER OF RANDOM POINTS TO BE GENERATED. THE SECOND C DIMENSION OF P MUST BE >=N IN CALLING PROGRAM [INTEGER] C C * R1 = THE RADIUS OF THE FIRST SPHERE (R1>=0) [REAL] C C * R2 = THE RADIUS OF THE SECOND SPHERE (R2>=0) [REAL] C C IFLAG = ERROR INDICATOR [INTEGER] ON OUTPUT: C 0 -> NO ERRORS DETECTED, N POINTS GENERATED. C 1 -> EITHER L OR N IS NONPOSITIVE. C 2 -> L>LDP. C 3 -> EITHER R1<0 OR R2<0. C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION P(LDP,*) C C--- CHECK FOR ILLEGAL VALUES OF INPUT PARAMETERS AND SET ERROR FLAG C IF (MIN0(L,N).LT.1) THEN IFLAG = 1 ELSEIF (L.GT.LDP) THEN IFLAG = 2 ELSEIF (AMIN1(R1,R2).LT.0.0) THEN IFLAG = 3 ELSE IFLAG = 0 C C--- COMPUTE THREE CONSTANTS C CON1 = R1**L CON2 = R2**L POWER = 1.0/REAL(L) C C--- BEGIN LOOP FOR POINTS C DO 30 J = 1, N C C--- COMPUTE ONE UNIFORM DEVIATE (U) FOR JTH POINT, THEN COMPUTE A C--- DEVIATE R IN THE INTERVAL [R1,R2] C U = RDUNI(0) R = (U*CON2+(1.0-U)*CON1)**POWER C C--- BEGIN LOOP FOR COORDINATES OF JTH POINT C S = 0.0 DO 10 I = 1, L C C--- COMPUTE ONE NORMAL DEVIATE FOR ITH COORDINATE OF JTH POINT. S IS C--- THE SQUARED DISTANCE OF THE JTH POINT FROM THE ORIGIN. C Z = RDNOR(0) S = S+Z*Z P(I,J) = Z 10 CONTINUE C C--- COMPUTE 'COMBINED SCALE FACTOR', THEN SCALE EACH COORDINATE OF THE C--- JTH POINT SO THAT THE POINT LIES IN THE REGION BETWEEN SPHERES OF C--- RADIUS R1 AND R2. C SCALE = R/SQRT(S) DO 20 I = 1, L P(I,J) = SCALE*P(I,J) 20 CONTINUE 30 CONTINUE ENDIF RETURN C END *RDTORS SUBROUTINE RDTORS (P,LDP,N,R1,R2,IFLAG) C C----------------------------------------------------------------------- C RDTORS WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS ON THE C SURFACE OF THE TORUS WHOSE MAJOR RADIUS IS R1 AND MINOR RADIUS C IS R2. NO POINTS WILL BE GENERATED UNLESS R1>0, R2>0, R1>=R2, C AND N>0. THE TORUS LIES IN THE (X,Y) PLANE CENTERED AT THE C ORIGIN. ON OUTPUT, EACH OF THE N COLUMNS OF MATRIX P CONTAINS C THE X, Y, AND Z COORDINATES OF THE GENERATED POINTS. THE C LEADING DIMENSION OF P (LDP) MUST BE 3 OR MORE IN THE CALLING C PROGRAM. C C THE ANGLE THETA (ASSOCIATED WITH R2) IS GENERATED BY A C REJECTION METHOD REQUIRING AN AVERAGE OF 2+2*(R2/R1) UNIFORM C RANDOM DEVIATES. THE ANGLE PHI (ASSOCIATED WITH R1) REQUIRES C A SINGLE UNIFORM RANDOM DEVIATE. THEREFORE, EACH GENERATED C POINT REQUIRES AN AVERAGE OF 3+2*(R2/R1) UNIFORM DEVIATES. C C SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES C C CURRENT VERSION COMPLETED MARCH 23, 1989 C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C P(LDP,*) = MATRIX (SIZE 3 BY N) WHOSE COLUMNS ARE THE N RANDOM C POINTS GENERATED ON THE TORUS SURFACE. THE LEADING C DIMENSION OF P IS LDP (>=3) AND THE SECOND DIMENSION, C DEFINED IN THE CALLING PROGRAM, MUST BE AT LEAST N. C [REAL] C C * LDP = LEADING DIMENSION OF MATRIX P (>=3) [INTEGER] C C * N = NUMBER OF POINTS TO BE GENERATED [INTEGER] C C * R1 = MAJOR RADIUS OF THE TORUS (>0) [REAL] C C * R2 = MINOR RADIUS OF THE TORUS (>0) [REAL] C C IFLAG = ERROR FLAG ON OUTPUT. INTERPRETATION: C 0 -> NO ERRORS DETECTED. C 1 -> R1<=0 OR R2<=0 OR R1 N<1. C 3 -> LDP<3. C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION P(LDP,*) IF (AMIN1(R1,R2).LE.0.0.OR.R1.LT.R2) THEN IFLAG = 1 ELSEIF (N.LT.1) THEN IFLAG = 2 ELSEIF (LDP.LT.3) THEN IFLAG = 3 ELSE IFLAG = 0 TWOPI = 8.0*ATAN(1.0) RSUM = R1+R2 DO 20 I = 1, N 10 THETA = TWOPI*RDUNI(0) U = RSUM*RDUNI(0) R = R1+R2*COS(THETA) IF (U.GT.R) GO TO 10 PHI = TWOPI*RDUNI(0) P(1,I) = R*COS(PHI) P(2,I) = R*SIN(PHI) P(3,I) = R2*SIN(THETA) 20 CONTINUE ENDIF RETURN C END *RDTORV SUBROUTINE RDTORV (P,LDP,N,R1,R2,IFLAG) C C----------------------------------------------------------------------- C RDTORV WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: GENERATING UNIFORMLY-SPACED PSEUDO-RANDOM POINTS WITHIN THE C TORUS WHOSE MAJOR RADIUS IS R1 AND MINOR RADIUS IS R2. NO C POINTS WILL BE GENERATED UNLESS R1>0, R2>=0, R1>=R2 AND N>0. C THE TORUS LIES IN THE (X,Y) PLANE CENTERED AT THE ORIGIN. C ON OUTPUT, EACH OF THE N COLUMNS OF MATRIX P CONTAINS THE C X, Y, AND Z COORDINATES OF THE GENERATED POINTS. THE LEADING C DIMENSION OF P (LDP) MUST BE 3 OR MORE IN THE CALLING PROGRAM. C C A REJECTION PROCEDURE IS USED WHEREBY POINTS ARE FIRST C GENERATED UNIFORMLY INSIDE A TORUS-LIKE FIGURE WHOSE C CROSS-SECTIONAL AREA IS SQUARE (EACH SIDE LENGTH 2*R2). C POINTS WHICH LIE OUTSIDE THE INSCRIBED TORUS ARE REJECTED. C REGARDLESS OF R1 AND R2, THE PROPORTION OF POINTS REJECTED C IS 0.2146 = 1 - PI/4. THEREFORE, THE EXPECTED NUMER OF C UNIFORM DEVIATES REQUIRED TO GENERATE EACH POINT IS 3.5465. C C SUBPROGRAMS CALLED: RDUNI (STSPAC) - U(0,1) RANDOM DEVIATES C C CURRENT VERSION COMPLETED MARCH 8, 1989 C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C P(LDP,*) = MATRIX (SIZE 3 BY N) WHOSE COLUMNS ARE THE N RANDOM C POINTS GENERATED INSIDE THE TORUS. THE LEADING C DIMENSION OF P IS LDP (>=3) AND THE SECOND DIMENSION, C DEFINED IN THE CALLING PROGRAM, MUST BE AT LEAST N. C [REAL] C C * LDP = LEADING DIMENSION OF MATRIX P (>=3) [INTEGER] C C * N = NUMBER OF POINTS TO BE GENERATED [INTEGER] C C * R1 = MAJOR RADIUS OF THE TORUS (>0) [REAL] C C * R2 = MINOR RADIUS OF THE TORUS (>0) [REAL] C C IFLAG = ERROR FLAG ON OUTPUT. INTERPRETATION: C 0 -> NO ERRORS DETECTED. C 1 -> R1<=0 OR R2<0 OR R1 N<1. C 3 -> LDP<3. C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION P(LDP,*) IF (R1.LE.0.0.OR.R2.LT.0.0.OR.R1.LT.R2) THEN IFLAG = 1 ELSEIF (N.LT.1) THEN IFLAG = 2 ELSEIF (LDP.LT.3) THEN IFLAG = 3 ELSE IFLAG = 0 SUM = R1+R2 DIFF = R1-R2 R2SQ = R2**2 TWOPI = 8.0*ATAN(1.0) DO 20 I = 1, N 10 U = RDUNI(0) R = SQRT(U*SUM**2+(1.0-U)*DIFF**2) Z = R2*(2.0*RDUNI(0)-1.0) IF ((R-R1)**2+Z**2.GT.R2SQ) GO TO 10 PHI = TWOPI*RDUNI(0) P(1,I) = R*COS(PHI) P(2,I) = R*SIN(PHI) P(3,I) = Z 20 CONTINUE ENDIF RETURN C END *REJ1 SUBROUTINE REJ1 (Y,N,NPASS,NREJ,SMEAN,SSD,IFLAG) C C----------------------------------------------------------------------- C REJ1 WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MARYLAND 20899 C C FOR: COMPUTING THE MEAN AND STANDARD DEVIATION OF A SAMPLE OF C 'NORMAL' DATA IN WHICH OUTLIERS MAY BE PRESENT. OUTLIERS ARE C FIRST REJECTED BY A PROCEDURE BASED ON THE SHORTEST INTERVAL C COVERING HALF THE POINTS. THE PROCEDURE IS ITERATED UNTIL C C 1) A USER-SPECIFIED NUMBER OF PASSES OCCURS, OR C 2) THE PROPORTION OF VALUES REJECTED IN A GIVEN PASS IS C 0.01 OR LESS. C C SIMULATION STUDIES ON NORMAL DATA WERE USED TO DETERMINE C THE APPROPRIATE VALUES OF CONSTANTS IN THIS PROGRAM. THEY C WERE CHOSEN SO THAT, ON THE FIRST PASS, THE EXPECTED PROPOR- C TION OF 'GOOD' VALUES REJECTED WAS 0.01 REGARDLESS OF SAMPLE C SIZE. WHEN THE NUMBER OF PASSES ARE NOT LIMITED, THE ACTUAL C PROPORTION OF VALUES REJECTED WAS FOUND TO BE 0.010 TO 0.012 C FOR ALL SAMPLE SIZES. C C THE PROCEDURE WAS ORIGINALLY DESIGNED FOR USE ON LARGE SETS C OF 'NORMAL' DATA ASYMMETRICALLY CONTAMINATED WITH UP TO 50% C OUTLIERS. ITS BEHAVIOR WAS EXAMINED FOR SAMPLE SIZES OF 15 C TO 10,000 AND IT APPEARS TO WORK WELL. WHEN THE SAMPLE SIZE C IS 25 OR LESS, HOWEVER, THE USER MAY WANT TO CONSIDER THE C WELL-ESTABLISHED DIXON TEST AS AN ALTERNATIVE. THAT TEST IS C DISCUSSED IN MANY STATISTICS BOOKS. C C NOTE: ON INPUT, THE SAMPLE SIZE MUST BE AT LEAST 15 AND THE C USER-SPECIFIED MAXIMUM NUMBER OF PASSES MUST BE AT LEAST ONE. C ON RETURN, THE NUMBER OF 'GOOD' VALUES CAN BE COMPUTED AS: C C N - NREJ(1) - NREJ(2) - ... - NREJ(NPASS). C C SUBPROGRAMS CALLED: SORT1 (STSPAC) C C CURRENT VERSION COMPLETED OCTOBER 2, 1989 C C REFERENCE: REEVE, CHARLES P., 'AN OUTLIER REJECTION PROCEDURE, C BASED ON RANGES, FOR ASYMMETRICALLY-CONTAMINATED NORMAL C DATA', STATISTICAL ENGINEERING DIVISION NOTE 89-2, C DECEMBER 1989. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * Y(*) = VECTOR (LENGTH N) OF SAMPLE VALUES. MUST BE C DIMENSIONED AT LEAST N IN CALLING PROGRAM. [REAL] C C * N = SAMPLE SIZE (N>=15) [INTEGER] C C * NPASS = ON INPUT, MAXIMUM NUMBER OF PASSES (ITERATIONS) OF C THE OUTLIER REJECTION PROCEDURE. ON OUTPUT, THE C ACTUAL NUMBER OF PASSES. [INTEGER] C C NREJ(*) = VECTOR (LENGTH NPASS) CONTAINING, ON OUTPUT, THE C NUMBER OF VALUES REJECTED AT EACH PASS. IT MUST C BE DIMENSIONED AT LEAST NPASS IN THE CALLING C PROGRAM. [INTEGER] C C SMEAN = THE COMPUTED SAMPLE MEAN OF THE NON-REJECTED VALUES C [REAL] C C SSD = THE COMPUTED SAMPLE STANDARD DEVIATION OF THE C NON-REJECTED VALUES [REAL] C C IFLAG = ERROR INDICATOR ON OUTPUT [INTEGER] INTERPRETATION: C 0 -> NO ERRORS DETECTED. C 1 -> N<15 ON INPUT. C 2 -> NPASS<1 ON INPUT. C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION NREJ(*),Y(*) LOGICAL LG DATA CNORM / 1.349 / DATA C1,C2,C3,C4 / 2.576,9.573,-3.013,-0.6989 / DATA C5,C6,C7,C8 / 2.576,7.889,1.687,-0.6729 / C C--- CHECK FOR VALIDITY OF INPUT VALUES C IF (N.LT.15) THEN IFLAG = 1 ELSEIF (NPASS.LT.1) THEN IFLAG = 2 ELSE IFLAG = 0 C C--- SORT Y-VALUES C CALL SORT1 (Y,1,N) C C--- DETERMINE OUTLIERS BY FINDING THE SHORTEST INTERVAL COVERING C--- HALF THE POINTS C NADJ = N NIT = 0 10 NIT = NIT+1 RNA = REAL(NADJ) IF (MOD(NADJ,2).EQ.0) THEN SIGMLT = C1+C2*(RNA+C3)**C4 ELSE SIGMLT = C5+C6*(RNA+C7)**C8 ENDIF L = (NADJ+1)/2 RMIN = Y(N)-Y(1) DO 20 I = 1, N-L R = Y(I+L)-Y(I) IF (R.LE.RMIN) THEN RMIN = R K = I ENDIF 20 CONTINUE SMEAN = 0.5*(Y(K)+Y(K+L)) BOUND = SIGMLT*RMIN/CNORM C C--- TRIM OUTLIERS AT LOWER END C NLO = 1 30 IF (SMEAN-Y(NLO).GT.BOUND) THEN NLO = NLO+1 GO TO 30 C ENDIF C C--- TRIM OUTLIERS AT UPPER END C NHI = N 40 IF (Y(NHI)-SMEAN.GT.BOUND) THEN NHI = NHI-1 GO TO 40 C ENDIF NGOOD = NHI-NLO+1 NREJ(NIT) = NADJ-NGOOD LG = NIT.EQ.NPASS.OR.NGOOD.LT.15 IF (NREJ(NIT).LE.1+INT(0.01*RNA).OR.LG) THEN NPASS = NIT C C--- COMPUTE MEAN AND STANDARD DEVIATION OF NON-REJECTED VALUES C SMEAN = 0.0 DO 50 I = NLO, NHI SMEAN = SMEAN+Y(I) 50 CONTINUE SMEAN = SMEAN/REAL(NGOOD) SSD = 0.0 DO 60 I = NLO, NHI SSD = SSD+(Y(I)-SMEAN)**2 60 CONTINUE SSD = SQRT(SSD/REAL(NGOOD-1)) ELSE NADJ = NGOOD GO TO 10 C ENDIF ENDIF RETURN C END *SKEKUR SUBROUTINE SKEKUR (Y,NLO,NHI,YSKEW,YKURT,IOPT) C C----------------------------------------------------------------------- C SKEKUR WRITTEN BY CHARLES P. REEVE C C FOR: COMPUTING SKEWNESS AND KURTOSIS FOR ENTRIES NLO THROUGH NHI C IN VECTOR Y. THE VALUES MAY BE CENTERED ABOUT EITHER THE C MEAN (IOPT <> 0) OR ABOUT ZERO (IOPT = 0). THE TRADITIONAL C DIVISIOR OF N (NOT N-1) IS USED WHEN THE MEAN IS ESTIMATED. C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED FEBRUARY 28, 1986 C----------------------------------------------------------------------- C DIMENSION Y(*) RN = REAL(NHI-NLO+1) IF (IOPT.EQ.0) THEN S = 0. ELSE S = 0. DO 10 I = NLO, NHI S = S+Y(I) 10 CONTINUE S = S/RN ENDIF T2 = 0. T3 = 0. T4 = 0. DO 20 I = NLO, NHI D = Y(I)-S T2 = T2+D**2 T3 = T3+D**3 T4 = T4+D**4 20 CONTINUE YSKEW=SQRT(RN)*T3/T2**1.5 YKURT=RN*T4/T2**2 RETURN C END *SORT1 SUBROUTINE SORT1 (X,M,N) C C----------------------------------------------------------------------- C SORT1 OBTAINED BY CHARLES P. REEVE FROM DR. D.A. ZAHN AT THE C FLORIDA STATE UNIVERSITY, TALLAHASSEE, FLORIDA UNDER THE C NAME 'FTASORT'. SEVERAL MINOR CORRECTIONS WERE MADE TO C THE ORIGINAL VERSION AFTER LINE NUMBER 90. C C FOR: SORTING THE SEGMENT OF A REAL ARRAY BETWEEN ENTRIES M AND N C FROM SMALLEST TO LARGEST. THE ARRAYS LP(K) AND LQ(K) PERMIT C SORTING UP TO 2**(K+1)-1 ELEMENTS, THUS K=25 PERMITS SORTING C UP TO 67,108,863 ELEMENTS. C C SUBROUTINES CALLED: -NONE- C C CURRENT VERSION COMPLETED JANUARY 30, 1978 C----------------------------------------------------------------------- C DIMENSION LP(25),LQ(25),X(*) IF (M.GE.N) RETURN I = 1 KM = M KN = N IF (KM.GE.KN) GO TO 70 10 J = KM K = (KN+KM)/2 A = X(K) IF (X(KM).LE.A) GO TO 20 X(K) = X(KM) X(KM) = A A = X(K) 20 L = KN IF (X(KN).GE.A) GO TO 40 X(K) = X(KN) X(KN) = A A = X(K) IF (X(KM).LE.A) GO TO 40 X(K) = X(KM) X(KM) = A A = X(K) GO TO 40 30 X(L) = X(J) X(J) = B 40 L = L-1 IF (X(L).GT.A) GO TO 40 B = X(L) 50 J = J+1 IF (X(J).LT.A) GO TO 50 IF (J.LE.L) GO TO 30 IF (L-KM.LE.KN-J) GO TO 60 LP(I) = KM LQ(I) = L KM = J I = I+1 GO TO 80 60 LP(I) = J LQ(I) = KN KN = L I = I+1 GO TO 80 70 I = I-1 IF (I.EQ.0) RETURN KM = LP(I) KN = LQ(I) 80 IF (KN-KM.GE.11) GO TO 10 KM = KM-1 90 KM = KM+1 IF (KM.EQ.KN) GO TO 70 A = X(KM+1) IF (X(KM).LE.A) GO TO 90 J = KM 100 X(J+1) = X(J) J = J-1 IF (J.EQ.M-1) GO TO 110 IF (A.LT.X(J)) GO TO 100 110 X(J+1) = A GO TO 90 END *SORT2 SUBROUTINE SORT2 (X,Y,M,N) C C----------------------------------------------------------------------- C SORT2 OBTAINED BY CHARLES P. REEVE FROM DR. D. A. ZAHN AT C THE FLORIDA STATE UNIVERSITY, TALLAHASSEE, FLORIDA C UNDER THE NAME *FTASORT*. SEVERAL MINOR CORRECTIONS C WERE MADE TO THE ORIGINAL VERSION AFTER LINE NUMBER 90. C C FOR: SORTING THE SEGMENT OF REAL ARRAY X BETWEEN ENTRIES M AND C N FROM SMALLEST TO LARGEST. A SECOND ARRAY Y IS CARRIED C ALONG AND SORTED IDENTICALLY AS THE FIRST. THIS FEATURE C IS A MODIFICATION TO THE SUBROUTINE SORT1. C C SUBROUTINES CALLED: -NONE- C C CURRENT VERSION COMPLETED OCTOBER 30, 1978 C C NOTE: ARRAYS LP(K) AND LQ(K) PERMIT SORTING UP TO 2**(K+1)-1 C ELEMENTS, I.E., FOR K=25 YOU MAY SORT 67,108,863 ELEMENTS. C----------------------------------------------------------------------- C DIMENSION LP(25),LQ(25),X(1),Y(1) I = 1 KM = M KN = N IF (KM.GE.KN) GO TO 70 10 J = KM K = (KN+KM)/2 A = X(K) B = Y(K) IF (X(KM).LE.A) GO TO 20 X(K) = X(KM) Y(K) = Y(KM) X(KM) = A Y(KM) = B A = X(K) B = Y(K) 20 L = KN IF (X(KN).GE.A) GO TO 40 X(K) = X(KN) Y(K) = Y(KN) X(KN) = A Y(KN) = B A = X(K) B = Y(K) IF (X(KM).LE.A) GO TO 40 X(K) = X(KM) Y(K) = Y(KM) X(KM) = A Y(KM) = B A = X(K) B = Y(K) GO TO 40 C 30 X(L) = X(J) Y(L) = Y(J) X(J) = G Y(J) = H 40 L = L-1 IF (X(L).GT.A) GO TO 40 G = X(L) H = Y(L) 50 J = J+1 IF (X(J).LT.A) GO TO 50 IF (J.LE.L) GO TO 30 IF (L-KM.LE.KN-J) GO TO 60 LP(I) = KM LQ(I) = L KM = J I = I+1 GO TO 80 C 60 LP(I) = J LQ(I) = KN KN = L I = I+1 GO TO 80 C 70 I = I-1 IF (I.EQ.0) RETURN KM = LP(I) KN = LQ(I) 80 IF (KN-KM.GE.11) GO TO 10 KM = KM-1 90 KM = KM+1 IF (KM.EQ.KN) GO TO 70 A = X(KM+1) B = Y(KM+1) IF (X(KM).LE.A) GO TO 90 J = KM 100 X(J+1) = X(J) Y(J+1) = Y(J) J = J-1 IF (J.EQ.M-1) GO TO 110 IF (A.LT.X(J)) GO TO 100 110 X(J+1) = A Y(J+1) = B GO TO 90 C END C *TRISYS SUBROUTINE TRISYS (SUB,DIAG,SUP,B,N,IFLAG) C C----------------------------------------------------------------------- C TRISYS WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL INSTITUTE FOR STANDARDS AND TECHNOLOGY, C GAITHERSBURG, MD 20899 C C FOR: SOLVING THE TRIDIAGONAL SYSTEM OF LINEAR EQUATIONS C C SUB X + DIAG X + SUP X = B C I I-1 I I I I+1 I C C WHERE I = 1,2,...,N. SUB(1) AND SUP(N) ARE ZERO BY DEFINITION, C THEREFORE THEIR PASSED VALUES ARE IGNORED. THE CONTENTS OF C DIAG ARE ALTERED AND THE SOLUTION VECTOR (X) IS RETURNED IN B. C SEE ALGORITHM 3.3 IN THE REFERENCE. C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED SEPTEMBER 23, 1988 C C REFERENCE: CONTE, S.D. AND DEBOOR, CARL, 'ELEMENTARY NUMERICAL C ANALYSIS', MCGRAW-HILL BOOK COMPANY, 1972, PP. 121-122. C----------------------------------------------------------------------- C DEFINITION OF PASSED PARAMETERS: C C * SUB(*) = VECTOR (LENGTH N) OF SUB-DIAGONAL ELEMENTS. THE VALUE C IN SUB(1) IS NOT USED. [REAL] C C * DIAG(*) = VECTOR (LENGTH N) OF DIAGONAL ELEMENTS. ALL VALUES C ARE CHANGED ON OUTPUT. [REAL] C C * SUP(*) = VECTOR (LENGTH N) OF SUPER-DIAGONAL ELEMENTS. THE C VALUE IN SUP(N) IS NOT USED. [REAL] C C * B(*) = VECTOR (LENGTH N) CONTAINING THE RIGHT HAND SIDE OF C THE SYSTEM ON INPUT, AND THE SOLUTION VECTOR ON OUTPUT. C [REAL] C C * N = THE NUMBER OF EQUATIONS AND UNKNOWNS [INTEGER] C C IFLAG = ERROR INDICATOR ON OUTPUT. INTERPRETATION: C -1 -> N<1, THEREFORE NO SOLUTION RETURNED. C 0 -> NO ERRORS DETECTED. C K -> AT SOME POINT DIAG(K)=0.0, THEREFORE FAILURE. C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C DIMENSION SUB(*),DIAG(*),SUP(*),B(*) IFLAG = 0 C C--- CHECK FOR N<1 C IF (N.LT.1) THEN IFLAG = -1 RETURN C ENDIF C C--- SOLVE TRIDIAGONAL SYSTEM C DO 10 K = 2, N IF (DIAG(K-1).EQ.0.0) THEN IFLAG = K-1 RETURN C ELSE RATIO = -SUB(K)/DIAG(K-1) DIAG(K) = DIAG(K)+RATIO*SUP(K-1) B(K) = B(K)+RATIO*B(K-1) ENDIF 10 CONTINUE IF (DIAG(N).EQ.0.0) THEN IFLAG = N RETURN C ELSE B(N) = B(N)/DIAG(N) ENDIF DO 20 K = N-1, 1, -1 IF (DIAG(K).EQ.0.0) THEN IFLAG = K RETURN C ELSE B(K) = (B(K)-SUP(K)*B(K+1))/DIAG(K) ENDIF 20 CONTINUE RETURN C END *UPDATE SUBROUTINE UPDATE (OLDMN,OLDSD,NEWMN,NEWSD,YNEW,N) C C----------------------------------------------------------------------- C UPDATE WRITTEN BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 C C FOR: TAKING A PREVIOUSLY COMPUTED MEAN (OLDMN) AND STANDARD C DEVIATION (OLDSD) BASED ON N-1 VALUES AND COMPUTING AN C UPDATED MEAN (NEWMN) AND STANDARD DEVIATION (OLDSD) BASED C ON THE N-1 PREVIOUS VALUES PLUS THE NTH VALUE (YNEW). C C NOTE: ALL PASSED PARAMETERS OTHER THAN N ARE OF TYPE REAL. C C SPECIAL CASES: C C IF N < 1 THEN A MESSAGE IS PRINTED AND EXECUTION STOPPED. C IF N = 1 THEN INPUT VALUES OF OLDMN AND OLDSD ARE IGNORED C WHILE COMPUTING NEWMN = YNEW AND NEWSD = 0.0 . C C SUBPROGRAMS CALLED: -NONE- C C CURRENT VERSION COMPLETED MARCH 2, 1986 C----------------------------------------------------------------------- C REAL NEWMN,NEWSD IF(N.GE.2) THEN RN = REAL(N) RN1 = REAL(N-1) NEWMN = (RN1*OLDMN+YNEW)/RN NEWSD = SQRT(REAL(N-2)*OLDSD**2/RN1+(YNEW-OLDMN)**2/RN) ELSE IF(N.EQ.1) THEN NEWMN=YNEW NEWSD=0.0 ELSE PRINT *,' *** N MUST BE A POSITIVE INTEGER' PRINT *,' *** EXECUTION STOPPED IN SUBROUTINE UPDATE' STOP END IF RETURN END *ZEROBR SUBROUTINE ZEROBR (FUNC,X1,X2,TOL,X0,IFLAG) C C----------------------------------------------------------------------- C ZEROBR COPIED BY CHARLES P. REEVE, STATISTICAL ENGINEERING C DIVISION, NATIONAL BUREAU OF STANDARDS, GAITHERSBURG, C MARYLAND 20899 FROM REFERENCE 1 BELOW WITH SLIGHT C MODIFICATIONS. C C FOR: FINDING THE REAL ZERO X0 OF THE REAL FUNCTION FUNC KNOWN TO C LIE BETWEEN X1 AND X2 USING BRENT'S METHOD AS DESCRIBED IN C THE REFERENCES. THE ZERO IS REFINED UNTIL ITS ABSOLUTE C ACCURACY IS TOL. C C NOTE: THE CONSTANT EPS FOR MACHINE FLOATING POINT PRECISION IS C MACHINE DEPENDENT. C C SUBPROGRAMS CALLED: FUNC (USER DEFINED) 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 FUNC = AN EXTERNAL FUNCTION WITH ONE ARGUMENT (REAL) C C * X1 = THE LEFT ENDPOINT OF AN INTERVAL CONTAINING A ZERO (REAL) C C * X2 = THE RIGHT ENDPOINT OF AN INTERVAL CONTAINING A ZERO (REAL) C C * TOL = THE ABSOLUTE TOLERANCE ON THE COMPUTED ZERO (REAL) C C X0 = THE COMPUTED ZERO OF FUNC IN THE INTERVAL (X1,X2) (REAL) C C IFLAG = ERROR INDICATOR ON OUTPUT (INTEGER) INTERPRETATION: C 0 -> NO ERRORS DETECTED C 1 -> FUNC(X1) AND FUNC(X2) HAVE THE SAME SIGN - THE ZERO C IS NOT DEFINED C 2 -> MAXIMUM ITERATIONS EXCEEDED - CURRENT VALUE OF ZERO C RETURNED C C * INDICATES PARAMETERS REQUIRING INPUT VALUES C----------------------------------------------------------------------- C C--- DEFINE MAXIMUM ITERATIONS AND MACHINE FLOATING POINT PRECISION C DATA ITMAX,EPS / 50,1.0E-14 / A = X1 B = X2 FA = FUNC(A) FB = FUNC(B) C C--- CHECK FOR ZERO BEING BRACKETED C IF (FB*FA.GT.0.0) THEN IFLAG = 1 RETURN C ENDIF FC = FB DO 10 ITER = 1, ITMAX IF (FB*FC.GT.0.0) THEN C C--- RENAME A,B,C AND ADJUST BOUNDING INTERVAL D C C = A FC = FA D = B-A E = D ENDIF IF (ABS(FC).LT.ABS(FB)) THEN A = B B = C C = A FA = FB FB = FC FC = FA ENDIF C C--- CONVERGENCE CHECK C TOL1 = 2.0*EPS*ABS(B)+0.5*TOL XM = 0.5*(C-B) IF (ABS(XM).LE.TOL1.OR.FB.EQ.0.0) THEN X0 = B RETURN C ENDIF IF (ABS(E).GE.TOL1.AND.ABS(FA).GT.ABS(FB)) THEN C C--- ATTEMPT INVERSE QUADRATIC INTERPOLATION C S = FB/FA IF (A.EQ.C) THEN P = 2.0*XM*S Q = 1.0-S ELSE Q = FA/FC R = FB/FC P = S*(2.0*XM*Q*(Q-R)-(B-A)*(R-1.0)) Q = (Q-1.0)*(R-1.0)*(S-1.0) ENDIF C C--- CHECK WHETHER IN BOUNDS C IF (P.GT.0.0) Q = -Q P = ABS(P) IF (2.0*P.LT.AMIN1(3.0*XM*Q-ABS(TOL1*Q),ABS(E*Q))) THEN C C--- ACCEPT INTERPLATION C E = D D = P/Q ELSE C C--- INTERPOLATION FAILED, USE BISECTION C D = XM E = D ENDIF ELSE C C--- BOUNDS DECREASING TOO SLOWLY, USE BISECTION C D = XM E = D ENDIF C C--- MOVE LAST BEST GUESS TO A C A = B FA = FB C C--- EVALUATE NEW TRIAL ZERO C IF (ABS(D).GT.TOL1) THEN B = B+D ELSE B = B+SIGN(TOL1,XM) ENDIF FB = FUNC(B) 10 CONTINUE IFLAG = 2 X0 = B RETURN C END