在R中用EM算法计算Gumbel与非结构Student-t混合Copula模型的混合系数

核心思路 #

混合Copula的EM算法分两大核心步骤：E步计算每个样本归属各Copula成分的后验概率（即混合权重的期望），M步通过最大化加权似然更新混合系数与各Copula的参数。由于mixtools主要适配常规分布混合场景，对Copula的支持有限，因此需要手动搭建EM框架，结合copula包完成Copula密度的计算与参数优化。

步骤分解与代码实现 #

1. 环境准备与模拟数据生成 #

先加载依赖包并生成测试用的混合Copula数据：

library(copula)
library(ggplot2)

# 模拟混合Copula数据
set.seed(123)
n <- 500
w_true <- 0.6  # 真实混合权重：Gumbel占60%，Student-t占40%

# 定义两个Copula成分
gumbel_cop <- gumbelCopula(theta = 1.8, dim = 2)
t_cop <- tCopula(param = c(0.7), dim = 2, df = 5)

# 生成混合数据
u1 <- rCopula(n * w_true, gumbel_cop)
u2 <- rCopula(n * (1 - w_true), t_cop)
u <- rbind(u1, u2)
set.seed(NULL)

2. EM算法迭代实现 #

参数初始化

# 初始化混合权重、Gumbel的theta、Student-t的rho和自由度
w <- 0.5
theta_g <- 1.5
rho_t <- 0.5
df_t <- 4

# 定义对数似然函数（用于监控收敛）
log_lik <- function(w, theta_g, rho_t, df_t) {
  gumbel_d <- dCopula(u, gumbelCopula(theta_g, dim=2))
  t_d <- dCopula(u, tCopula(rho_t, dim=2, df=df_t))
  sum(log(w * gumbel_d + (1 - w) * t_d))
}

# 迭代设置
max_iter <- 1000
tol <- 1e-6
log_lik_history <- numeric(max_iter)
log_lik_history[1] <- log_lik(w, theta_g, rho_t, df_t)

EM循环迭代

for (iter in 2:max_iter) {
  # E步：计算每个样本的后验概率（责任度）
  gumbel_d <- dCopula(u, gumbelCopula(theta_g, dim=2))
  t_d <- dCopula(u, tCopula(rho_t, dim=2, df=df_t))
  numerator_g <- w * gumbel_d
  numerator_t <- (1 - w) * t_d
  total <- numerator_g + numerator_t
  
  gamma_g <- numerator_g / total  # 样本属于Gumbel的后验概率
  gamma_t <- numerator_t / total  # 样本属于Student-t的后验概率
  
  # M步：更新参数
  # 更新混合权重
  w <- mean(gamma_g)
  
  # 更新Gumbel Copula的theta：最大化加权似然
  opt_g <- optim(par = theta_g, 
                 fn = function(x) -sum(gamma_g * log(dCopula(u, gumbelCopula(x, dim=2)))),
                 method = "L-BFGS-B", lower = 1.0001, upper = 10)
  theta_g <- opt_g$par
  
  # 更新Student-t Copula的rho和df：最大化加权似然
  opt_t <- optim(par = c(rho_t, df_t),
                 fn = function(par) -sum(gamma_t * log(dCopula(u, tCopula(par[1], dim=2, df=par[2])))),
                 method = "L-BFGS-B", lower = c(-0.9999, 2.0001), upper = c(0.9999, 30))
  rho_t <- opt_t$par[1]
  df_t <- opt_t$par[2]
  
  # 检查收敛
  log_lik_history[iter] <- log_lik(w, theta_g, rho_t, df_t)
  if (abs(log_lik_history[iter] - log_lik_history[iter-1]) < tol) {
    cat("收敛于迭代次数：", iter, "\n")
    break
  }
}

3. 结果验证与可视化 #

cat("估计的混合权重w：", round(w, 3), "\n")
cat("真实混合权重w：", w_true, "\n")
cat("估计的Gumbel theta：", round(theta_g, 3), "\n")
cat("估计的Student-t rho：", round(rho_t, 3), "\n")
cat("估计的Student-t df：", round(df_t, 3), "\n")

# 绘制对数似然收敛曲线
ggplot(data.frame(iter=1:iter, log_lik=log_lik_history[1:iter]), aes(x=iter, y=log_lik)) +
  geom_line(color="blue") +
  labs(x="迭代次数", y="对数似然") +
  theme_minimal()

关键注意事项 #

• 参数约束：Gumbel Copula的theta必须大于1，Student-t的rho需在(-1,1)区间内，自由度需大于2，优化时必须设置合理的上下界。
• 初始值敏感性：EM算法易陷入局部最优，建议多尝试几组不同初始值对比结果。
• 数值稳定性：当Copula密度极小时可能出现数值下溢，可直接在对数空间计算加权似然避免该问题。
• 高维扩展：若处理高维非结构Student-t Copula，参数为正定相关矩阵，可借助copula包中ellipCopula的相关工具完成约束优化。

参考资料 #

• copula包官方文档：涵盖各类Copula的密度计算、参数优化与模拟函数细节
• 混合模型EM算法标准推导：统计教材中混合模型章节的核心逻辑，即E步计算责任度、M步最大化加权似然

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