ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Main program: computes the powers of the unconditional test (E-test) and c conditional test (C-test) for testing about the means of two Poisson c distributions. c c INPUT: c iside = side of the test; 1 for right-sided, 2 for two-sided c alpha = nominal level of the test c eli = mean of the ith population, i = 1,2 c ni = sample size from the ith lot, i=1,2 c c OUTPUT: c power1 = power of the unconditional test c power2 = power of the conditional test (Przyborowski and Wilenski 1940) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc implicit doubleprecision (a-h,o-z) 1 print*,'Enter 1 for right-tail test or' print*,' 2 for two-tail test' print*,' ' read*, iside if(iside .ne. 1 .and. iside .ne. 2) then goto 1 end if 2 print*, 'Enter the size of the test, alpha' read*, alpha if(alpha .le. 0.0d0 .or. alpha .ge. 1.0d0) then print*, '0 < alpha < 1' goto 2 end if print*, ' ' 3 print*,'Enter the population means lambda1 and lambda2' read *, el1, el2 if(el1 .le. 0.0d0 .or. el2 .le. 0.0d0) then print*, 'el1 > 0, el2 > 0' goto 3 end if print*,' ' if(iside .eq. 1) then if(el1 .lt. el2) then print*,'For the right-tail test lambda1 must be > lambda2' goto 3 end if end if 4 print*,'Enter the sample sizes n1 and n2' read*, n1, n2 if(n1 .le. 0.0 .or. n2 .le. 0.0) then print*, 'n1 > 0, n2 > 0' goto 4 end if print*,' ' 5 print*,'Enter the value of Mean1 - Mean2 under H0:' read*, d if(d .lt. 0.0d0) then print*, 'd >= 0' goto 5 end if print*,' ' ratio = 1.0d0*n1/n2 + (n1*d)/(n2*el2) call poistest(iside, n1, el1, n2, el2, d, ratio, alpha, power1, & power2) print*, 'The power of the E-test is', power1 print*, 'The power of the C-test is', power2 end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Program for computing the power of the E-test and C-test c In the first subroutine, the sum over k1 is carried out cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine poistest(iside, n1, el1, n2, el2, d, ratio, alpha, $ power1, power2) implicit doubleprecision(a-h,o-z) c mean of the Y_1 = sum_1^n1 X_1i is el1*n1 c mean of the Y_2 = sum_1^n2 X_2i is el2*n2 el1n1 = el1*n1 el2n2 = el2*n2 c compute the modes k1mode = int(el1n1) k2mode = int(el2n2) c computing the probability at the k1mode pk1mode = poipr(k1mode, el1n1) pk1 = pk1mode c computing the probability at the k2mode pk2mode = poipr(k2mode, el2n2) power1 = 0.0d0 power2 = 0.0d0 c summing over k1 do k1 = k1mode, 1000 if(pk1 .lt. 1.0d-07) goto 3 call sumk2(iside, n1, k1, pk1, n2, el2n2, k2mode, pk2mode, & d, ratio, alpha, power1, power2) c compute the probability using the forwrad recurrence relation pk1 = el1n1*pk1/(k1+1.0d0) end do c initializing the probablity for backward iteeration 3 k1 = k1mode-1 pk1 = pk1mode pk1 = k1mode*pk1/el1n1 do k1 = k1mode-1, 0, -1 if(pk1 .lt. 1.0d-07) return c call sumk2(iside, n1, k1, pk1, n2, el2n2, k2mode, pk2mode, & d, ratio, alpha, power1, power2) c compute the probability using the backward recurrence relation pk1 = k1*pk1/el1n1 end do end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Here, we carry out the sum over k2 to compute the power of the C-test c and the E-test cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine sumk2(iside, n1, k1, pk1, n2, el2n2, k2mode, pk2mode, & d, ratio, alpha, power1, power2) implicit doubleprecision(a-h,o-z) c computing the probability at the k2mode pk2 = pk2mode sprob = ratio/(1.0d0+ratio) c summing over k2 do k2 = k2mode, 1000 if(pk2 .lt. 1.0d-07) goto 2 c computing estimates of the means el1 and el2 elhatk1 = 1.0d0*k1/n1 elhatk2 = 1.0d0*k2/n2 c computing estimate of the common mean under H_0 elhatk = 1.0d0*(k1+k2)/(n1+n2)-d*n1/(n1+n2) c computing the observed value of the standardized difference; t_k1k2 if(iside .eq. 1) then if(1.0d0*k1/n1 - 1.0d0*k2/n2 .le. d) then pvalue1 = 1.0d0 else var = elhatk1/n1 + elhatk2/n2 t_k1k2 = (elhatk1-elhatk2-d)/sqrt(var) call sumi1(iside, n1, n2, elhatk, t_k1k2, d, alpha, pvalue1) end if c p-value of the conditional test if(k1 .eq. 0) then pvalue2 = 1.0d0 else pvalue2 = betadf(sprob, 1.0d0*k1, k2 + 1.0d0) end if if(pvalue1 .le. alpha) power1 = power1 + pk1*pk2 if(pvalue2 .le. alpha) power2 = power2 + pk1*pk2 else if(iside .eq. 2) then if(dabs(1.0d0*k1/n1 - 1.0d0*k2/n2) .le. d) then pvalue1 = 1.0d0 else var = elhatk1/n1 + elhatk2/n2 t_k1k2 = (elhatk1-elhatk2-d)/sqrt(var) call sumi1(iside, n1, n2, elhatk, t_k1k2, d, alpha, pvalue1) end if c p-value of the conditional test if(k1 .eq. 0 .and. k2 .eq. 0) then pv_rs = 1.0d0 pv_ls = 1.0d0 else if(k1 .ne. 0 .and. k2 .eq. 0) then pv_ls = 1.0d0 pv_rs = betadf(sprob, 1.0d0*k1, k2+1.0d0) else if (k1 .eq. 0 .and. k2 .ne. 0) then pv_rs = 1.0d0 pv_ls = betadf(sprob, 1.0d0*k2, k1+1.0d0) else if (k1 .ne. 0 .and. k2 .ne. 0) then pv_rs = betadf(sprob, 1.0d0*k1, k2+1.0d0) pv_ls = betadf(sprob, 1.0d0*k2, k1+1.0d0) end if pvalue2 = min(pv_rs, pv_ls) if(pvalue1. le. alpha) power1 = power1 + pk1*pk2 if(pvalue2 .le. alpha/2) power2 = power2 + pk1*pk2 end if pk2 = el2n2*pk2/(k2+1.0d0) end do c sum over k2 backward 2 pk2 = pk2mode k2 = k2mode-1 pk2 = k2mode*pk2/el2n2 do k2 = k2mode-1, 0, -1 if(pk2 .lt. 1.0d-07) return elhatk1 = 1.0d0*k1/n1 elhatk2 = 1.0d0*k2/n2 elhatk = 1.0d0*(k1+k2)/(n1+n2)-d*n1/(n1+n2) c computing the observed value of the standardized difference; t_k1k2 if(iside .eq. 1) then if(1.0d0*k1/n1 - 1.0d0*k2/n2 .le. d) then pvalue1 = 1.0d0 else var = elhatk1/n1 + elhatk2/n2 t_k1k2 = (elhatk1-elhatk2-d)/sqrt(var) call sumi1(iside, n1, n2, elhatk, t_k1k2, d, alpha, pvalue1) end if c p-value of the conditional test if(k1 .eq. 0) then pvalue2 = 1.0d0 else pvalue2 = betadf(sprob, 1.0d0*k1, k2 + 1.0d0) end if if(pvalue1 .le. alpha) power1 = power1 + pk1*pk2 if(pvalue2 .le. alpha) power2 = power2 + pk1*pk2 else if(iside .eq. 2) then if(dabs(1.0d0*k1/n1 - 1.0d0*k2/n2) .le. d) then pvalue1 = 1.0d0 else var = elhatk1/n1 + elhatk2/n2 t_k1k2 = (elhatk1-elhatk2-d)/sqrt(var) call sumi1(iside, n1, n2, elhatk, t_k1k2, d, alpha, pvalue1) end if c p-value of the conditional test if(k1 .eq. 0 .and. k2 .eq. 0) then pv_rs = 1.0d0 pv_ls = 1.0d0 else if(k1 .ne. 0 .and. k2 .eq. 0) then pv_ls = 1.0d0 pv_rs = betadf(sprob, 1.0d0*k1, k2+1.0d0) else if (k1 .eq. 0 .and. k2 .ne. 0) then pv_rs = 1.0d0 pv_ls = betadf(sprob, 1.0d0*k2, k1+1.0d0) else if (k1 .ne. 0 .and. k2 .ne. 0) then pv_rs = betadf(sprob, 1.0d0*k1, k2+1.0d0) pv_ls = betadf(sprob, 1.0d0*k2, k1+1.0d0) end if pvalue2 = min(pv_rs, pv_ls) if(pvalue1 .le. alpha) power1 = power1 + pk1*pk2 if(pvalue2 .le. alpha/2) power2 = power2 + pk1*pk2 end if pk2 = k2*pk2/el2n2 end do end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Here, we carry out the sum over i1 to compute the power of the E-test c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine sumi1(iside, n1, n2, elhatk, t_k1k2, d, alpha, pvalue1) implicit doubleprecision(a-h,o-z) pvalue1 = 0.0d0 c computing estimates of el1*n1 and el2*n2 under H_0 elhat1 = n1*(elhatk+d) elhat2 = n2*elhatk c computing the modes i1mode = int(elhat1) i2mode = int(elhat2) c initializing the probability at the i1mode pi1mode = poipr(i1mode, elhat1) pi1 = pi1mode c initializing the probability at the i2mode pi2mode = poipr(i2mode, elhat2) do i1 = i1mode, 1000 if(pi1 .lt. 1.0d-07) goto 1 call sumi2(iside, n1, n2, elhat2, t_k1k2, i1, pi1, i2mode, & pi2mode, d, pvalue1) if(pvalue1 .gt. alpha) return pi1 = elhat1*pi1/(i1+1.0d0) end do 1 i1 = i1mode-1 pi1 = pi1mode pi1 = i1mode*pi1/elhat1 do i1 = i1mode-1, 0, -1 if(pi1 .lt. 1.0d-07) return call sumi2(iside, n1, n2, elhat2, t_k1k2, i1, pi1, i2mode, & pi2mode, d, pvalue1) if(pvalue1 .gt. alpha) return pi1 = i1*pi1/elhat1 end do end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Here, we carry out the sum over i2 to compute the power of the E-test c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine sumi2(iside, n1, n2, elhat2, t_k1k2, i1, pi1, i2mode, & pi2mode, d, pvalue1) implicit doubleprecision(a-h,o-z) pi2 = pi2mode c do i2 = i2mode, 1000 if(pi2 .lt. 1.0d-07) goto 1 elhati1 = 1.0d0*i1/n1 elhati2 = 1.0d0*i2/n2 elhati = 1.0d0*(i1+i2)/(n1+n2)-n1*d/(n1+n2) diffi = elhati1 - elhati2 - d var = elhati1/n1 + elhati2/n2 if(iside .eq. 1) then if(t_k1k2 .gt. 0.0d0 .and. diffi .lt. 0.0d0) goto 1 if(1.0d0*i1/n1 - 1.0d0*i2/n2 .le. d) then t_i1i2 = 0.0d0 else t_i1i2 = diffi/sqrt(var) end if if(t_i1i2 .ge. t_k1k2) pvalue1 = pvalue1 + pi1*pi2 else if(iside .eq. 2) then if(dabs(1.0d0*i1/n1 - 1.0d0*i2/n2) .le. d) then t_i1i2 = 0.0d0 else t_i1i2 = diffi/sqrt(var) end if if(dabs(t_i1i2) .ge. dabs(t_k1k2)) pvalue1 = pvalue1 + pi1*pi2 end if pi2 = elhat2*pi2/(i2+1.0d0) end do c 1 i2 = i2mode-1 pi2 = pi2mode pi2 = i2mode*pi2/elhat2 do i2 = i2mode-1, 0, -1 if(pi2 .lt. 1.0d-07) return elhati1 = 1.0d0*i1/n1 elhati2 = 1.0d0*i2/n2 elhati = 1.0d0*(i1+i2)/(n1+n2)-n1*d/(n1+n2) diffi = elhati1 - elhati2 - d var = elhati1/n1 + elhati2/n2 if(iside .eq. 1) then if(1.0d0*i1/n1 - 1.0d0*i2/n2 .le. d) then t_i1i2 = 0.0d0 else t_i1i2 = diffi/sqrt(var) end if if(t_i1i2 .ge. t_k1k2) pvalue1 = pvalue1 + pi1*pi2 else if(iside .eq. 2) then if(dabs(1.0d0*i1/n1 - 1.0d0*i2/n2) .le. d) then t_i1i2 = 0.0d0 else t_i1i2 = diffi/sqrt(var) end if if(dabs(t_i1i2) .ge. dabs(t_k1k2)) pvalue1 = pvalue1 + pi1*pi2 end if pi2 = i2*pi2/elhat2 end do end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c This program computes the P(X = k), where X is a Poisson random c variable with mean # of defective items = el, and observed # of c defective items = k c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc double precision function poipr(k, el) implicit doubleprecision(a-h,o-z) ek = k*1.0d0 poipr = dexp(-el+ek*dlog(el)-alng(ek+1.0d0)) end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Logarithmic gamma function = alng(x), x > 0 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc double precision function alng(x) implicit doubleprecision (a-h, o-z) double precision b(8) logical indx data b/8.33333333333333d-2, 3.33333333333333d-2, + 2.52380952380952d-1, 5.25606469002695d-1, + 1.01152306812684d0, 1.51747364915329d0, + 2.26948897420496d0, 3.00991738325940d0/ if(x .lt. 8.0d0) then xx = x + 8.0d0 indx = .true. else indx = .false. xx = x end if c fterm = (xx-0.5d0)*dlog(xx) - xx + 9.1893853320467d-1 sum = b(1)/(xx+b(2)/(xx+b(3)/(xx+b(4)/(xx+b(5)/(xx+b(6) +/(xx+b(7)/(xx+b(8)))))))) alng = sum + fterm if(indx) alng = alng-dlog(x+7.0d0)-dlog(x+6.0d0)-dlog + (x+5.0d0)-dlog(x+4.0d0)-dlog(x+3.0d0)-dlog(x+2.0d0) + -dlog(x+1.0d0)-dlog(x) end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c This program computes the P(X .ge. k) where X is a beta random variable c with a = alpha and b = beta c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc double precision function betadf(x, a, b) implicit double precision (a-h, o-z) logical check eps = 1.0d-14 if (x .ge. 1.0) then betadf = 1.0d0 return end if if (x .le. 0.0) then betadf = 0.0d0 return end if if (x .eq. 0.5d0 .and. a .eq. b) then betadf = 0.5d0 return end if betf = alng(a)+alng(b)-alng(a+b) aplusb = a + b omx = 1.0d0 - x if(a .lt. aplusb*x) then y = omx omx = x p = b q = a check = .TRUE. else y = x p = a q = b check= .FALSE. end if ensteps = q + omx*aplusb xovomx = y/omx i = 1 term = 1.0d0 ai = 1.0d0 ans = 1.0d0 1 etermq = q - ai if(ensteps .eq. 0.0d0) xovomx = y 2 term = term*etermq*xovomx/(p+ai) ans = ans + term absterm = dabs(term) if(absterm .le. eps .and. absterm .le. eps*ans & .or. i .gt. 1000) goto 3 ai = ai + 1.0d0 ensteps = ensteps-1.0d0 i = i + 1 if(ensteps .ge. 0.0d0) goto 1 etermq = aplusb aplusb = aplusb + 1.0d0 goto 2 3 ans = ans*dexp(p*dlog(y)+(q-1.0d0)*dlog(omx)-betf)/p if(check) ans = 1.0d0-ans betadf = ans end