R语言循环构建置信区间仅最后行有值问题排查求助

设置n=100、B=20、S=50，循环生成Gamma分布原始样本，抽取Bootstrap样本构建三种置信区间（渐近正态CI、渐近Bootstrap CI、分位数Bootstrap CI），但最终仅置信区间矩阵的最后一行有有效数据，需修正代码实现每个原始样本对应一行的结果。

原代码 #

n <- 100
B <- 20
S <- 50
alpha <- 0.3
beta <- 1.2
theta <- alpha*beta

for (i in 1:S) {
  ###
  sim_original_samples <- rgamma(n, alpha, beta)                      # for each S, we have a sample of 100 observations
  sim_original_samples_X_bar <- mean(sim_original_samples)            # for each dataset, compute the sample mean and input it 
  sim_bs_samples_X_bar <- matrix(0,B,1)
  # in the same loop we are going to compute the sample mean per bootstrap per original sample i
  ####
  ####
  for (j in 1:B) {
    sim_bs_samples <- sample(sim_original_samples,n,replace=TRUE)
    # for each original sample, we are going to draw B times a bootstrap sample
    sim_bs_samples_X_bar[j] <- mean(sim_bs_samples)
    # all the elements of this matrix should be the bootstrap sample mean
    var_sim_bs_samples <- matrix(0,B,1)
    var_sim_bs_samples[j] <- (sim_bs_samples_X_bar[j] - sim_original_samples_X_bar)^2
    se_sim_bs_samples <- sqrt((1/B*sum(var_sim_bs_samples)))
  }
  ####
  ####
  # now we want to compute the asymptotic CI of i)
  z <- 1.96                              
  var_gamma <- alpha*beta^2/n      
  CI_sim_asy_norm <- matrix(ncol = 3, nrow = S)         # create a vector for the CI
  names <- c("Lower bound", "Upper bound", "teta covered")
  colnames(CI_sim_asy_norm) <- names
  #
  CI_sim_asy_norm[i,1] <- theta - z*sqrt(var_gamma)
  CI_sim_asy_norm[i,2] <- theta + z*sqrt(var_gamma)
  CI_sim_asy_norm[i,3] <- theta >= CI_sim_asy_norm[i,1] & theta <= CI_sim_asy_norm[i,2]
  # check whether the true parameter of interest is covered 
  ####
  ####
  # do the same for the asymptotic BS CI of ii)
  CI_sim_asy_bs <- matrix(ncol = 3, nrow = S)
  colnames(CI_sim_asy_bs) <- names
  CI_sim_asy_bs[i,1] <- sim_original_samples_X_bar - z*se_sim_bs_samples
  CI_sim_asy_bs[i,2] <- sim_original_samples_X_bar + z*se_sim_bs_samples
  CI_sim_asy_bs[i,3] <- theta >= CI_sim_asy_bs[i,1] & theta <= CI_sim_asy_bs[i,2]
  ####
  ####
  # do the same for the percentile BS CI of iii) assuming B = 1000 for simplicity
  sim_bs_samples_X_bar_sorted <- sort(sim_bs_samples_X_bar, decreasing=FALSE)
  CI_sim_percentile <- matrix(ncol = 3, nrow = S)
  colnames(CI_sim_percentile) <- names
  CI_sim_percentile[i,1] <-  sim_bs_samples_X_bar_sorted[1000*(0.05/2)]
  CI_sim_percentile[i,2] <-  sim_bs_samples_X_bar_sorted[1000*((1-0.05)/2)]
  CI_sim_percentile[i,3] <-  theta >= CI_sim_percentile[i,1] & theta <= CI_sim_percentile[i,2]
  ####
}

问题原因分析 #

• 置信区间矩阵重复初始化：三个结果矩阵CI_sim_asy_norm、CI_sim_asy_bs、CI_sim_percentile被放在for (i in 1:S)循环内部创建，每次循环都会生成新的空矩阵，覆盖之前循环填充的数据，最终仅保留最后一次循环的结果。
• Bootstrap标准误计算逻辑错误：var_sim_bs_samples在for (j in 1:B)循环内部重复初始化，导致每次仅能计算当前j对应的方差项，无法累积所有Bootstrap样本的方差，标准误计算结果不准确。
• 分位数Bootstrap索引错误：代码中错误使用1000*(0.05/2)作为索引，但实际设置的Bootstrap次数B=20，索引超出范围，会导致取值为NA或错误值。

修正后的代码 #

n <- 100
B <- 20
S <- 50
alpha <- 0.3
beta <- 1.2
theta <- alpha*beta
z <- 1.96
var_gamma <- alpha*beta^2/n      
names <- c("Lower bound", "Upper bound", "theta covered")

# 提前初始化所有结果矩阵，放在循环外部
CI_sim_asy_norm <- matrix(ncol = 3, nrow = S)
colnames(CI_sim_asy_norm) <- names

CI_sim_asy_bs <- matrix(ncol = 3, nrow = S)
colnames(CI_sim_asy_bs) <- names

CI_sim_percentile <- matrix(ncol = 3, nrow = S)
colnames(CI_sim_percentile) <- names

for (i in 1:S) {
  # 生成当前原始样本并计算均值
  sim_original_samples <- rgamma(n, alpha, beta)
  sim_original_X_bar <- mean(sim_original_samples)
  
  # 存储Bootstrap样本均值的向量，提前初始化
  sim_bs_X_bar <- numeric(B)
  
  # 计算Bootstrap样本均值，累积方差项
  var_bs <- numeric(B)
  for (j in 1:B) {
    sim_bs_sample <- sample(sim_original_samples, n, replace = TRUE)
    sim_bs_X_bar[j] <- mean(sim_bs_sample)
    var_bs[j] <- (sim_bs_X_bar[j] - sim_original_X_bar)^2
  }
  # 计算Bootstrap标准误
  se_bs <- sqrt(mean(var_bs))
  
  # 填充渐近正态置信区间
  CI_sim_asy_norm[i, 1] <- theta - z * sqrt(var_gamma)
  CI_sim_asy_norm[i, 2] <- theta + z * sqrt(var_gamma)
  CI_sim_asy_norm[i, 3] <- theta >= CI_sim_asy_norm[i,1] & theta <= CI_sim_asy_norm[i,2]
  
  # 填充渐近Bootstrap置信区间
  CI_sim_asy_bs[i, 1] <- sim_original_X_bar - z * se_bs
  CI_sim_asy_bs[i, 2] <- sim_original_X_bar + z * se_bs
  CI_sim_asy_bs[i, 3] <- theta >= CI_sim_asy_bs[i,1] & theta <= CI_sim_asy_bs[i,2]
  
  # 填充分位数Bootstrap置信区间，用quantile函数更可靠
  bs_percentiles <- quantile(sim_bs_X_bar, c(0.025, 0.975))
  CI_sim_percentile[i, 1] <- bs_percentiles[1]
  CI_sim_percentile[i, 2] <- bs_percentiles[2]
  CI_sim_percentile[i, 3] <- theta >= CI_sim_percentile[i,1] & theta <= CI_sim_percentile[i,2]
}

# 查看结果示例
head(CI_sim_asy_norm)
head(CI_sim_asy_bs)
head(CI_sim_percentile)

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