R for Statistical Inference Question 6(b)

set.seed(123)

# Function to calculate the coverage probability of the credible interval 
coverage_probability <- function(n, p, alpha, num_simulations) {
  z_alpha <- qnorm(1 - alpha / 2)
  covered <- 0
  
  for (i in 1:num_simulations) {
    # Simulate binomial data 
    x <- rbinom(1, n, p)
    
    # Compute the quantiles for the credible interval
    # Use non-informative priors (1, 1)
    lower_bound <- qbeta(alpha/2, x+1, n-x+1)
    upper_bound <- qbeta(1-alpha/2, x+1, n-x+1)

    # 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 Credible Interval for Binomial Data")
abline(h = 1 - alpha, col = "red", lty = 2)