41
$\begingroup$

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$?

$\endgroup$
50
$\begingroup$
  • 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.

$\endgroup$
  • 3
    $\begingroup$ (+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. $\endgroup$ – whuber Jan 12 '11 at 2:49
  • $\begingroup$ @whuber thanks for the feedback. I have removed the other copy. $\endgroup$ – David LeBauer Jan 12 '11 at 2:53
  • $\begingroup$ In the first formula, what are z1 and alpha? $\endgroup$ – Cirdec Nov 22 '13 at 17:44
  • $\begingroup$ 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 $\endgroup$ – Cirdec Nov 22 '13 at 17:50
  • $\begingroup$ Should that be $3/n$ on the confidence interval for the second bullet point? $\endgroup$ – Juan A. Navarro Dec 3 '13 at 14:52
5
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy