OddsRatioCBs <- function(Table,alpha=0.05)
# OddsRatioCBs(Table,alpha= 0.05)
# computes for a 2x2 contingency table the boundaries 
# of the (1-alpha)-confidence intervals
#    [a1,inf],  [0,b1],  [a2,b2]
# as well as the standard estimator rhohat for an 
# unknown odds ratio.
# 
# Lutz Duembgen, May 22, 2009
{
	# Row and column totals:
	z <- sum(Table[1,]);
	s <- sum(Table[,1]);
	n = sum(Table);
	
	# MLE:
	if (sum(Table == 0) > 0)
	{
		rhohat <- (Table[1,1]+0.5)*(Table[2,2]+0.5)/(Table[1,2]+0.5)/(Table[2,1]+0.5)
	}
	else
	{
		rhohat <- Table[1,1]*Table[2,2] / Table[1,2] / Table[2,1]
	}
	
	# Test statistic for exact confidence bounds:
	T = Table[1,1];
	
	# Exact confidence bounds:
	a1 <- 1/OddsRatioUCB(z-T,n,z,n-s,alpha)
	b1 <- OddsRatioUCB(T,n,z,s,alpha)
	a2 <- 1/OddsRatioUCB(z-T,n,z,n-s,alpha/2)
	b2 <- OddsRatioUCB(T,n,z,s,alpha/2)
	
	# Asymptotic confidence bounds:
	kappa1 <- qnorm(1 - alpha)
	kappa2 <- qnorm(1 - alpha/2)
	if (sum(Table == 0) > 0)
	{
		shut <- sqrt(sum(1/(Table + 0.5)))
	}
	else
	{
		shut <- sqrt(sum(1/Table))
	}
	c1 <- rhohat*exp(- shut * kappa1)
	d1 <- rhohat*exp(shut * kappa1)
	c2 <- rhohat*exp(- shut * kappa2)
	d2 <- rhohat*exp(shut * kappa2)
	
	return(list(estimator=rhohat,lower.bound=a1,upper.bound=b1,conf.int=c(a2,b2),lower.bound.appr=c1,upper.bound.appr=d1,conf.int.appr=c(c2,d2)))
}


OddsRatioUCB <- function(T,n,z,s,beta)
# returns a number b such that in case of  
# max(z+s-n,0) <= T < min(z,s):
#    0  <=  beta - F(b,T,n,z,s)  <=  delta ,
#    F(b - delta*b, T,n,z,s) >= beta .
# Here delta is a given precision (at present 0.0001), 
# and
#    F(rho,j,n,z,s)  =  C  sum_{j=0}^k 
#                     {z choose j} {n-z choose s-j} rho^j 
# with C > 0 such that F(rho,min(z,s),n,z,s) = 1.
# The underlying method is binary search.
# 
# Lutz Duembgen, May 22, 2009
{
	delta <- 0.00001
	if (T >= min(z,s))
	{
		b <- Inf
		# only trivial upper bound
		return(b)
	}
	a  <- 0
	pa <- 1
	b  <- 1
	pb <- EWHygeCDF(T,n,z,s,b)
	
	# Look for interval [a,b] such that pa > beta >= pb:
	while (pb > beta)
	{
		a  <-  b
		pa <- pb
		b  <- 2*b
		pb <- EWHygeCDF(T,n,z,s,b)
	}
	
	# Split the interval [a,b] into two halves until 
	# the desired precision is reached:
	while ((b - a > delta) || (pa - pb > delta))
	{
		tmp  <- (a+b)/2
		ptmp <- EWHygeCDF(T,n,z,s,tmp)
		if (ptmp > beta)
		{
			a  <-  tmp
			pa <- ptmp
		}
		else
		{
			b  <-  tmp
			pb <- ptmp
		}
	}
	return(b)
}


EWHygeCDF <- function(x,n,z,s,rho=1)
# computes the exponentially weighted hypergeometric
# distribution function with parameter rho for each component
# of x.
# 
# Lutz Duembgen, May 22, 2009
{
	x2 <- floor(x)
	F  <- rep(0,times=length(x))
	F[x2 >= min(z,s)] <- 1
	jj <- which(x2 >= 0 & x2 < min(z,s))
	if (length(jj) > 0)
	{
		gv <- EWHygePDF(n,z,s,rho)
		gv <- cumsum(gv)
		F[jj] <- gv[x2[jj]+1]
	}
	return(F)
}


EWHygePDF <- function(n,z,s,rho=1)
# returns a vector gv of length min(z,s)+1 with 
#    gv(i+1)  =  C {L choose i} {n-L choose m-i} rho^i ,
# where C > 0 is chosen such that sum(gv) = 1. 
# 
# Lutz Duembgen, May 22, 2009
{
	kmax <- min(z,s)
	kmin <- max(0,z+s-n)
	if (kmin == kmax)
	{
		return(rep(0,times=kmin),1)
	}
	if (rho == 0)
	{
		return(c(rep(0,times=kmin),1,rep(0,times=kmax-kmin)))
	}
	xv <- ((kmin+1):kmax)
	gv <- rep(0,times=kmax-kmin)
	gv <- log(rho) + log(z-xv+1) + log(s-xv+1) - log(xv) - log(n-z-s+xv)
	gv <- c(0,cumsum(gv))
	gv <- exp(gv - max(gv))
	gv <- gv/sum(gv)
	return(c(rep(0,times=kmin),gv))
}