计算PRESS统计量遇阻，寻求R代码调试与完成协助

原代码存在的问题 #

• 变量初始化错误：press <- Null应为press <- numeric(r)，且混用大小写变量（press/Press）
• 误差项生成逻辑错误：未按题目要求用伯努利变量选择误差分布，直接将两个正态分布相加不符合模型定义
• 留一法实现错误：未排除目标观测就拟合模型，且混淆了r次重复与n个观测的循环层级
• 语法错误：predict中data.frame[i,]应为data[j,]，调用函数时未正确设置随机种子位置
• 直方图代码错误：使用未定义的press变量，颜色参数rainbow(press)不符合要求

修正后的完整代码 #

# 定义PRESS统计量计算函数
pressftn <- function(r, n, pi) {
  # 初始化存储r次PRESS结果的向量
  press_results <- numeric(r)
  
  for (rep in 1:r) {
    # (i) 生成符合模型要求的数据
    x1i <- runif(n, 0.5, 1.5)
    x2i <- runif(n, 0.5, 1.5)
    # 生成伯努利变量δ_i
    delta_i <- rbinom(n, 1, pi)
    # 根据δ_i生成混合误差项e_i
    ei <- ifelse(delta_i == 1, rnorm(n, 0, 0.5), rnorm(n, 0, 2))
    # 生成响应变量y_i
    yi <- 0.5 + 2*x1i - 3*x2i + ei
    data <- data.frame(x1 = x1i, x2 = x2i, y = yi)
    
    # (ii)-(iii) 执行留一法计算单次重复的PRESS值
    press_single <- 0
    for (j in 1:n) {
      # 排除第j个观测拟合回归模型
      model <- lm(y ~ x1 + x2, data = data[-j, ])
      # 预测第j个观测的y值
      pred_y <- predict(model, newdata = data[j, ])
      # 累加预测误差的平方
      press_single <- press_single + (data$y[j] - pred_y)^2
    }
    
    # 存储本次重复的PRESS结果
    press_results[rep] <- press_single
  }
  
  return(press_results)
}

# 设置随机种子保证结果可复现
set.seed(32816)
# 运行函数：50次重复，每次30个观测，伯努利参数0.5
results <- pressftn(50, 30, 0.5)
# 将结果写入文本文件
write.table(results, file="press.txt", sep=" ", row.names = FALSE, col.names = FALSE, quote=FALSE)

# (v) 绘制PRESS值的直方图
hist(results, 
     breaks = 10, 
     probability = TRUE, 
     col = rainbow(10), 
     xlab = "PRESS统计量", 
     ylab = "比例", 
     main = "50次重复计算的PRESS统计量直方图")

关键修正说明 #

1. 误差模型修复：用rbinom生成伯努利变量，通过ifelse选择对应正态误差，匹配题目中的混合误差要求
2. 留一法逻辑完善：新增内层循环遍历所有观测，每次排除当前观测拟合模型，累加误差平方得到单次PRESS值
3. 变量与语法规范：统一变量命名，修正索引错误，将随机种子设置在函数调用前保证可复现
4. 直方图优化：使用函数返回的有效变量，设置对应10个区间的彩虹配色，完善标签与标题

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