R语言模拟直方图叠加F分布：ggplot或循环问题排查

下面是代码里的几个关键问题，以及修正后的版本：

• 核心参数未提前定义：你没给Nsim（要求是10000）、sigmax、sigma、b0、b1赋值，R要么报错，要么乱用环境中存在的其他值，直接导致模拟结果失真。按照要求需添加：

  Nsim <- 10000
  b0 <- 0
  b1 <- 0
  sigmax <- 1
  sigma <- 1
  

• 变量名重复冲突：循环里先将残差赋值给res，随后又把模型汇总结果覆盖到res上，虽然本次未影响F统计量提取，但极易引发逻辑混乱，建议将残差变量改名为epsilon。

• 理论分布曲线生成粗糙：

  • 你用seq(1,10,1)生成x轴序列，步长1过大，画出的曲线呈锯齿状，无法准确匹配直方图的分布形态，应改为seq(0,10,0.1)缩小步长。
  • x序列起始值设为1，遗漏了F分布0到1的区间，而模拟得到的F统计量大量落在该区间，导致理论曲线与直方图完全错位。

• 绘图范围未对齐：直方图的x范围需与理论曲线一致，从0开始才能直观对比分布拟合效果。

修正后的完整代码：

library(tibble)
library(ggplot2)

# 定义所有模拟参数
n <- 30
Nsim <- 10000
b0 <- 0
b1 <- 0
sigmax <- 1
sigma <- 1

F1 <- array(NA, dim = Nsim)
for(i in 1:Nsim){
  X <- rnorm(n, 0, sd = sigmax) 
  epsilon <- rnorm(n, 0, sd = sigma)  # 更换变量名避免冲突
  Y <- b0 + b1*X + epsilon 
  mod <- lm(Y~X)
  mod_summary <- summary(mod)  # 重命名汇总结果，逻辑更清晰
  F1[i] <- mod_summary$fstatistic[1] 
}

# 整理模拟数据与理论分布数据
df <- tibble(F1 = F1)
x_seq <- seq(0, 10, 0.1)
y_density <- df(x_seq, df1 = 1, df2 = n-2)  # df()是R中F分布的密度函数
df_theory <- tibble(x = x_seq, y = y_density)

# 绘图对比
ggplot() +  
  geom_histogram(data = df, aes(x = F1, y = ..density..), 
                 binwidth = 0.1, color = "black", fill = "white") + 
  geom_line(data = df_theory, aes(x = x, y = y), color = "red", linewidth = 1) +
  xlab("F统计量") +
  xlim(c(0, 10)) +
  ggtitle("n=30时F统计量分布 vs F(1,28)理论分布") +
  theme_minimal()

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