R for Statistical Inference Question 6(a)

set.seed(99528)

# Function to calculate the coverage probability of the confidence interval 
coverage_probability <- function(n, p, alpha = 0.05, num_simulations = 1000) {
  z_alpha <- qnorm(1 - alpha / 2)
  covered <- 0
  
  for (i in 1:num_simulations) {
    # Simulate binomial data 
    x <- rbinom(1, n, p)
    p_hat <- x / n
    
    # Calculate the standard error
    se <- sqrt(p_hat * (1 - p_hat) / n)
    
    # Calculate the confidence interval
    lower_bound <- p_hat - z_alpha * se
    upper_bound <- p_hat + z_alpha * se
    
    # Check if the true probability is within the interval
    if (p >= lower_bound && p <= upper_bound) {
      covered <- covered + 1
    }
  }
  
  return(covered / num_simulations)
}

# Set parameters
p <- 0.8
alpha <- 0.05
num_simulations <- 1000

# Calculate coverage probabilities for different sample sizes
sample_sizes <- seq(5, 100, by = 5)
coverage_probs <- sapply(sample_sizes, coverage_probability, p = p, alpha = alpha, num_simulations = num_simulations)

# Plot the coverage probability against sample size
plot(sample_sizes, coverage_probs, type = "l", col = "blue", xlab = "Sample Size", ylab = "Coverage Probability",
     main = "Coverage Probability of Confidence Interval for Binomial Data")
abline(h = 1 - alpha, col = "red", lty = 2)