ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  Name: FinitePs.for 
c
c This program computes the powers of an unconditional test for testing the
c difference between proportions in finite populations. 
c
c INPUT:
c iside = side of the test; 1 for right-sided; 2 for two-sided
c alpha = nominal level of the test; takes values in (0, 0.5)
c loti = lot size of the ith lot, i=1,2
c pi = proportion of defective items in the ith lot, i=1,2
c ni = sample size from the ith lot, i=1,2
c
c OUTPUT:
c power = power of the test
c 
c Reference: Krishnamoorthy, K. and Thomson, J. (2002) Hypothesis testing about proportions 
c in two finite populations.  The American Statistician,  56, 215-222.
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	implicit doubleprecision (a-h,o-z)

	integer n1, lot1, n2, lot2
	common //iside

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 level of the test; level must be in (0, 0.5)'
	read *, alpha
	if(alpha .ge. 0.5) goto 2
	print *,' '

3     print *,'Enter the population proportions p1 and p2'
	read *, p1, p2
	if(p1 .le. 0.0d0 .or. p1 .ge. 1.0d0 .or. 
     &   p2 .le. 0.0d0 .or. p2 .ge. 1.0d0) then
      print *, '0 < p1 < 1, 0 < p2 < 1'
	goto 3
	end if
	print *,' '

	if(iside .eq. 1) then
	if(p1 .lt. p2) then
	print*,'For the right-tail test p1 must be > p2'
	goto 3
	end if
	end if

      print *,'Enter the sizes of the lot1 and lot2'
	read *, lot1, lot2
	print *,' '


4     print *,'Enter the sample sizes n1 and n2'
	read *, n1, n2
	if(n1 .le. 0.0 .or. n1 .ge. lot1 .or. n2 .le. 0.0 .or. n2.ge.lot2) 
     &   then
      print *, '0 < n1 < lot1, 0 < n2 < lot2'
	goto 4
	end if
	print *,' '
	
      call etest(n1, p1, lot1, n2, p2, lot2, alpha, power)
	print *, 'Power of the E-test = ',power

	end 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Program for computing the power of the E-test
c In the first subroutine, the sum over k1 is carried out
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	subroutine etest(n1, p1, lot1, n2, p2, lot2, alpha, power)
	implicit doubleprecision(a-h,o-z)

	integer lot1, m1, n1, lot2, m2, n2
	integer L1, IU1, L2, IU2

c setup the range for k1 and k2

	m1 = int(p1*lot1)
	m2 = int(p2*lot2)
	L1 = max(0, m1-lot1+n1)
	IU1 = min(n1,m1)
	L2 = max(0, m2-lot2+n2)
	IU2 = min(n2,m2)

c determine the mode for summing over k1 and k2

	k1mode = int((n1+1.0d0)*(m1+1.0)/(lot1+2.0d0))
	k2mode = int((n2+1.0d0)*(m2+1.0)/(lot2+2.0d0)) 

c finite population correction

	fc1 = (lot1-n1)/(lot1-1.0d0)/n1
	fc2 = (lot2-n2)/(lot2-1.0d0)/n2

c computing the probability at the k1mode
	pk1mode = hyprob(k1mode, n1, m1, lot1)
	pk1 = pk1mode

c computing the probability at the k2mode
	pk2mode = hyprob(k2mode, n2, m2, lot2)

	power = 0.0d0

c summing over k1 (forward)
      do k1 = k1mode, IU1
		if(pk1 .lt. 1.0d-14) goto 3

		call sumk2(k1, n1, lot1, pk1, fc1, k2mode, n2, m2, lot2, 
     &                 pk2mode, fc2, L2, IU2, alpha, power)

c compute the probability using the forwrad recurrence relation
		pk1 = (n1-k1)*(m1-k1)*pk1/((k1+1.0d0)*(lot1-m1-n1+k1+1.0d0))
	end do

c initializing the probablity for backward iteration
3     k1 = k1mode-1
	pk1 = pk1mode
      pk1 = k1mode*(lot1-m1-n1+k1mode)*pk1
     &    /((n1-k1mode+1.0d0)*(m1-k1mode+1.0d0))
	do k1 = k1mode-1, L1, -1
		if(pk1 .lt. 1.0d-14) return
c 
		call sumk2(k1, n1, lot1, pk1, fc1, k2mode, n2, m2, lot2, 
     &                 pk2mode, fc2, L2, IU2, alpha, power)

c compute the probability using the backward recurrence relation
		pk1 = k1*(lot1-m1-n1+k1)*pk1/((n1-k1+1.0d0)*(m1-k1+1.0d0))
	end do
	end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Here, we carry out the sum over k2 to compute the power of the E-test 
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine sumk2(k1, n1, lot1, pk1, fc1, k2mode, n2, m2, lot2, 
     &                 pk2mode, fc2, L2, IU2, alpha, power)
	implicit doubleprecision(a-h,o-z)

	integer IUi1, IUi2, Li1, Li2, IU2, L2, k1, n1, m1, lot1
	integer k2mode, n2, m2, lot2

c computing the probability at the k2mode
	pk2 = pk2mode

	do k2 = k2mode, IU2
	if(pk2 .lt. 1.0d-14) goto 2

	phatk1 = 1.0d0*k1/n1
	phatk2 = 1.0d0*k2/n2

      phat = 1.0*(k1+k2)/(n1+n2)
	qhat=1.0-phat

	m1hat = int((lot1)*phat)
	m2hat = int((lot2)*phat)

	if(k1+k2 .eq. 0 .or. k1+k2 .eq. n1+n2) then
	t_k1k2 = 0.0d0
	else
	t_k1k2 = (phatk1-phatk2)/sqrt((fc1+fc2)*phat*qhat)
	end if

	Li1 = max(0, m1hat-lot1+n1)
	IUi1 = min(n1,m1hat)
	Li2 = max(0, m2hat-lot2+n2)
	IUi2 = min(n2,m2hat)
	call sumi1(n1, lot1, m1hat, Li1, IUi1, fc1, n2, lot2, m2hat, 
     &           Li2, IUi2, fc2, t_k1k2, pvalue)

	if(pvalue .le. alpha) power = power + pk1*pk2
	
	pk2 = (n2-k2)*(m2-k2)*pk2/((k2+1.0d0)*(lot2-m2-n2+k2+1.0d0))
 	end do

c sum over k2 backward

2	pk2 = pk2mode
	k2 = k2mode-1
      pk2 = k2mode*(lot2-m2-n2+k2mode)*pk2
     &    /((n2-k2mode+1.0d0)*(m2-k2mode+1.0d0))
	
	do k2 = k2mode-1, L2, -1
	if(pk2 .lt. 1.0d-14) return
	
	phatk1 = 1.0d0*k1/n1
	phatk2 = 1.0d0*k2/n2
      phat = 1.0d0*(k1+k2)/(n1+n2)
	qhat=1.0d0-phat

c estimte the number of defs in the lots
	m1hat = int((lot1)*phat)
	m2hat = int((lot2)*phat)


c compute the observed value of the std difference 
	if(k1+k2 .eq. 0 .or. k1+k2 .eq. n1+n2) then
	t_k1k2 = 0.0d0
	else
	t_k1k2 = (phatk1-phatk2)/sqrt((fc1+fc2)*phat*qhat)
	end if

	Li1 = max(0, m1hat-lot1+n1)
	IUi1 = min(n1,m1hat)
	Li2 = max(0, m2hat-lot2+n2)
	IUi2 = min(n2,m2hat)

	call sumi1(n1, lot1, m1hat, Li1, IUi1, fc1, n2, lot2, m2hat, 
     &           Li2, IUi2, fc2, t_k1k2, pvalue)

	if(pvalue .le. alpha) power = power + pk1*pk2
	pk2 = k2*(lot2-m2-n2+k2)*pk2/((n2-k2+1.0d0)*(m2-k2+1.0d0))
	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(n1, lot1, m1hat, Li1, IUi1, fc1, n2, lot2, m2hat, 
     &           Li2, IUi2, fc2, t_k1k2, pvalue)
	implicit doubleprecision(a-h,o-z)
	
	integer Li1, Li2, IUi1, IUi2, lot1, lot2, m1hat, m2hat
	integer n1, n2
	pvalue = 0.0d0

	i1mode = int((n1+1.0d0)*(m1hat+1.0d0)/(lot1+2.0d0))
	i2mode = int((n2+1.0d0)*(m2hat+1.0d0)/(lot2+2.0d0))

	pi1mode = hyprob(i1mode, n1, m1hat, lot1)
	pi1 = pi1mode

	pi2mode = hyprob(i2mode, n2, m2hat, lot2)

	do i1 = i1mode, IUi1
		if(pi1 .lt. 1.0d-14) goto 1
		call sumi2(n1, fc1, n2, lot2, m2hat, fc2, Li2, IUi2, 
     &			  t_k1k2, i1, pi1, i2mode, pi2mode, pvalue)
		pi1 = (n1-i1)*(m1hat-i1)*pi1/((i1+1.0d0)*(lot1-m1hat-n1+i1+1.0d0))
	end do

1	i1 = i1mode-1
	pi1 = pi1mode
	pi1 = i1mode*(lot1-m1hat-n1+i1mode)*pi1
     &     /((n1-i1mode+1.0d0)*(m1hat-i1mode+1.0d0))
	
	do i1 = i1mode-1, Li1, -1
		if(pi1 .lt. 1.0d-14) return
		call sumi2(n1, fc1, n2, lot2, m2hat, fc2, Li2, IUi2, 
     &              t_k1k2, i1, pi1, i2mode, pi2mode, pvalue)
		pi1 = i1*(lot1-m1hat-n1+i1)*pi1/((n1-i1+1.0d0)*(m1hat-i1+1.0d0))
	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(n1, fc1, n2, lot2, m2hat, fc2, Li2, 
     &                IUi2, t_k1k2, i1, pi1, i2mode, pi2mode, pvalue)
	implicit doubleprecision(a-h,o-z)

c
	integer Li2, IUi2
	common //iside

c	
	pi2 = pi2mode
c
	do i2 = i2mode, IUi2
		if(pi2 .lt. 1.0d-14) goto 1
		phati1 = 1.0d0*i1/n1
		phati2 = 1.0d0*i2/n2
		phati = 1.0d0*(i1+i2)/(n1+n2)
		qhati = 1.0d0-phati
		if(i1+i2 .eq. 0 .or. i1+i2 .eq. n1+n2) then
			t_i1i2 = 0.0d0
		else
			t_i1i2 = (phati1-phati2)/sqrt((fc1+fc2)*phati*qhati)
		end if
	    if(iside .eq. 1) then
			if(t_i1i2 .ge. t_k1k2) pvalue = pvalue + pi1*pi2
		else
			if(dabs(t_i1i2) .ge. dabs(t_k1k2)) pvalue = pvalue + pi1*pi2
		end if
		pi2 = (n2-i2)*(m2hat-i2)*pi2/((i2+1.0d0)*(lot2-m2hat-n2+i2+1.0d0))
	end do
c
1     i2 = i2mode-1 
	pi2 = pi2mode
	pi2 = i2mode*(lot2-m2hat-n2+i2mode)*pi2
     &   /((n2-i2mode+1.0d0)*(m2hat-i2mode+1.0d0))

	do i2 = i2mode-1, Li2, -1
		if(pi2 .lt. 1.0d-14) return
		phati1 = 1.0d0*i1/n1
		phati2 = 1.0d0*i2/n2
		phati = 1.0d0*(i1+i2)/(n1+n2)
		qhati = 1.0d0-phati
		if(i1+i2 .eq. 0 .or. i1+i2 .eq. n1+n2) then
			 t_i1i2 = 0.0d0
		else
			 t_i1i2 = (phati1-phati2)/sqrt((fc1+fc2)*phati*qhati)
		end if
	    if(iside .eq. 1) then
			if(t_i1i2 .ge. t_k1k2) pvalue = pvalue + pi1*pi2
		else
			if(dabs(t_i1i2) .ge. dabs(t_k1k2)) pvalue = pvalue + pi1*pi2
		end if
		pi2 = i2*(lot2-m2hat-n2+i2)*pi2/((n2-i2+1.0d0)*(m2hat-i2+1.0d0))
	end do
	end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c END: of the program for computing the power of the E-test 
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c This program computes the P(X = k), where X is a hypergeometric random
c variable with lot size = lot, # of defective items = m, sample size = n
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	double precision function hyprob(k, n, m, lot)
	implicit doubleprecision(a-h,o-z)
	integer L, U

	L = max(0, m-lot+n)
	U = min(n,m)

	if(k .lt. L .or. k .gt. U) then
	hyprob = 0.0d0 
	return
	end if
	  
	ek = k
	en = n
	em = m
	elot = lot

	term1 = alng(em+1.0d0)-alng(ek+1.0d0)-alng(em-ek+1.0d0)
	term2 = alng(elot-em+1.0d0)-alng(en-ek+1.0d0)
     &   -alng(elot-em-en+ek+1.0d0)
	term3 = alng(elot+1.0d0)-alng(en+1.0d0)-alng(elot-en+1.0d0)
	hyprob = dexp(term1+term2-term3)

	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