C THIS PROGRAM PERFORMS YOO AND SPOTH'S 1993 AND MIETTENEN'S 1968
C PROCEDURES FOR DETERMINING SAMPLE SIZES FOR BINOMIALLY
C DISTRIBUTED OUTCOMES
C
        OPEN(UNIT=1,FILE='/CONBIN.OUT',STATUS='UNKNOWN')
C
C       INITIALIZE VARIABLES
C
1       A=0.0
        AL=0.0
        APC1=0.
        APC2=0.
        APT1=0.
        DELTA=0.0
        EFF=0.
        EXPT2CON=0.0
        EXPT2TRT=0.0
        H=0.
        QUANTA=0.0
        QZ1=0.0
        QZ1I=0.0
        SI=0.0
        VARDIFFI=0.
        VARDIFFM=0.
        VARIND1=0.
        VARIND2=0.
        VARMAT1=0.
        VARMAT2=0.
        X2I=0.0
        XDEN=0.0
        XDENOM=0.0
        XDENOMI=0.0
        XL1=0.0
        XL2=0.0
        XMIET=0.0
        XNBIMATCH=0.0
        XNBIMATCHI=0.0
        XNIND=0.0
        XNUM=0.0
        XNUM1=0.0
        XNUM2=0.0
        XNUMER=0.0
        XNUMERI=0.0
        XPB=0.0
        XPBI=0.0
        XP1=0.0
        XP2=0.0
        XPC=0.0
        XPCONAD=0.0
        XPCONADJ=0.0
        XPCONADJI=0.0
        XPOW=0.0
        XPTRTAD=0.0
        XPTRTADJ=0.0
        XPTRTADJI=0.0
        XRATIO=0.0
        XRATIOI=0.0
        YCON=0.0
        YTRT=0.0
        ZA=0.0
        ZAL=0.0
        ZB=0.0
        ZBI=0.0
        ZBN=0.0
        ZBNI=0.0
        ZP=0.0
C
       WRITE(5,2)
2      FORMAT(/' ENTER THE PROPORTION OF CONTROL GROUP  2 QUITS')
       READ(5,*)XP1
       IF(XP1.EQ.2)GO TO 29
       WRITE(5,3)
3      FORMAT(' ENTER THE PROPORTION OF TREATMENT GROUP ')
       READ(5,*)XP2
       WRITE(5,4)
4      FORMAT(' ENTER THE PROPORTION OF CHANGE ')
       READ(5,*)XPC
       WRITE(5,5)
5      FORMAT(' ENTER YOUR ALPHA LEVEL')
       READ(5,*)AL
       WRITE(5,6)
6      FORMAT(' ENTER YOUR DESIRED POWER')
       READ(5,*)XPOW
C
C      CALCULATE EFFECT SIZE INDEX
C
        EFF = XPC * (XP2-XP1)
C
C EXPECTED TIME 2 PROPORTION FOR THE TIME 1 CONTROL GROUP
C
        EXPT2CON=XP2
C
C EXPECTED TIME 2 PROPORTION FOR THE TIME 1 TREATMENT GROUP
C
        EXPT2TRT=XP2-EFF
C
       WRITE(5,7)
7      FORMAT(' FOR INDEPENDENT GROUPS, ENTER 1 ',
     +/,  ' FOR MATCHED PAIRS, ENTER 2 ')
       READ(5,*)A
       IF(A.EQ.2.) GOTO 18
C
       WRITE(1,8) XP1
       WRITE(5,8) XP1
8      FORMAT(' THE PROPORTION OF THE CONTROL GROUP = ', F10.4)
       WRITE(1,9) XP2
       WRITE(5,9) XP2
9      FORMAT(' THE PROPORTION OF THE TREATMENT GROUP = ', F8.4)
       WRITE(1,10) XPC
       WRITE(5,10)XPC
10     FORMAT(' THE PROPORTION OF CHANGE = ' ,F10.4)
       WRITE(1,11)AL
       WRITE(5,11)AL
11     FORMAT(' THE ALPHA = ',F10.4)
       WRITE(1,12)XPOW
       WRITE(5,12)XPOW
12     FORMAT(' THE POWER = ',F10.4)
       WRITE(1,13)EFF
       WRITE(5,13)EFF
13     FORMAT(' THE EFFECT SIZE = ',F10.4)
       WRITE(1,14)EXPT2CON
       WRITE(5,14)EXPT2CON
14     FORMAT(' EXPECTED TIME 2 PROPORTION FOR THE ',
     +  'CONTROL GROUP = ',F10.4)
       WRITE(1,15)EXPT2TRT
       WRITE(5,15)EXPT2TRT
15     FORMAT(' EXPECTED TIME 2 PROPORTION FOR THE ',
     +  'TREATMENT GROUP = ',F8.4)
C
C       ESTIMATED n PER GROUP FOR TWO INDEPENDENT PROPORTIONS
C
        YCON=ASIN(SQRT(EXPT2CON))
        YTRT=ASIN(SQRT(EXPT2TRT))
        H=2*(ABS(YCON-YTRT))
        ZAL=1-AL
        ZA=PPND7(ZAL,XPA)
        QZ1=PPND7(XPC,XXP)
        QZ1I=-1*QZ1
        ZP=PPND7(XPOW,XP)
        XNIND = (2*(ZP+ZA/H))**2
C
        WRITE(1,16)XNIND
        WRITE(5,16)XNIND
16      FORMAT(' THE N PER GROUP FOR TWO INDEPENDENT PROPORTIONS ',
     +  'USING THE ARCSINE '
     +/,' TRANSFORMATION METHOD =  ',F10.4)
C
C       ADJUSTED PROPORTION RATE FOR THE CONTROL GRP (INDEP GRPS)
C
        XL1=XNIND*(XP2-XP1)
        XL2=XNIND*XP1
        XPCONADJI=XL1/(XNIND-XL2)
C
C       ADJUSTED PROPORTION RATE FOR THE TREATMENT GRP (INDEP GRPS)
C
        X2I=XNIND*((XP2-XP1)-EFF)
        XPTRTADJI=X2I/(XNIND-XL2)
C
C       OTHER QUANTITIES COMPUTED FOR YOO AND SPOTH (INDEP GRPS)
C
        XPBI=XL2/XNIND
        VARIND1=(XPCONADJI*(1-XPCONADJI))/(XNIND-XL2)
        VARIND2=(XPTRTADJI*(1-XPTRTADJI))/(XNIND-XL2)
        VARDIFFI=VARIND1+VARIND2
C
C       COMPUTE Z BETA FOR INDEPENDENT GROUPS
C
        ZBI=ZA-((XPCONADJI-XPTRTADJI)/SQRT(VARDIFFI))
        ZBNI=-1*ZBI
C
C       COMPUTE N PER GRP VIA CONDITIONAL BINOMIAL METHOD
C       FOR INDEPENDENT PROPORTIONS
C
        XNUMERI=(XPCONADJI*(1-XPCONADJI))+(XPTRTADJI*(1-XPTRTADJI))
        XDENOMI=((XPCONADJI-XPTRTADJI)**2)*(1-XPBI)
        XRATIOI=XNUMERI/XDENOMI
        XNBIMATCHI=((QZ1I+ZBNI)**2)*XRATIOI
        WRITE(1,17)XNBIMATCHI
        WRITE(5,17)XNBIMATCHI
17      FORMAT(' THE N PER GROUP FOR TWO INDEPENDENT PROPORTIONS ',
     +  'USING THE CONDITIONAL '
     +/,' BINOMIAL METHOD = ',F10.4)
        GO TO 1
C
18     WRITE(1,19) XP1
       WRITE(5,19) XP1
19     FORMAT(' THE PROPORTION OF THE CONTROL GROUP = ', F10.4)
       WRITE(1,20) XP2
       WRITE(5,20) XP2
20     FORMAT(' THE PROPORTION OF THE TREATMENT GROUP = ', F8.4)
       WRITE(1,21) XPC
       WRITE(5,21) XPC
21     FORMAT(' THE PROPORTION OF CHANGE = ' ,F10.4)
       WRITE(1,22)AL
       WRITE(5,22)AL
22     FORMAT(' THE ALPHA = ',F10.4)
       WRITE(1,23)XPOW
       WRITE(5,23)XPOW
23     FORMAT(' THE POWER = ',F10.4)
       WRITE(1,24)EFF
       WRITE(5,24)EFF
24     FORMAT(' THE EFFECT SIZE = ',F10.4)
       WRITE(1,25)EXPT2CON
       WRITE(5,25)EXPT2CON
25     FORMAT(' EXPECTED TIME 2 PROPORTION FOR THE',
     + ' CONTROL GROUP = ',F10.4)
       WRITE(1,26)EXPT2TRT
       WRITE(5,26)EXPT2TRT
26     FORMAT(' EXPECTED TIME 2 PROPORTION FOR THE',
     + ' TREATMENT GROUP = ',F8.4)
C
C       ESTIMATED n PER GROUP USING MIETTINEN'S (1968) PROCEDURE
C
        DELTA = EXPT2TRT-EXPT2CON
        ZP=PPND7(XPOW,XP)
        SI=(EXPT2TRT*(1-EXPT2CON))+((1-EXPT2TRT)*EXPT2CON)
        QUANTA=1-AL
        ZAL=PPND7(QUANTA,XPA1)
        XNUM1=(ZAL*SI)
        XNUM2=(ZP*SQRT(SI**2-(0.25*DELTA**2*(3+SI))))
        XNUM=(XNUM1+XNUM2)**2
        XDEN=(SI*DELTA**2)
        XMIET=XNUM/XDEN
C
        WRITE(1,27)XMIET
        WRITE(5,27)XMIET
27      FORMAT(' THE N PER GROUP FOR MATCHED PAIRS VIA THE ',
     +  'MIETTINEN PROCEDURE = ',F10.4)
C
C      ADJUSTED PROPORTION RATE FOR CONTROL GRP (MATCHED PAIRS)
C
        APC1=XMIET*(XP2-XP1)
        APC2=XMIET*XP1
        XPCONAD=APC1/(XMIET-APC2)
C
C      ADJUSTED PROPORTION RATE FOR TRT GROUP (MATCHED PAIRS)
C
        APT1=XMIET*((XP2-XP1)-EFF)
        XPTRTAD=APT1/(XMIET-APC2)
C
C       OTHER QUANTITIES COMPUTED FOR YOO AND SPOTH (MATCHED PAIRS)
C
        XPB=APC2/XMIET
        VARMAT1=(XPCONAD*(1-XPCONAD))/(XMIET-APC2)
        VARMAT2=(XPTRTAD*(1-XPTRTAD))/(XMIET-APC2)
        VARDIFFM=VARMAT1+VARMAT2
C
C       COMPUTE Z BETA FOR MATCHED PAIRS
C
        ZB=ZAL-((XPCONAD-XPTRTAD)/SQRT(VARDIFFM))
        ZBN=-1*ZB
C
C
C       COMPUTE N PER GRP VIA CONDITIONAL BINOMIAL METHOD
C       FOR THE MATCHED PAIRS DESIGN
C
        XNUMER=(XPCONAD*(1-XPCONAD))+(XPTRTAD*(1-XPTRTAD))
        XDENOM=((XPCONAD-XPTRTAD)**2)*(1-XPB)
        XRATIO=XNUMER/XDENOM
        QZ=PPND7(XPC,AVV)
        QZ1=-1*QZ
        XNBIMATCH=((QZ1+ZBN)**2)*XRATIO
        WRITE(1,28)XNBIMATCH
        WRITE(5,28)XNBIMATCH
28      FORMAT(' THE N PER GROUP FOR MATCHED PAIRS VIA THE ',
     +  'CONDITIONAL BINOMIAL METHOD = ',F8.4)
C
        GO TO 1
29      WRITE(5,30)
30      FORMAT(' A SUMMARY IS CONTAINED IN CONBIN.OUT')
        STOP
        END
C
C     -------------------------------------------------------------
C

C
       REAL FUNCTION PPND7 (P, IFAULT)
C
C	ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3, 477-
C   484.
C
C	Produces the normal deviate Z corresponding to a given lower
C	tail area of P; Z is accurate to about 1 part in 10**7.
C
C	The hash sums below are the sums of the mantissas of the
C	coefficients.   They are included for use in checking
C	transcription.
C
       REAL ZERO, ONE, HALF, SPLIT1, SPLIT2, CONST1, CONST2, A0, A1,
     +      A2, A3, B1, B2, B3, C0, C1, C2, C3, D1, D2, E0, E1, E2,
     +      E3, F1, F2, P, Q, R
       PARAMETER (ZERO = 0.0, ONE = 1.0, HALF = 0.5,
     +      SPLIT1 = 0.425, SPLIT2 = 5.0,
     +      CONST1 = 0.180625, CONST2 = 1.6)
C
C	Coefficients for P close to 0.5
C
       PARAMETER (A0 = 3.38713 27179E+00, A1 = 5.04342 71938E+01,
     +         A2 = 1.59291 13202E+02, A3 = 5.91093 74720E+01,
     +         B1 = 1.78951 69469E+01, B2 = 7.87577 57664E+01,
     +         B3 = 6.71875 63600E+01)
C	HASH SUM AB    32.31845 77772
C
C	Coefficients for P not close to 0, 0.5 or 1.
C
       PARAMETER (C0 = 1.42343 72777E+00, C1 = 2.75681 53900E+00,
     +         C2 = 1.30672 84816E+00, C3 = 1.70238 21103E-01,
     +         D1 = 7.37001 64250E-01, D2 = 1.20211 32975E-01)
C	HASH SUM CD    15.76149 29821
C
C	Coefficients for P near 0 or 1.
C
       PARAMETER (E0 = 6.65790 51150E+00, E1 = 3.08122 63860E+00,
     +         E2 = 4.28682 94337E-01, E3 = 1.73372 03997E-02,
     +         F1 = 2.41978 94225E-01, F2 = 1.22582 02635E-02)
C	HASH SUM EF    19.40529 10204
C
       IFAULT = 0
       Q = P - HALF
       IF (ABS(Q) .LE. SPLIT1) THEN
       R = CONST1 - Q * Q
       PPND7 = Q * (((A3 * R + A2) * R + A1) * R + A0) /
     +            (((B3 * R + B2) * R + B1) * R + ONE)
       RETURN
       ELSE
       IF (Q .LT. ZERO) THEN
	    R = P
	  ELSE
	    R = ONE - P
	  END IF
	  IF (R .LE. ZERO) THEN
	    IFAULT = 1
	    PPND7 = ZERO
	    RETURN
	  END IF
	  R = SQRT(-LOG(R))
	  IF (R .LE. SPLIT2) THEN
	    R = R - CONST2
	    PPND7 = (((C3 * R + C2) * R + C1) * R + C0) /
     +           ((D2 * R + D1) * R + ONE)
	  ELSE
	    R = R - SPLIT2
	    PPND7 = (((E3 * R + E2) * R + E1) * R + E0) /
     +           ((F2 * R + F1) * R + ONE)
	  END IF
	  IF (Q .LT. ZERO) PPND7 = - PPND7
	  RETURN
      END IF
      END


