#Poisson E-test to calculate p-value for difference in two Poisson means
#See #http://www.ucs.louisiana.edu/~kxk4695/JSPI-04.pdf
#A more powerful test for comparing two Poisson means
#Krishamoorthy & Thomson, 2004, Journal of Stat Planning and Inf 119, 23-35

#Questions? email apwheele@gmail.com
#Andy Wheeler, https://andrewpwheeler.wordpress.com/

#this function finds the integer for poisson
#PMF that has a density of less than eps
#where to terminate the sum
minPMF <- function(r,eps=1e-20){
	p <- 1
	int <- trunc(r) #only evaluate right tail
	while(p > eps){
		int <- int + 1
		p <- dpois(int,r)
	}
	return(int)
}

#Tests the null k1/n1 = k2/n2 + d
#When d = 0, can switch k1/n1 & k2/n2 and still get the same p-value
#with d != 0 though it makes a difference
poisson.etest <- function(k1,k2,n1=1,n2=1,d=0,eps=1e-20){
	r1 <- k1/n1 #rate first process
	r2 <- k2/n2 #rate second process
	lhat <- (k1 + k2)/(n1 + n2) - (d*n1)/(n1 + n2) #mean rate
	Tk <- abs((r1 - r2 - d)/sqrt(lhat)) #test statistic
	d1 <- n2*lhat #updated so termination of sum did not happen too early
	d2 <- n1*(lhat+d) #problem when using lhat instead of d's
	int1 <- minPMF(r=d1,eps=eps) #where to terminate sum
	if (d==0){int2 <- int1} else {int2 <- minPMF(r=d2,eps=eps)}
	df <- expand.grid(x1 = 0:int1,x2 = 0:int2) #for really large lambdas
	df$p1 <- dpois(df$x1,d2)          #this might not fit in memory 
	df$p2 <- dpois(df$x2,d1)     
	df$rp <- (df$x1 + df$x2)/(n1 + n2)
	df$Tx <- abs((df$x1/n1 - df$x2/n2 - d))/sqrt(df$rp)
	df$I <- 1*(df$Tx >= Tk)
	df$ptot <- df$p1 * df$p2 * df$I
	return(sum(df$ptot, na.rm=TRUE)) #x1 = 0 and x2 = 0 produces nulls for Tx, rp=0
}
#Need to update for 0,0 case, which is undefined so should NOT return 0
	
###################################################################################################
#poisson.etest(3,0)
#poisson.etest(3,0,2,2)
#poisson.etest(3,0,100,100) #updated so stopping sums earlier is not a problem
#
#poisson.etest(6,2) #second example from article
#poisson.etest(2,6) #when d=0, can switch arguments and get the same p-value
#
##I would think these two should be the same but they arent
##Might be my error
#poisson.etest(6,2,d=1)
#poisson.etest(2,6,d=-1)
###################################################################################################