BinoCBs <- function(y, n, alpha=0.05)
# Computes exact one- and two-sided confidence bounds 
# for the unknown parameter p of the binomial distribution 
# Bin(n,p), given an observation y in {0,1,...,n}.
# For a random variable Y with distribution Bin(n,p), 
#   Pr{p >= lower.bound(Y,n,alpha)}  >=  1 - alpha ,
#   Pr{p <= upper.bound(Y,n,alpha)}  >=  1 - alpha ,
# and
#   Pr{p within confidence.interval(Y,n,alpha)}  >=  alpha .
# 
# Lutz Duembgen, August 28, 2008
{
	if ((alpha <= 0) || (y < 0) || (alpha >= 1) || (y > n))
	{
		print("Input is nonsense!")
		return()
	}

	# Precision for BinoUCB:
	delta <- 0.00001
	
	a1 <- 0
	b1 <- 1
	a2 <- 0
	b2 <- 1
	if (y > 0)
	{
		a1 <- 1 - BinoUCB(n-y,n,alpha,delta)
		a2 <- 1 - BinoUCB(n-y,n,alpha/2,delta)	}
	if (y < n)
	{
		b1 <-     BinoUCB(y,n,alpha,delta)
		b2 <-     BinoUCB(y,n,alpha/2,delta)
	}
	return(list(
		observation=y,
		sample.size=n,
		test.level=alpha,
		estimator=y/n,
		lower.bound=a1, upper.bound=b1,
		confidence.interval=c(a2,b2)))
}


#====================
# Auxiliary programs:
#====================

BinoUCB <- function(x, n, alpha, delta)
# b = BinoUCB(x, n, alpha, delta) computes for 
# x in {0,1,2,...,n-1} an approximation b = b(x) for the 
# unique number p in ]0,1[ such that 
#    BinoCDF(x, n, p)  =  alpha.
# Namely, p <= b <= p + delta and
# alpha >= BinoCDF(x, n, b) >= alpha - delta.
# 
# If X is a random variable with distribution Bin(n,p), then 
#    Pr{p > b(X)}  <=  alpha .
# 
# Lutz Duembgen, August 28, 2008
{
	a <- 0
	b <- 1
	pa <- 1
	pb <- 0
	while ((b - a > delta) || (pa - pb > delta))
	{
		t <- (a+b)/2
		pt <- pbinom(q=x,size=n,prob=t)
		if (pt >= alpha)
		{
			a <- t
			pa <- pt
		}
		else
		{
			b <- t
			pb <- pt
		}
	}
	return(b)
}