Comparing Exact and Approximate Confidence Intervals

STAT/MATH 355: Statistical Theory

Author

Solutions

Published

April 2, 2026

Setup

Let \(X_1\),…,\(X_n\) be from a Bernoulli distribution, so \[p_X(x_i; p) = p^{x_i}(1-p)^{1-x_i} \quad \text{for} \quad x_i = 0,1 \quad \text{and} \quad 0 < p < 1\]

We will use two different approaches to construct a confidence interval for \(p\):

  • establishing an approximate pivotal quantity for \(p\)
  • using our estimator’s exact dsitribution

Approximate CI

See the course website for exercises walking you through the derivation of an approximate 95% CI for \(p\).

Write down the formula for that CI here:

\[\hat{p} \pm 1.96 \sqrt{\frac{\hat{p}(1-\hat{p})}{n}}\]

Exact CI

Since \(X_1, \dots, X_n \overset{iid}{\sim} \text{Bernoulli}(p)\), we know \(X = \sum_{i=1}^n X_i \sim \text{Binomial}(n, p)\).

The Clopper-Pearson Interval uses this result to construct an exact CI for \(p\).

Suppose we observe \(X = k\) successes. The Clopper-Pearson CI for \(p\) is constructed as follows:

  • Find the largest \(p\) such that \(P(X \leq k) \geq \frac{\alpha}{2}\). Call this largest value \(p_U\).
  • Find the smallest \(p\) such that \(P(X \geq k) \geq \frac{\alpha}{2}\). Call this smallest value \(p_L\).
  • Define the \(100(1-\alpha)\%\) confidence interval for \(p\) to be: \((p_L, p_U)\).


Exercises

Part 1

  1. Read through the CI construction process described above. Turn to a neighbor and try to articulate the logic behind this process. Why does this yield a valid CI?

Confidence intervals provide us with a range of plausible values for our unknown parameter that are compatible with the data we observed. Thus, to construct a confidence interval for \(p\) we need to think about how big we could make \(p\) and still plausibly see as few as \(k\) successes, and how small we could make \(p\) and still plausibly see as many as \(k\) successes. So, to find the lower bound of our CI (i.e., the smallest \(p\) that still looks compatible with our data): find the smallest \(p\) such that \(P(X \ge k) \ge \frac{\alpha}{2}\) (i.e., there’s still a reasonable probability of seeing as many as \(k\) successes). Similarly, for the upper bound of our CI we want to find the largest \(p\) such that \(P(X \le k) \ge \frac{\alpha}{2}\) (i.e, there’s still a reasonable probability of seeing as few as \(k\) successes).

  1. Suppose we observe 3 successes in 10 trials. Find the Clopper-Pearson 95% CI for \(p\). Use R’s pbinom function to calculate \(P(X \le k)\) = pbinom(k, size = n, prob = p) and \(P(X \ge k) = 1 - P(X < k) = 1 - P(X \le k - 1)\) = ... You can use a trial-and-error/grid search approach, or define a function that reflects equality in those inequalities and use R’s uniroot to find the root of your function. Starter code for the latter is provided below.
## to use uniroot, we need to set up a function f
## where f(p) = 0 when p is a soln to the eqn we want to solve

## use uniroot to solve for upper bound
## find p such that P(X <= k) = alpha/2 or P(X <= k) - alpha/2 = 0
f_upper <- function(p){
  pbinom(3, size = 10, prob = p) - 0.05/2
}

pU <- uniroot(f_upper, interval = c(0,1))$root ## search over [0,1] since p is probability


## now use uniroot to solve for lower bound
## find p such that P(X >= k) = alpha/2 or P(X >= k) - alpha/2 = 0
f_lower <- function(p){
  1 - pbinom(3 - 1, size = 10, prob = p) - 0.05/2 ## use pbinom(k - 1), not pbinom(k)
}

pL <- uniroot(f_lower, interval = c(0,1))$root

## print out CI
CI <- c(pL, pU)
print(CI)
[1] 0.06675044 0.65246828
  1. Check your answer to the previous question by running the binom.test function, which is a pre-existing function in R that can be used to get an exact confidence interval for \(p\).
## get exact 95% CI for p using binom.test
binom.test(x = 3, n = 10, conf.level = 0.95)

    Exact binomial test

data:  3 and 10
number of successes = 3, number of trials = 10, p-value = 0.3438
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.06673951 0.65245285
sample estimates:
probability of success 
                   0.3
  1. Compare the Clopper-Pearson CI to the approximate CI you derived earlier. Comment on any similarities or differences that you notice.
# setup
alpha <- 0.05
phat <- 3/10
n <- 10

# get quantiles of standard normal
zLower <- qnorm(alpha/2, mean = 0, sd = 1) # -1.96
zUpper <- qnorm(1 - alpha/2, mean = 0, sd = 1) # 1.96

# get lower and upper bounds of CI
pLower <- phat - zUpper * sqrt(phat*(1-phat)/n)
pUpper <- phat - zLower * sqrt(phat*(1-phat)/n)

# print out CI
CI <- c(pLower, pUpper)
print(CI)
[1] 0.01597423 0.58402577

The approximate CI is narrower than the Clopper-Pearson interval!

  1. How might you go about deciding which of these CIs to report?

We’ll talk more about this next class!


Part 2

With confidence intervals, our goal is to construct an interval that contains the true parameter as an interior point \(100(1-\alpha)\)% of the time, ie \[P(LB \le \theta \le UB) = 1-\alpha\]

It’s worth confirming whether this is actually true: when we construct a 95% CI according to some procedure, do those intervals actually contain the true parameter 95% of the time?

To check this, we need to design and implement a simulation study:

  1. generate some data where we know the true \(\theta\)
  2. construct a CI based on this data
  3. check if \(\theta\) is in this interval
  4. repeat steps 1–3 many times
  5. across all these reps, how often does the CI actually cover the true \(\theta\)?

This simulation approach yields an empirical estimate of the coverage of the CI—we’d like coverage to be close to \(100(1-\alpha)\)%.

Let’s perform a simulation study to evaluate the coverage of the exact and approximate CIs we constructed above!

  1. Insert code into the do_one_p function, below, to calculate an exact and approximate \(100(1-\alpha)\) percent confidence interval for \(p\) using your work from above.

    NOTES:

    • I’ve given you most of the code here. Take your time looking through this to make sure you understand what each line is doing. ASK QUESTIONS about any parts that don’t make sense.
    • Feel free to update as you see fit, especially if you see ways to make this code more efficient, easier to read, etc.!
## function to do one simulation replicate
do_one_p <- function(n, p, alpha = 0.05){
  # generate sample of size n from Binom(1,p)
  x <- rbinom(n, size = 1, prob = p)
  phat <- mean(x)
        
  # exact CI
  exactCI <- binom.test(x = sum(x), n = n,
                        conf.level = 1-alpha)$conf.int

  # approximate CI
  zL <- qnorm(alpha/2)
  zU <- qnorm(1 - alpha/2)
  approxCI <- c(phat - zU*sqrt(phat*(1-phat)/n), 
                phat - zL*sqrt(phat*(1-phat)/n))
        
  # check if true p is in CIs
  inExact <- (p >= exactCI[1]) & (p <= exactCI[2])
  inApprox <- (p >= approxCI[1]) & (p <= approxCI[2])
        
  # return results
  return(c(inExact, inApprox))
}

## run simulation once to check what code is doing
## (using n = 10, p = 0.3)
do_one_p(n = 10, p = 0.3)
[1] TRUE TRUE
## repeat simulation 1000 times
## (using n = 10, p = 0.3)
set.seed(455)
simres <- replicate(1000, do_one_p(n = 10, p = 0.3))

## calculate coverage for each type of CI
apply(simres, 1, mean)
[1] 0.984 0.838

Notice that the coverage of the approximate confidence interval is quite a bit lower than 95%! This is bad. Although we are saying we are “95% confident” in this CI construction process, we should not be. (We’re actually closer to 85% confident.)

It also seems like exact CI might be a bit wider than it needs to be, since our coverage is over 95%. But this is less of a concern: better to under-promise and over-deliver than the reverse!

  1. What if we increased \(n\)? Would you expect the coverage of either interval to be impacted? Write down your hypothesis, and then try running your simulation again with different values of \(n\) to test your hypothesis.

Since the approximate CI relies on the CLT, which is an asymptotic (\(n \rightarrow \infty\)) result, I would hypothesize that the coverage of the approximate CI should improve as we get larger samples. Let’s test this out with n = 100 and n = 1000:

## repeat the process above, with different values of n

## try simulation using n = 10 and 10k reps
set.seed(455)
simres <- replicate(10000, do_one_p(n = 10, p = 0.3))
apply(simres, 1, mean)
[1] 0.9884 0.8395
## try simulation using n = 100 and 10k reps
set.seed(455)
simres <- replicate(10000, do_one_p(n = 100, p = 0.3))
apply(simres, 1, mean)
[1] 0.9618 0.9486
## try simulation using n = 1000 and 10k reps
set.seed(455)
simres <- replicate(10000, do_one_p(n = 1000, p = 0.3))
apply(simres, 1, mean)
[1] 0.9544 0.9503
## loop over many values of n
set.seed(455)
ns <- seq(1, 100, by = 5)
coverages_approx <- c()
coverages_exact <- c()
for(i in 1:length(ns)){
  simres <- replicate(10000, do_one_p(n = ns[i], p = 0.3))
  coverages_exact[i] <- apply(simres, 1, mean)[1]
  coverages_approx[i] <- apply(simres, 1, mean)[2]
}
## make a plot summarizing how coverage changes with n
library(tidyverse)
plotdat <- data.frame(n = ns, 
                      exact = coverages_exact, 
                      approx = coverages_approx) %>%
  pivot_longer(cols = exact:approx, names_to = "CI Method", 
               values_to = "coverage")

plotdat %>%
  ggplot(aes(x = n, y = coverage, color = `CI Method`)) + 
  geom_point() + 
  geom_line() + 
  theme_bw() + 
  geom_hline(yintercept = 0.95, linetype = 2)

As \(n\) gets bigger, the coverage of the approximate method approaches the desired 95%!

  1. Repeat the process above, this time keeping \(n\) fixed but changing \(p\). What do you learn about how the true probability of success impacts the performance of these confidence intervals?
## loop over many values of n
set.seed(455)
ps <- seq(0, 1, length = 50)
coverages_approx <- c()
coverages_exact <- c()
for(i in 1:length(ps)){
  simres <- replicate(10000, do_one_p(n = 100, p = ps[i]))
  coverages_exact[i] <- apply(simres, 1, mean)[1]
  coverages_approx[i] <- apply(simres, 1, mean)[2]
}
## plot
plotdat <- data.frame(p = ps, 
                      exact = coverages_exact, 
                      approx = coverages_approx) %>%
  pivot_longer(cols = exact:approx, names_to = "CI Method", 
               values_to = "coverage")

plotdat %>%
  ggplot(aes(x = p, y = coverage, color = `CI Method`)) + 
  geom_point() + 
  geom_line() + 
  theme_bw() + 
  geom_hline(yintercept = 0.95, linetype = 2)

It seems like both methods perform better when \(p\) is closer to 0.5.