使用R语言分步实现矩估计法推导Gamma分布参数

问题背景 #

Gamma分布的概率密度函数为：
$$g(x) = \frac{1}{b^a \Gamma(a)} x^{a-1} e^{-x/b}$$
注：这里的参数定义与R中rgamma函数对应，shape=a，scale=b；当shape=1时，Gamma分布退化为指数分布。

你当前的蒙特卡洛模拟代码存在逻辑问题：循环内每次都会覆盖gamma1到gamma4变量，最终仅保留第1000次循环的样本，无法获取1000次模拟的完整结果。

接下来将手动推导矩估计步骤，并给出R语言实现代码，无需依赖gmm包。

矩估计推导步骤 #

Gamma分布的前两阶总体矩为：

• 一阶原点矩（均值）：$\mu_1' = a \cdot b$
• 二阶中心矩（方差）：$\mu_2 = a \cdot b^2$

用样本矩替代总体矩：

• 样本均值：$\bar{x} = \frac{1}{n}\sum_{i=1}^n x_i$
• 样本方差：$s^2 = \frac{1}{n-1}\sum_{i=1}^n (x_i - \bar{x})^2$

联立方程求解参数：

1. 由均值公式可得：$b = \frac{\bar{x}}{a}$
2. 将$b$代入方差公式：$s^2 = a \cdot (\frac{\bar{x}}{a})^2 = \frac{\bar{x}^2}{a}$
3. 解得参数估计值：
   • $\hat{a} = \frac{\bar{x}^2}{s2}$
   • $\hat{b} = \frac{s^2}{\bar{x}}$

R语言手动实现代码 #

第一步：修正蒙特卡洛模拟代码，存储所有模拟结果 #

# 模拟次数
number_sims <- 10000
# 待测试的样本量
sample_sizes <- c(10, 100, 1000, 5000)
# 创建列表存储不同样本量的模拟均值和方差
sim_results <- lapply(sample_sizes, function(n) {
  data.frame(
    mean = numeric(number_sims),
    var = numeric(number_sims)
  )
})
names(sim_results) <- paste0("n_", sample_sizes)

# 执行蒙特卡洛模拟
for (sim in 1:number_sims) {
  for (i in seq_along(sample_sizes)) {
    n <- sample_sizes[i]
    sample_data <- rgamma(n, shape = 1)  # 真实参数a=1，b=1
    sim_results[[i]]$mean[sim] <- mean(sample_data)
    sim_results[[i]]$var[sim] <- var(sample_data)
  }
}

第二步：编写矩估计函数计算参数 #

# 自定义Gamma分布矩估计函数
gamma_mme <- function(sample_data) {
  x_bar <- mean(sample_data)
  s_sq <- var(sample_data)
  a_hat <- x_bar^2 / s_sq
  b_hat <- s_sq / x_bar
  return(list(a_estimate = a_hat, b_estimate = b_hat))
}

# 对单一样本计算矩估计值（以样本量1000为例）
test_sample <- rgamma(1000, shape = 1)
estimates <- gamma_mme(test_sample)
cat("单样本矩估计结果：a =", round(estimates$a_estimate, 3), "，b =", round(estimates$b_estimate, 3), "\n")

# 查看样本量1000时所有模拟的估计值均值
n1000_estimates <- data.frame(
  a = sim_results$n_1000$mean^2 / sim_results$n_1000$var,
  b = sim_results$n_1000$var / sim_results$n_1000$mean
)
cat("n=1000时，a的估计均值：", round(mean(n1000_estimates$a), 3), "\n")
cat("n=1000时，b的估计均值：", round(mean(n1000_estimates$b), 3), "\n")

说明 #

• 完全无需额外包，自行编写简单函数即可实现矩估计，逻辑清晰易理解
• 样本量越大，矩估计结果越接近真实参数（本例中真实参数为a=1，b=1）

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