R语言线性模型蒙特卡洛模拟循环实现技术求助

首先得指出你之前代码里的一个关键问题：你在lm()公式里重新生成了x变量，这会导致y中的x和回归用的x不是同一组数据，这会让回归结果完全错误！所以第一步要先在每次循环里统一生成所有自变量和误差项，再计算y，然后用这些变量做回归。

接下来是具体的实现步骤：

步骤1：初始化结果存储容器 #

我们需要存储10000次模拟中每个回归系数（β₀、β₁、β₂、β₃）的估计值，所以可以创建一个10000行×4列的矩阵，列名对应四个系数：

# 设置模拟次数
n_sim <- 10000
# 初始化存储系数的矩阵
coef_results <- matrix(NA, nrow = n_sim, ncol = 4)
colnames(coef_results) <- c("beta0_hat", "beta1_hat", "beta2_hat", "beta3_hat")

步骤2：编写循环执行模拟 #

用for循环遍历10000次，每次循环完成以下操作：

1. 生成样本量为50的自变量x1、x2、x3
2. 生成误差项ε
3. 计算因变量y
4. 拟合线性回归模型
5. 提取回归系数存入结果矩阵

代码示例：

# 设置真实参数
true_beta <- c(150, 1, 4, -70)
sample_size <- 50

set.seed(123) # 设置随机种子保证结果可重复
for (i in 1:n_sim) {
  # 生成自变量：注意rnorm的sd是标准差，题目给的是方差，所以要开平方
  x1 <- rnorm(sample_size, mean = 400, sd = sqrt(250))
  x2 <- rnorm(sample_size, mean = 5, sd = sqrt(6))
  x3 <- rbinom(sample_size, size = 1, prob = 0.6)
  # 生成误差项
  eps <- rnorm(sample_size, mean = 0, sd = sqrt(250))
  # 计算y
  y <- true_beta[1] + true_beta[2]*x1 + true_beta[3]*x2 + true_beta[4]*x3 + eps
  # 拟合回归
  model <- lm(y ~ x1 + x2 + x3)
  # 提取系数存入矩阵
  coef_results[i, ] <- coef(model)
}

步骤3：计算偏差、方差与MSE #

现在我们有了所有系数的估计值，接下来计算每个系数的三个核心指标：

• 偏差：mean(估计值) - 真实值
• 方差：var(估计值)
• MSE：mean((估计值 - 真实值)^2)（等价于偏差平方 + 方差）

代码示例：

# 逐个计算参数的指标
beta0_bias <- mean(coef_results[, "beta0_hat"]) - true_beta[1]
beta0_var <- var(coef_results[, "beta0_hat"])
beta0_mse <- mean((coef_results[, "beta0_hat"] - true_beta[1])^2)

beta1_bias <- mean(coef_results[, "beta1_hat"]) - true_beta[2]
beta1_var <- var(coef_results[, "beta1_hat"])
beta1_mse <- mean((coef_results[, "beta1_hat"] - true_beta[2])^2)

beta2_bias <- mean(coef_results[, "beta2_hat"]) - true_beta[3]
beta2_var <- var(coef_results[, "beta2_hat"])
beta2_mse <- mean((coef_results[, "beta2_hat"] - true_beta[3])^2)

beta3_bias <- mean(coef_results[, "beta3_hat"]) - true_beta[4]
beta3_var <- var(coef_results[, "beta3_hat"])
beta3_mse <- mean((coef_results[, "beta3_hat"] - true_beta[4])^2)

# 整理成数据框方便查看
results_df <- data.frame(
  参数 = c("β₀", "β₁", "β₂", "β₃"),
  真实值 = true_beta,
  估计均值 = round(c(mean(coef_results[,1]), mean(coef_results[,2]), mean(coef_results[,3]), mean(coef_results[,4])), 3),
  偏差 = round(c(beta0_bias, beta1_bias, beta2_bias, beta3_bias), 3),
  方差 = round(c(beta0_var, beta1_var, beta2_var, beta3_var), 3),
  MSE = round(c(beta0_mse, beta1_mse, beta2_mse, beta3_mse), 3)
)

步骤4：用表格呈现结果 #

运行上述代码后，你可以得到类似如下的markdown表格（数值为模拟近似结果）：

优化建议（可选） #

如果觉得for循环速度慢，可以用purrr包的map_dfr()函数实现向量化操作，代码更简洁且效率更高：

library(purrr)

set.seed(123)
simulate_one <- function() {
  x1 <- rnorm(sample_size, 400, sqrt(250))
  x2 <- rnorm(sample_size, 5, sqrt(6))
  x3 <- rbinom(sample_size, 1, 0.6)
  eps <- rnorm(sample_size, 0, sqrt(250))
  y <- true_beta[1] + true_beta[2]*x1 + true_beta[3]*x2 + true_beta[4]*x3 + eps
  coef(lm(y ~ x1 + x2 + x3))
}

# 生成所有模拟结果
coef_results <- map_dfr(1:n_sim, ~simulate_one())

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