Note: You may work with a partner but you must write up and submit your own assignment. Submit your Rmarkdown (.Rmd) and the compiled pdf on Gradescope in a zip file. Include any additional files (e.g. scanned handwritten solutions) in zip file with the pdf.

Problem 1. Frequentist Coverage of The Bayesian Posterior Interval.

In quiz 1 we explored the importance and difficulty of well-calibrated prior distributions by examining the calibration of subjective intervals. Suppose that \(y_1,..,y_n\) is an IID sample from a \(Normal(\mu, 1)\). We wish to estimate \(\mu\).

1a. (20 points) For Bayesian inference, we will assume the prior distribution \(\mu \sim Normal(0,\frac{1}{\kappa_0})\) for all parts below. Remember, from lecture that we can interpret \(\kappa_0\) as the pseudo-number of prior observations with sample mean \(\mu_0 = 0\). State the posterior distribution of \(\mu\) given \(y_1,..,y_n\). Report the lower and upper bounds of the \(95\%\) quantile-based posterior credible interval for \(\mu\), using the fact that for a normal distribution with standard deviation \(\sigma\), approximately \(95\%\) of the mass is between \(\pm 1.96\sigma\).

SEE PDF

1b. (20 points) Plot the length of the posterior credible interval as a function of \(\kappa_0\), for \(\kappa_0 = 1, 2, ..., 25\) assuming \(n=10\). Report how this prior parameter effects the length of the posterior interval and why this makes intuitive sense.

n <- 10
k0_values <- 1:25

# posterior variance for each k0
posterior_variances <- 1 / (k0_values + n)

# length of 95% credible intervals
credible_interval_lengths <- 2 * 1.96 * sqrt(posterior_variances)

# data frame 
data <- data.frame(k0 = k0_values, IntervalLength = credible_interval_lengths)

# plot
ggplot(data, aes(x = k0, y = IntervalLength)) +
  geom_line() +
  geom_point() +
  theme_minimal() +
  labs(title = "Length of 95% Posterior Credible Interval vs. k0",
       x = "k0 (Strength of Prior)",
       y = "Length of 95% Credible Interval")

This makes intuitive sense because larger k0 value (stronger prior) makes the sample data less influential which leads to more confidence in the central value. This explains the narrower confidence interval.

1c. (20 points) Now we will evaluate the frequentist coverage of the posterior credible interval on simulated data. Generate 1000 data sets where the true value of \(\mu=0\) and \(n=10\). For each data set, compute the posterior \(95\%\) interval endpoints (from the previous part) and see if it the interval covers the true value of \(\mu = 0\). Compute the frequentist coverage as the fraction of these 1000 posterior 95% credible intervals that contain \(\mu=0\). Do this for each value of \(\kappa_0 = 1, 2, ..., 25\). Plot the coverage as a function of \(\kappa_0\).

set.seed(101)
n <- 10
k0_values <- 1:25
num_simulations <- 1000

# simulate data and compute coverage
get_coverage <- function(k0) {
  coverage_count <- 0
  
  for (i in 1:num_simulations) {
    # simulate data
    sample_data <- rnorm(n, mean = 0, sd = 1)
    
    # posterior mean and variance
    mu_n <- sum(sample_data) / (k0 + n)
    sigma_n_squared <- 1 / (k0 + n)
    
    # 95% credible interval
    lower <- mu_n - 1.96 * sqrt(sigma_n_squared)
    upper <- mu_n + 1.96 * sqrt(sigma_n_squared)
    
    # check if the true mu = 0 is within the interval
    if (lower <= 0 && upper >= 0) {
      coverage_count <- coverage_count + 1
    }
  }
  
  return(coverage_count / num_simulations)
}

# coverage for each k0
coverage_rates <- sapply(k0_values, get_coverage)

# data frame
coverage_data <- data.frame(k0 = k0_values, Coverage = coverage_rates)

# plot
ggplot(coverage_data, aes(x = k0, y = Coverage)) +
  geom_line(color = "blue") +
  geom_point(color = "red") +
  theme_minimal() +
  labs(title = "Frequentist Coverage",
       x = "k0 (Strength of Prior)",
       y = "Coverage Rate")

1d. (20 points) Repeat the process as in 1c, but now generate data assuming the true \(\mu=1\).

set.seed(101)

n <- 10
k0_values <- 1:25
num_simulations <- 1000
true_mu <- 1  

# simulate data and compute coverage
get_coverage <- function(k0) {
  coverage_count <- 0
  
  for (i in 1:num_simulations) {
    # simulate data
    sample_data <- rnorm(n, mean = true_mu, sd = 1) 
    
    # posterior mean and variance
    mu_n <- sum(sample_data) / (k0 + n)
    sigma_n_squared <- 1 / (k0 + n)
    
    # compute 95% credible interval
    lower <- mu_n - 1.96 * sqrt(sigma_n_squared)
    upper <- mu_n + 1.96 * sqrt(sigma_n_squared)
    
    # check if the true mu = 1 is within the interval
    if (lower <= true_mu && upper >= true_mu) {
      coverage_count <- coverage_count + 1
    }
  }
  
  return(coverage_count / num_simulations)
}

# coverage for each k0
coverage_rates <- sapply(k0_values, get_coverage)

# data frame
coverage_data <- data.frame(k0 = k0_values, Coverage = coverage_rates)

# plot
ggplot(coverage_data, aes(x = k0, y = Coverage)) +
  geom_line(color = "blue") +
  geom_point(color = "red") +
  theme_minimal() +
  labs(title = "Frequentist Coverage With True μ",
       x = "k0 (Strength of Prior)",
       y = "Coverage Rate")

1e. (20 points) Explain the differences between the coverage plots when the true \(\mu\) = 0 and the true \(\mu = 1\). For what values of \(\kappa_0\) do you see closer to nominal coverage (i.e. 95%)? For what values does your posterior interval tend to over-cover (the interval covers the true value more than 95% of the time)? Under-cover (the interval covers the true value less than 95% of the time)? Why does this make sense?

When \(\mu = 0\) the prior was centered around the around the true value of the parameter. This means that as the strength of the prior increases, the frequency rate should increase as well, as shown in the first graph. When \(\mu = 1\), the opposite happens. The prior is not centered around the true value of the parameter. This misalignment causes the coverage rate to decrease as the strength of the prior increases.

In both scenarios, a weaker prior (smaller k0) leads to nominal coverage. This is because the weak prior allows the sample data to dominate the posterior which leads to a 95% coverage rate.

The first posterior distribution tended to over-cover for larger values of k0. This is because the prior for this example was well aligned with the truth, so a more influential strength lead to a higher coverage rate of the true value. The first posterior distribution did not under-cover for any values of k0. The second posterior distribution did not over-cover for any values of k0, but it under-covered when the value of k0 increased. This is because the prior for this example was not well aligned with the true value of \(\mu\). As the influence of the prior increased, the posterior distribution moved closer and closer to the un-true value of \(\mu\), which led to under-coverage.