bin.int <- function(x,n,conf=.95) {

#
# Descrição e detalhes:
# Função constrói os seguintes intervalos para uma binomial:
# Wald: O intervalo obtido ao inverter o teste Wald (Agresti, 2002, pág.15) tende a cobrir sempre menos que o nível de
#	confiança estabelecido;
# score: O intervalo obtido ao inverter o teste escore pelo método de Wilson (Agresti, 2002, pág.16). Oscila cobrindo as
#	 vezes um pouco menos e as vezes um pouco mais que o nível de confiança estabelecido. Quando a proporção é
#	 próxima de 0 ou 1, irá cobrir sempre menos que o nível estabelecido;
# Agresti: O intervalo sugerido por Agresti é uma correção do intervalo de Wald motivada pelos pesos do intervalo de escore.
# 	   Agresti e Coull (1988) e Brown et alli (2001) mostraram que esse intervalo se comporta muito melhor que o Wald.
# 	   Ele contém o intervalo de escore (é sempre maior) e se comporta melhor que o de escore quando a proporção está
#	   próxima de 0 ou 1, veja Agresti (2002, págs.32-33);
# Clopper: Teste exato de Clopper-Pearson (Agresti, 2002, pág.18). Tende a cobrir sempre mais que o nível de confiança
#	   estabelecido;
# Blaker: Teste exato de Blaker (2000), Agresti (2002, pág.20). Ao invés de inverter um teste bilateral, usa dois testes
#	  unilaterais. Agresti (2002, pág.27) comenta que o intervalo de Blaker tem propriedades boas e contém sempre o
#	  intervalo de Clopper-Pearson.
#
# Todos as referências citadas e não listadas abaixo podem ser encontradas em Agresti (2002).
# Sobre a implementação computacional, veja Thompson (2003).
#
# Argumentos obrigatórios:
# x: número de sucessos observados;
# n: tamanho da amostra.
#
# Argumentos opcionais:
# conf: se não for estabelecido outro nível de confiança a função usa 95%.
#
# Autor: Frederico Zanqueta Poleto <fred@poleto.com>, arquivo disponível em http://www.poleto.com
#
# Referências:
# AGRESTI, A. (2002). Categorical Data Analysis. 2ª ed. John Wiley & Sons, New Jersey.
# THOMPSON, L. A. (2003). S-PLUS (and R) Manual to Accompany Agresti’s Categorical Data Analysis (2002). 2ª ed.
#    [Não publicado, disponível em http://math.cl.uh.edu/~thompsonla/Splusdiscrete2.pdf]
#
# Exemplos:
# bin.int(20,100)
# bin.int(20,100,conf=.99)
#

library(Hmisc, T) #para função binconf
phat <- x/n
z<- - qnorm(p = (1-conf)/2)

#wald<-phat + c(-1, 1)*z*sqrt((phat*(1 - phat))/n)
#ou equivalentemente
wald<-binconf(x=x, n=n, alpha=1-conf, method="asymptotic")[2:3]

#score<-(x + (z^2)/2)/(n + (z^2)) + c(-1, 1)*z*sqrt( (phat*(1 - phat))*n + .25*(z^2) ) / (n+z^2)
#ou equivalentemente
#score<-binconf(x=x, n=n, alpha=1-conf, method="wilson")[2:3]
#ou equivalentemente
score<-prop.test(x=x,n=n,conf.level=conf,correct=F)$conf.int[1:2]

phat2<-(x + (z^2)/2)/(n + (z^2))
agresti<-phat2 + c(-1, 1)*z*sqrt((phat2*(1 - phat2))/(n+z^2))

clopper<-binconf(x=x, n=n, alpha=1-conf, method="exact")[2:3]

blaker<-acceptinterval(x=x, n=n, level=conf)

intervalos<-rbind(wald,score,agresti,clopper,blaker)
dimnames(intervalos)<-list(c("Wald","score","Agresti","Clopper","Blaker"),c("Inferior","Superior"))
return(intervalos)

}

acceptbin<-function(x, n, p){
	# computes the acceptability of p when x is observed and X is Bin(n, p)
	# adapted from Blaker (2000)
	p1<-1-pbinom(x-1, n, p)
	p2<-pbinom(x, n, p)
	a1<-p1 + pbinom(qbinom(p1, n, p) - 1, n, p)
	a2<-p2 + 1 - pbinom(qbinom(1-p2, n, p), n, p)
	min(a1, a2)
}

acceptinterval<-function(x, n, level, tolerance=1e-04){
	# computes acceptability interval for p at 1 - alpha equal to level
	# (in (0,1)) when x is an observed value of X which is Bin(n, p)
	# adapted from Blaker (2000)
	lower<-upper<-0

	if(x) {
		lower<-qbeta((1-level)/2, x, n-x+1)
		while(acceptbin(x, n, lower) < (1-level)) lower<-lower+tolerance
	}

	if(x!=n) {
		upper<-qbeta(1-(1-level)/2, x + 1, n - x)
		while(acceptbin(x, n, upper) < (1-level)) upper<-upper-tolerance
	}

#	c(lower=lower, upper=upper)
	c(lower, upper)

}