FUNCTION PLN(N, S2, B2, S1, B1) PLN 10 C COMPUTES PROB FOR BOUNDARIES OF FORM S*Y+B. N IS C SAMPLE SIZE, S2,B2 ARE PARAMETERS FOR UPPER BOUNDARY, C CALLS SUBROUTINE RAKK. TO CALL DURB C CHANGE AS INDICATED BELOW. DIMENSION YY(2,202), S(2), B(2) C MAX SIZE FOR N IS 200 S(1) = S1 S(2) = S2 B(1) = B1 B(2) = B2 IF (N.GT.0) GO TO 10 PLN = 0. RETURN 10 IF (S(1).GE.0.) GO TO 20 S(1) = 0. 20 IF (S(2).GE.0.) GO TO 30 B(2) = S(2) + B(2) S(2) = 0. 30 N1 = 1 + N AN = N DO 50 I=1,2 DO 40 K=1,N1 C = FLOAT(K-1)/AN - B(I) YY(I,K) = -.5 IF (C.LE.0.) GO TO 40 YY(I,K) = 1.1 IF (S(I).GT.C) YY(I,K) = C/S(I) 40 CONTINUE 50 CONTINUE CALL RAKK(N, YY, ANS) C TO CALL DURB CHANGE PRECEDING STATEMENT TO C CALL DURB (N,YY,ANS) PLN = ANS RETURN END SUBROUTINE RAKK(N, YY, ANS) RAK 10 C N IS SAMPLE SIZE. YY(2,.),YY(1,.) C CONTAIN, RESPECTIVELY, UPPER AND LOWER BOUNDARY C LEVEL POINTS CORRESPONDING TO 0,1/N,2/N,...,1. MUST C BE N+1 OF EACH. ROUTINE ASSUMES THEY ARE NONDECREASING. DIMENSION W(201,2), YY(2,202), KBLOC(41), FRNG(41,2), SSS(41) DATA RNG, RRNG, SQRNG, SQRRNG, MAXNB /1.E292,1.E-292,1.E146, * 1.E-146,40/ C RNG AND RRNG ARE MACHINE DEPENDENT CONSTANTS SELECTED C SO RRNG=1./RNG G.E. SMALLEST POSITIVE FLOATING C NUMBER AND RNG L.E. LARGEST POSITIVE FLOATING C NUMBER OF MACHINE. SQRNG AND SQRRNG ARE SQUARE ROOTS C OF RNG AND RRNG. C MAXNB IS THE MAX NUMBER OF BLOCKS INTO WHICH C W(.,.) WILL BE PARTITIONED. IF MAXNB IS C INCREASED REVISE DIMENSION STATEMENT ACCORDINGLY. N1 = N + 1 TF = 1. NNN = 1 IEXT = 1 M2 = 1 M1 = 2 W(1,M2) = SQRNG NBLOC = 1 KBLOC(1) = 1 KBLOC(2) = 2 FRNG(1,M2) = 1. ANS = 0. IF (N1.LE.1) RETURN YY(2,N1+1) = 1. Y1 = 0. IF ((YY(1,1).LE.0.) .OR. (YY(2,N1).GE.1.) .OR. * (YY(1,N1).LE.1.) .OR. (YY(2,1).GE.0.)) RETURN JBOTP = 1 JTOPP = 1 JTOP1 = 1 C ENTRY POINT FOR SUBSEQUENT ITERATIONS. 10 JBOT0 = JBOTP JTOP0 = JTOP1 MT = M1 M1 = M2 M2 = MT Y0 = Y1 C COMPUTE SUMMATION INDICES, JBOTP,JTOPP FOR NEXT C ITERATION. 20 IF (YY(2,JTOPP+1).GT.Y0) GO TO 30 JTOPP = JTOPP + 1 GO TO 20 30 Y1 = AMIN1(YY(2,JTOPP+1),YY(1,JBOTP)) IF (Y1.GE.1.) GO TO 50 40 IF (YY(1,JBOTP).GT.Y1) GO TO 60 JBOTP = JBOTP + 1 GO TO 40 50 Y1 = 1. C SET EXIT FLAG. NEXT ITERATION IS FINAL ONE. IEXT = 2 JBOTP = N1 60 IF (JBOTP.GT.JTOPP) RETURN C RETURN WITH ANS=0 SINCE NO PATHS BETWEEN BOUNDARIES. JTOP1 = JTOPP P = Y1 - Y0 KBLOC(NBLOC+1) = JTOP1 + 1 DO 70 L=JBOT0,JTOP1 W(L,M2) = 0. 70 CONTINUE L0 = KBLOC(NBLOC) L1 = JTOP0 DO 80 L=L0,L1 IF (W(L,M1).GT.SQRRNG) GO TO 80 JTOP0 = L - 1 GO TO 90 80 CONTINUE 90 IF (JBOT0.LT.KBLOC(2)) GO TO 110 C THE NEXT FEW STATEMENTS ELIMINATE BLOCKS NO C LONGER USED. DO 100 I=1,NBLOC FRNG(I,M1) = FRNG(I+1,M1) KBLOC(I) = KBLOC(I+1) 100 CONTINUE NBLOC = NBLOC - 1 CALL ACCUM(RRNG/FRNG(1,M1), TF, NNN) GO TO 90 110 KBLOC(1) = JBOT0 MBAV = (MAXNB-NBLOC+1)/2 IF (P.LE.(1./FLOAT(4*N1))) MBAV = 2 MBAV = MIN0(MAXNB,NBLOC+MBAV) C MBAV IS UPPER LIMIT ON NUMBER OF BLOCKS DURING NEXT C ITERATION. MBAV IS SET TO LIMIT NUMBER OF NEW C BLOCKS WHEN P IS SMALL. CONVI = 1. I1 = 1 J1 = 1 120 DO 400 I=I1,NBLOC KBLI = KBLOC(I) IF (I.EQ.I1) GO TO 170 FRNG(I,M2) = AMIN1(FRNG(I,M1),1./(SQRNG*W(KBLI-1,M2))) CONVI = (CONVI/FRNG(I-J1+1,M1))*FRNG(I,M2) IF (J1.EQ.1) GO TO 170 K1 = KBLOC(I) - KBLOC(I-J1+2) + 1 K0 = KBLOC(I-1) - KBLOC(I-J1+1) + 1 IF (K1-K0) 150, 170, 130 130 K0P1 = K0 + 1 DO 140 K=K0P1,K1 CONVI = CONVI*P/FLOAT(K) 140 CONTINUE GO TO 170 150 K1P1 = K1 + 1 DO 160 K=K1P1,K0 CONVI = CONVI*FLOAT(K)/P 160 CONTINUE 170 IF (CONVI.GT.RRNG) GO TO 210 DO 190 J=J1,J11 C CHANGE UPPER LIMIT IF EARLY EXIT FROM PREVIOUS LOOP CONVI = FRNG(I,M2)*RNG K0 = MAX0(1,KBLOC(I-1)-MIN0(JTOP0+1,KBLOC(I-J+1))+2) K1 = KBLOC(I) - MIN0(JTOP0+1,KBLOC(I-J+1)) + 1 DO 180 K=K0,K1 CONVI = CONVI*P/FLOAT(K) 180 CONTINUE CONVI = SSS(J)*CONVI IF (CONVI.GT.RRNG) GO TO 200 190 CONTINUE GO TO (10, 440), IEXT 200 J1 = J + 1 210 DO 390 J=J1,I LFG = 1 C LFG=2 IF NEXT ITERATION ON J BECOMES NECESSARY. LB1 = KBLOC(I) K01 = KBLOC(I) - KBLOC(I-J+1) + 1 S = CONVI IF (J.EQ.J1) GO TO 240 IF (TMLT) 220, 400, 230 220 S = -TMLT*FRNG(I-J+2,M1) IF (S.GT.RRNG) GO TO 240 230 S = TMLT*(FRNG(I-J+2,M1)*RNG) 240 K0 = MAX0(KBLOC(I)-MIN0(KBLOC(I-J+2)-1,JTOP0)+1,1) K1 = KBLOC(I+1) - KBLOC(I-J+1) TMLT = 0. RRNGS = 1. IF (S.LT.0.) RRNGS = -RRNG SSS(J) = S DO 370 K=K0,K1 LDF = 1 - K LB = MAX0(KBLOC(I),KBLOC(I-J+1)-LDF,LB1) LT = MIN0(KBLOC(I+1)-1,MIN0(JTOP0,KBLOC(I-J+2)-1)-LDF) IF (LB.GT.LB1) LFG = 2 IF (LB.GT.LT) GO TO 270 DO 260 L=LB,LT IF (W(L,M2).EQ.((W(L+LDF,M1)*RRNGS)*S+W(L,M2))) GO TO * 250 LB1 = L C WHEN THE SUMMANDS BECOME SO SMALL THAT THE VALUE OF C W(L,M2) IS UNCHANGED LB1 IS INCREASED. THIS PREVENTS C FURTHER ADDITIONS TO SUCH W(L,M2). GO TO 290 250 IF (W(L,M2).GT.RRNG) GO TO 260 IF (K-K01) 340, 340, 280 260 CONTINUE 270 LB1 = LT + 1 IF ((LB1.NE.KBLOC(I+1)) .OR. (K.LE.K01)) GO TO 340 280 J11 = J IF (W(KBLI,M2).LE.RRNG) LFG = 2 GO TO (400, 380), LFG 290 IF (S) 300, 380, 320 300 DO 310 L=LB1,LT W(L,M2) = (W(L+LDF,M1)*RRNGS)*S + W(L,M2) 310 CONTINUE GO TO 340 320 DO 330 L=LB1,LT W(L,M2) = W(L+LDF,M1)*S + W(L,M2) 330 CONTINUE 340 T = S S = S*P/FLOAT(K) IF (S) 360, 350, 360 350 IF (T.LE.0.) GO TO 380 S = -(T*RNG)*P/FLOAT(K) RRNGS = -RRNG C THE NEGATIVE BIT IN S,TMLT IS USED TO C INDICATE ITS VALUE IS RNG TIMES THE USUAL VALUE. 360 IF (K.EQ.K01) TMLT = S C THE ITERATION ON K MUST CONTINUE UNTIL TMLT IS C SET. TMLT USED IN NEXT ITERATION ON J. 370 CONTINUE 380 J11 = J 390 CONTINUE 400 CONTINUE C CREATE NEW BLOCK IF NECESSARY. IF (NBLOC.GE.MBAV) GO TO (10, 440), IEXT L1 = KBLOC(NBLOC+1) - KBLOC(NBLOC) DO 410 L=1,L1 LL = KBLOC(NBLOC+1) - L IF (W(LL,M2).GE.SQRRNG) GO TO 420 410 CONTINUE 420 IF (LL.EQ.JTOP1) GO TO (10, 440), IEXT KBLOC(NBLOC+2) = KBLOC(NBLOC+1) KBLOC(NBLOC+1) = LL + 1 NBLOC = NBLOC + 1 L0 = KBLOC(NBLOC) DO 430 L=L0,JTOP1 W(L,M2) = 0. 430 CONTINUE I1 = NBLOC FRNG(NBLOC,M2) = 1. CONVI = 0. GO TO 120 440 CALL ACCUM(W(N1,M2), TF, NNN) CALL ACCUM(SQRRNG, TF, NNN) IF (NBLOC.EQ.1) GO TO 460 DO 450 I=2,NBLOC CALL ACCUM(RRNG/FRNG(I,M2), TF, NNN) 450 CONTINUE 460 ANS = TF IF (N-NNN) 470, 490, 500 470 DO 480 J=N1,NNN ANS = ANS/FLOAT(J) 480 CONTINUE 490 RETURN 500 NNN = NNN + 1 DO 510 J=NNN,N ANS = ANS*FLOAT(J) 510 CONTINUE RETURN END SUBROUTINE ACCUM(F, TF, NNN) ACC 10 C THIS SUBROUTINE ANCILLARY TO RAKK AND ACCUMULATES THE C MULTIPLIERS USED TO KEEP W(.,.) IN RANGE. TF = TF*F 10 IF ((TF.GT.1.) .OR. (TF.EQ.0.)) RETURN NNN = NNN + 1 TF = TF*FLOAT(NNN) GO TO 10 END SUBROUTINE DURB(N, YY, ANS) DUR 10 C DURBINS METHOD. J.APPL.PROB., 8(1971), PP431-453, C FORMULAS (36),(37),(40). C N IS SAMPLE SIZE. YY(2,.),YY(1,.) C CONTAIN, RESPECTIVELY, UPPER AND LOWER BOUNDARY C LEVEL POINTS CORRESPONDING TO 0,1/N,2/N,...,1. MUST C BE N+1 OF EACH. ROUTINE ASSUMES THEY ARE NONDECREASING. DIMENSION ALPHA(201), BETA(201), YY(2,202) C DIMENSION STATEMENT PERMITS MAX SAMPLE SIZE OF 200 N1 = N + 1 ANS = 0. IF ((YY(1,1).LE.0.) .OR. (YY(2,N1).GE.1.) .OR. * (YY(1,N1).LE.1.) .OR. (YY(2,1).GE.0.)) RETURN IP = 1 DO 20 I=1,N1 IF (YY(2,I).LT.0.) GO TO 10 IP = I GO TO 30 10 YY(2,I) = 0. 20 CONTINUE 30 DO 40 I=1,N1 IF (YY(1,I).LE.1.) GO TO 40 YY(1,I) = 1. 40 CONTINUE ALPHA(1) = 0. BETA(1) = 1. DO 120 K=2,N1 ALPHA(K) = YY(2,K)**(K-1) IF (IP.GT.(K-1)) GO TO 60 COMB = 1. I1 = K - IP DO 50 I=1,I1 COMB = COMB*FLOAT(K-I)/FLOAT(I) C COMB IS BINOMIAL COEFFICIENT K-1 OVER I-1. ALPHA(K) = ALPHA(K) - ALPHA(K-I)*(YY(2,K)-YY(2,K-I))**I* * COMB 50 CONTINUE 60 COMB = 1. DO 70 J=1,N1 IF (YY(2,K).LE.YY(1,J)) GO TO 80 ALPHA(K) = ALPHA(K) - BETA(J)*(YY(2,K)-YY(1,J))**(K-J)* * COMB C COMB IS BINOMIAL COEFFICIENT K-1 OVER J-1. COMB = COMB*FLOAT(K-J)/FLOAT(J) 70 CONTINUE 80 BETA(K) = YY(1,K)**(K-1) IF (IP.GT.K) GO TO 100 COMB = 1. I1 = K - IP + 1 DO 90 I=1,I1 BETA(K) = BETA(K) - ALPHA(K-I+1)*(YY(1,K)-YY(2,K-I+1))** * (I-1)*COMB C COMB IS BINOMIAL COEFFICIENT K-1 OVER I-1. COMB = COMB*FLOAT(K-I)/FLOAT(I) 90 CONTINUE COMB = 1. 100 KM1 = K - 1 DO 110 J=1,KM1 IF (YY(1,J).EQ.1.) GO TO 120 C COMB IS BINOMIAL COEFFICIENT K-1 OVER J-1. BETA(K) = BETA(K) - BETA(J)*(YY(1,K)-YY(1,J))**(K-J)*COMB COMB = COMB*FLOAT(K-J)/FLOAT(J) 110 CONTINUE 120 CONTINUE COMB = 1. I1 = N1 - IP + 1 DO 130 I=1,I1 ANS = ANS + (1.-YY(2,N1-I+1))**(I-1)*ALPHA(N1-I+1)*COMB COMB = COMB*FLOAT(N1-I)/FLOAT(I) 130 CONTINUE COMB = 1. DO 140 J=1,N1 IF (YY(1,J).EQ.1.) GO TO 150 ANS = ANS + (1.-YY(1,J))**(N1-J)*BETA(J)*COMB C COMB IS BINOMIAL COEFFICIENT N OVER J-1. COMB = COMB*FLOAT(N1-J)/FLOAT(J) 140 CONTINUE C FOR PROB OF COMPLEMENTARY EVENT CHANGE NEXT STATEMENT TO C 150 CONTINUE 150 ANS = 1. - ANS RETURN END