蒙特卡洛函数模拟Beta估计量一致性的代码问题排查

问题背景 #

目标是通过蒙特卡洛模拟展示：随着观测数量增加，OLS估计的Beta系数会逐渐趋近于真实值1。但初始函数在改变观测数（n=200和n=1000）时，未呈现预期的收敛效果，推测是代码编写错误而非统计问题。

初始代码的核心错误 #

初始函数存在两个关键问题：

• 随机误差项生成错误：e <- rnorm(1)*x 仅生成1个随机数，再广播到整个x向量。这使得每个观测的误差项是同一个随机数与对应x的乘积，并非独立的随机扰动。无论观测数如何变化，误差结构固定，Beta估计量的分布不会随n增大而收窄。
• 变量赋值位置错误：assign(var_name, betas, env=.GlobalEnv) 放在for循环内部，每次迭代都会覆盖全局环境中的目标变量。最终得到的obs_200或obs_1000仅保留最后几次迭代的结果，而非完整的iterations次模拟的Beta值。

修正后的代码 #

library(tidyverse)

func <- function(iterations, observations){
  set.seed(1001)
  
  betas <-c()
  for (i in 1:iterations) {
    c <- 1:observations
    x <- c/5
    e <- rnorm(observations)*x # 为每个观测生成独立随机误差
    y <- x - 10 + e
    #y[y < 0] <- 0 # 可选：添加截尾逻辑
    model <- lm(y ~ x)
    betas[i] <- model$coefficients[2]
  }
  # 循环结束后统一赋值到全局环境
  arg_name <- deparse(substitute(observations))
  var_name <- paste("obs", arg_name, sep="_") 
  assign(var_name, betas, env=.GlobalEnv) 
}

# 运行不同观测数的模拟
func(iterations = 2000, observations = 200)
func(iterations = 2000, observations = 400)
func(iterations = 2000, observations = 600)
func(iterations = 2000, observations = 800)
func(iterations = 2000, observations = 1000)

# 查看均值对比
mean(obs_200);mean(obs_1000)

# 整理绘图数据并可视化
plot_data <- tibble(a = obs_200,
       b = obs_400,
       c = obs_600,
       d = obs_800,
       e = obs_1000) |> 
  pivot_longer(cols = everything(),
               names_to = "n")

ggplot(plot_data) +
  geom_boxplot(aes(x = n,
                   y = value,
                   fill = n)) +
  scale_fill_grey(start = 0.5, labels = c("200", "400", "600", "800", "1000")) +
  theme_bw() +
  labs(x = NULL, 
       y = NULL, 
       fill = "Observations",
       title = "Monte Carlo simulation using OLS",
       subtitle = "(Heterosedastic errors)")

模拟结果 #

从图中可见，随着观测数量增加，Beta估计量的分布逐渐收窄，均值趋近于真实值1，符合预期的统计收敛特性。

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