# Confidence interval for Bernoulli sampling

I have a random sample of Bernoulli random variables $X_1 ... X_N$, where $X_i$ are i.i.d. r.v. and $P(X_i = 1) = p$, and $p$ is an unknown parameter.

Obviously, one can find an estimate for $p$: $\hat{p}:=(X_1+\dots+X_N)/N$.

My question is how can I build a confidence interval for $p$?

• If the average, $\hat{p}$, is not near $1$ or $0$, and sample size $n$ is sufficiently large (i.e. $n\hat{p}>5$ and $n(1-\hat{p})>5$, the confidence interval can be estimated by a normal distribution and the confidence interval constructed thus:

$$\hat{p}\pm z_{1-\alpha/2}\sqrt{\frac{\hat{p}(1-\hat{p})}{n}}$$

• If $\hat{p} = 0$ and $n>30$, the $95\%$ confidence interval is approximately $[0,\frac{3}{n}]$ (Javanovic and Levy, 1997); the opposite holds for $\hat{p}=1$. The reference also discusses using using $n+1$ and $n+b$ (the later to incorporate prior information).

• Else Wikipedia provides a good overview and points to Agresti and Couli (1998) and Ross (2003) for details about the use of estimates other than the normal approximation, the Wilson score, Clopper-Pearson, or Agresti-Coull intervals. These can be more accurate when above assumptions about $n$ and $\hat{p}$ are not met.

R provides functions binconf {Hmisc} and binom.confint {binom} which can be used in the following manner:

set.seed(0)
p <- runif(1,0,1)
X <- sample(c(0,1), size = 100, replace = TRUE, prob = c(1-p, p))
library(Hmisc)
binconf(sum(X), length(X), alpha = 0.05, method = 'all')
library(binom)
binom.confint(sum(X), length(X), conf.level = 0.95, method = 'all')


Agresti, Alan; Coull, Brent A. (1998). "Approximate is better than 'exact' for interval estimation of binomial proportions". The American Statistician 52: 119–126.

Jovanovic, B. D. and P. S. Levy, 1997. A Look at the Rule of Three. The American Statistician Vol. 51, No. 2, pp. 137-139

Ross, T. D. (2003). "Accurate confidence intervals for binomial proportion and Poisson rate estimation". Computers in Biology and Medicine 33: 509–531.

• (+1) Nice answer. This will become a reference for similar questions in the future, I think. However, cross-posting is unusual; in fact, I believe it is frowned on, because it screws up many aspects of the feedback/referencing/threading/commenting system. Please consider removing one of the copies and replacing it by a link in a comment. – whuber Jan 12 '11 at 2:49
• @whuber thanks for the feedback. I have removed the other copy. – David LeBauer Jan 12 '11 at 2:53
• In the first formula, what are z1 and alpha? – Cirdec Nov 22 '13 at 17:44
• I found the answer to my own question: $z_{1-\alpha/2}$ is the ${1-\alpha/2}$ percentile of the standard normal distribution and $\alpha$ is the error percentile. en.wikipedia.org/wiki/Binomial_proportion_confidence_interval – Cirdec Nov 22 '13 at 17:50
• Should that be $3/n$ on the confidence interval for the second bullet point? – Juan A. Navarro Dec 3 '13 at 14:52

## Maximum likelihood confidence intervals

The normal approximation to the Bernoulli sample relies on having a relatively large sample size and sample proportions far from the tails. The maximum likelihood estimate focuses on the log-transformed odds and this provides non-symmetric, efficient intervals for $$p$$ that should be used instead.

Define the log-odds as $$\hat{\beta}_0 = \log(\hat{p}/(1-\hat{p}))$$

A 1-$$\alpha$$ CI for $$\beta_0$$ is given by:

$$\text{CI}(\beta_0)_\alpha = \hat{\beta}_0 \pm \mathcal{Z}_{\alpha/2} \sqrt{1/(n\hat{p}(1-\hat{p})}$$

And this is back transformed into a (non-symmetric) interval for $$p$$ with:

$$\text{CI}(p)_\alpha = 1/(1+\exp(-\text{CI}(\beta_0)_\alpha)$$

This CI has the added benefit that proportions lie in the interval between 0 or 1, and the CI is always narrower than the normal interval while being of the correct level. You can get this very easily in R by specifying:

set.seed(123)
y <- rbinom(100, 1, 0.35)
plogis(confint(glm(y ~ 1, family=binomial)))

2.5 %    97.5 %
0.2795322 0.4670450


## Exact binomial confidence intervals

In small samples, the normal approximation to the MLE--while better than the normal approximation to the sample proportion--may not be reliable. That is okay. $$Y = n\hat{p}$$ can be taken to follow a binomial$$(n,p)$$ density. Bounds for $$\hat{p}$$ can be found taking the 2.5th and 97.5-th percentiles from this distribution.

$$\text{CI}_\alpha = (F^{-1}_{\hat{p}}(0.025), F^{-1}_{\hat{p}}(0.975))$$

Rarely possible by-hand, an exact binomial confidence interval can be obtained for $$p$$ using computational methods.

qbinom(p = c(0.025, 0.975), size = length(y), prob = mean(y))/length(y)
[1] 0.28 0.47


## Median unbiased confidence intervals

And if $$p$$ is 0 or 1 exactly, a median unbiased estimator can be used to obtain non-singular interval estimates based on the median unbiased probability function. You can trivially take the lower bound of the all-0 case as 0 WLOG. The upper bound is any proportion $$p_{1-\alpha/2}$$ that satisfies:

$$p_{1-\alpha/2} : P(Y = 0)/2 + P(Y > y) > 0.975$$

This is also a computational routine.

set.seed(12345)
y <- rbinom(100, 1, 0.01) ## all 0
cil <- 0
mupfun <- function(p) {
0.5*dbinom(0, 100, p) +
pbinom(1, 100, p, lower.tail = F) -
0.975
} ## for y=0 successes out of n=100 trials
ciu <- uniroot(mupfun, c(0, 1))\$root
c(cil, ciu)

[1] 0.00000000 0.05357998 ## includes the 0.01 actual probability


The last two methods are implemented in the epitools package in R.