39
$\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$
49
$\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
4
$\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 of proportions in the interval lying beyond 0 or 1, and it 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 acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.