如何在R语言中计算极大似然估计量的蒙特卡罗均方误差？

核心思路 #

要计算极大似然估计量（MLE）的蒙特卡罗均方误差，核心逻辑是重复生成样本→拟合模型得到MLE→计算估计值与真实参数的平方误差均值。以下是针对Gumbel分布（真实参数mu=0，sigma=1，样本量n=20，重复次数N=30）的具体实现方案：

1. 依赖包与基础设置 #

首先加载VGAM包（用于生成Gumbel样本和拟合分布），并设置随机种子保证结果可复现：

library(VGAM)
set.seed(101) # 对应首次重复的种子要求

2. 定义单次MLE计算函数 #

封装单次实验的逻辑：生成样本→拟合Gumbel分布→提取MLE参数：

get_gumbel_mle <- function(n = 20, mu_true = 0, sigma_true = 1) {
  # 从真实分布抽取样本
  x <- rgumbel(n, loc = mu_true, scale = sigma_true)
  # 拟合Gumbel分布，关闭拟合过程输出
  fit <- vglm(x ~ 1, gumbelff, trace = FALSE)
  # 提取位置参数（mu）和尺度参数（sigma）的MLE
  mu_hat <- coef(fit)[1]
  sigma_hat <- Coef(fit)[2] # 注意用Coef()提取尺度参数
  return(c(mu_hat = mu_hat, sigma_hat = sigma_hat))
}

如果需要自定义似然函数（不依赖VGAM的拟合工具），可以用optim优化实现：

# 自定义Gumbel负对数似然函数（供optim最小化）
gumbel_neg_loglik <- function(params, x) {
  mu <- params[1]
  sigma <- params[2]
  ll <- sum( -log(sigma) - (x - mu)/sigma - exp( -(x - mu)/sigma ) )
  return(-ll)
}

get_gumbel_mle_custom <- function(n = 20, mu_true = 0, sigma_true = 1) {
  x <- rgumbel(n, loc = mu_true, scale = sigma_true)
  # 设置初始参数值（用样本均值和标准差初始化）
  start_params <- c(mu = mean(x), sigma = sd(x))
  # 约束sigma为正，用L-BFGS-B方法优化
  fit <- optim(start_params, gumbel_neg_loglik, x = x, method = "L-BFGS-B", lower = c(-Inf, 1e-6))
  return(c(mu_hat = fit$par[1], sigma_hat = fit$par[2]))
}

3. 执行蒙特卡罗重复实验 #

用replicate函数批量执行N次实验，结果会返回一个每行对应单次实验的参数估计矩阵：

N <- 30
# 使用VGAM拟合的版本
mle_results <- replicate(N, get_gumbel_mle())
# 转置为数据框，方便后续计算
mle_df <- as.data.frame(t(mle_results))

# 如果用自定义似然版本，替换为：
# mle_results <- replicate(N, get_gumbel_mle_custom())
# mle_df <- as.data.frame(t(mle_results))

4. 计算蒙特卡罗均方误差 #

均方误差（MSE）的计算公式为：MSE = E[(估计值 - 真实值)^2]，用样本均值近似期望：

# 计算位置参数mu的MSE
mu_mse <- mean( (mle_df$mu_hat - 0)^2 )
# 计算尺度参数sigma的MSE
sigma_mse <- mean( (mle_df$sigma_hat - 1)^2 )

# 输出结果
cat("位置参数mu的蒙特卡罗均方误差：", round(mu_mse, 4), "\n")
cat("尺度参数sigma的蒙特卡罗均方误差：", round(sigma_mse, 4), "\n")

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