基于百分位法的Bootstrap置信区间覆盖率与宽度精度对比研究

Got it, let's break this down into a practical, reproducible R solution. You want to evaluate how well the 95% percentile Bootstrap confidence interval performs for your scaled t₅ distribution parameter—specifically, its coverage accuracy (how often it captures the true parameter) and average width. Here's how to implement this:

Step-by-Step Implementation #

First, let's clarify the true value of θ we're targeting. For your scaled t₅ distribution with mean=0 and sd=2, the expected value of your estimator thetahat can be derived using the moments of the t-distribution:

• The variance of a t₅ distribution is df/(df-2) = 5/3, so scaling it by sd=2 gives a variance of 4, meaning E[X²] = Var(X) + (E[X])² = 4.
• The true θ is sqrt( (3/500)*n*E[X²] ) = sqrt( (3/500)*100*4 ) = sqrt(2.4) ≈ 1.549.

Full R Code #

# Load required package for scaled t-distribution samples
library("metRology", lib.loc="~/R/win-library/3.4")

# Define the estimator functions
thetahatsq <- function(x){(3/500)*sum(x^2)}
thetahat <- function(x){sqrt(thetahatsq(x))}

# Calculate the true value of theta
n <- 100
true_theta <- sqrt( (3/500)*n*4 )  # Theoretical true value (~1.549)

# Number of simulation repetitions to assess coverage and width
n_sims <- 1000

# Initialize vectors to store results
coverage_results <- logical(n_sims)
interval_widths <- numeric(n_sims)

# Set seed for reproducibility
set.seed(42)

# Run the simulation loop
for (i in 1:n_sims) {
  # 1. Generate original data sample from scaled t₅
  data <- rt.scaled(n, df=5, mean=0, sd=2)
  
  # 2. Compute theta hat for the original sample
  theta_hat <- thetahat(data)
  
  # 3. Generate 1000 Bootstrap samples and compute theta hat for each
  bootstrap_thetas <- replicate(1000, thetahat(rt.scaled(n, df=5, mean=0, sd=theta_hat)))
  
  # 4. Construct 95% percentile Bootstrap confidence interval
  ci_percentile <- quantile(bootstrap_thetas, probs=c(0.025, 0.975))
  
  # 5. Check if true theta is inside the interval
  coverage_results[i] <- (true_theta >= ci_percentile[1]) & (true_theta <= ci_percentile[2])
  
  # 6. Calculate and store interval width
  interval_widths[i] <- ci_percentile[2] - ci_percentile[1]
}

# Compute summary metrics
average_coverage <- mean(coverage_results)
average_width <- mean(interval_widths)

# Print results
cat(sprintf("Average 95%% Percentile Bootstrap Coverage Rate: %.4f\n", average_coverage))
cat(sprintf("Average Interval Width: %.4f\n", average_width))

# Optional: Visualize interval width distribution
hist(interval_widths, breaks=20, col="lightblue", 
     main="Distribution of Bootstrap Interval Widths",
     xlab="Interval Width")

Key Explanations #

• Simulation Loop: We repeat the entire Bootstrap process 1000 times to get reliable estimates of coverage and width.
• Bootstrap Sample Generation: For each original sample, we generate 1000 Bootstrap samples using the scaled t₅ distribution fitted with our θ̂ estimate.
• Confidence Interval Construction: The percentile interval uses the 2.5th and 97.5th quantiles of the Bootstrap θ̂ values.
• Coverage Check: We track whether the true θ falls within each interval, then take the mean to get the overall coverage rate.
• Interval Width: We calculate the width of each interval and average them to understand the interval's precision.

内容的提问来源于stack exchange，提问作者Hans Christensen