如何在R语言中为样本拟合Erlang混合模型？

Erlang分布是Gamma分布的特例，其shape参数必须为正整数，这是拟合Erlang混合模型时需要重点约束的核心条件。以下是针对你生成的样本数据，实现Erlang混合模型拟合的具体步骤：

1. 样本数据准备（你的原始代码） #

首先运行你提供的代码生成样本并可视化：

set.seed(123)
n <- 500
data1 <- rgamma(n / 2, shape = 9, rate = 10) 
data2 <- mean(data1) + rgamma(n / 2, shape = 5, rate = 5)
sample <- c(data1, data2)

# 绘制直方图与核密度曲线
hist(sample, breaks = 50, prob = T)
lines(density(sample), lwd = 2, col = "red")

2. 手动构建似然函数拟合（最灵活的方式） #

由于多数混合模型包默认支持Gamma分布（无整数约束），我们可以手动定义对数似然函数，通过优化器约束shape为整数来拟合：

步骤2.1 定义对数似然函数 #

# 混合Erlang分布的对数似然函数
log_likelihood <- function(params, x) {
  pi1 <- params[1]
  k1 <- round(params[2])  # 强制shape为整数
  lambda1 <- params[3]
  k2 <- round(params[4])
  lambda2 <- params[5]
  
  # 参数合法性校验，避免无效值
  if (pi1 <= 0 || pi1 >= 1 || k1 < 1 || k2 < 1 || lambda1 <= 0 || lambda2 <= 0) {
    return(-Inf)
  }
  
  # 计算混合密度并求和对数
  dens <- pi1 * dgamma(x, shape = k1, rate = lambda1) + (1 - pi1) * dgamma(x, shape = k2, rate = lambda2)
  sum(log(dens))
}

步骤2.2 设置初始参数 #

基于你之前的分组策略，估计初始参数（保证优化收敛）：

# 分组获取初始值
sample1 <- sample[sample < 1.3]
sample2 <- sample[sample > 1.3]

pi_init <- length(sample1)/length(sample)
# Gamma分布中shape ≈ mean²/var，取整作为Erlang的初始k
k1_init <- round(mean(sample1)^2 / var(sample1))
lambda1_init <- k1_init / mean(sample1)  # Gamma的mean=shape/rate → rate=shape/mean

k2_init <- round(mean(sample2)^2 / var(sample2))
lambda2_init <- k2_init / mean(sample2)

init_params <- c(pi_init, k1_init, lambda1_init, k2_init, lambda2_init)

步骤2.3 最大化对数似然 #

使用optim函数进行优化，约束参数范围：

fit <- optim(init_params, log_likelihood, x = sample, method = "L-BFGS-B", 
             lower = c(0.01, 1, 0.01, 1, 0.01), upper = c(0.99, 20, 20, 20, 20),
             control = list(fnscale = -1))  # fnscale=-1表示最大化似然

# 提取拟合结果
pi1_est <- fit$par[1]
k1_est <- round(fit$par[2])
lambda1_est <- fit$par[3]
k2_est <- round(fit$par[4])
lambda2_est <- fit$par[5]

# 打印参数
cat("拟合得到的混合Erlang模型参数：\n")
cat("组1：权重=", round(pi1_est, 3), "，shape(k)=", k1_est, "，rate(λ)=", round(lambda1_est, 3), "\n")
cat("组2：权重=", round(1-pi1_est, 3), "，shape(k)=", k2_est, "，rate(λ)=", round(lambda2_est, 3), "\n")

步骤2.4 可视化拟合结果 #

将拟合的混合密度曲线与原始数据对比：

x_seq <- seq(min(sample), max(sample), length.out = 1000)
mix_dens <- pi1_est * dgamma(x_seq, shape = k1_est, rate = lambda1_est) + 
  (1 - pi1_est) * dgamma(x_seq, shape = k2_est, rate = lambda2_est)

hist(sample, breaks = 50, prob = T, main = "混合Erlang模型拟合结果")
lines(density(sample), lwd = 2, col = "red", lty = 2, label = "核密度")
lines(x_seq, mix_dens, lwd = 2, col = "blue", label = "拟合Erlang混合")
legend("topright", legend = c("核密度", "拟合Erlang混合"), col = c("red", "blue"), lwd = 2, lty = c(2,1))

3. 使用flexmix包拟合（更简洁的封装方式） #

flexmix支持自定义分布组件，我们可以封装Erlang分布进行拟合：

install.packages("flexmix")
library(flexmix)

# 定义Erlang密度函数
dErlang <- function(x, k, rate) {
  dgamma(x, shape = k, rate = rate)
}

# 定义Erlang混合组件
erlang_dist <- FLXMRdist(dist = "Erlang", d = dErlang, 
                         start = function(x) {
                           k <- round(mean(x)^2 / var(x))
                           lambda <- k / mean(x)
                           c(k, lambda)
                         })

# 拟合2组分的混合模型
fit_flex <- flexmix(sample ~ 1, k = 2, model = erlang_dist)

# 查看结果
cat("组件权重：", prior(fit_flex), "\n")
cat("组件参数：\n")
print(parameters(fit_flex))

关键注意事项 #

• Erlang分布的shape参数必须为正整数，拟合时必须通过取整或整数约束实现，这是与普通Gamma混合模型的核心差异。
• 初始参数的选择会影响优化收敛结果，建议基于数据分组先获取合理初始值，避免陷入局部最优。
• 可以尝试多组初始参数，对比对数似然值选择最优拟合结果。

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