在R中拟合Gamma混合分布：适配性检验问题求助

我发现双Gamma混合分布似乎能较好拟合我的数据，但为确保适配性，需要通过Kolmogorov-Smirnov（KS）检验与卡方拟合优度检验验证，同时希望了解其他验证混合分布适配性的方法。

拟合实现与可视化 #

我使用R语言mixtools包的gammamixEM函数估计了两个Gamma分布的参数，从直方图来看拟合效果尚可：

foo_CE <- gammamixEM(data_cexio$diff)
x1 <- seq(0, 70, 1)
hist(data_cexio$diff, freq= FALSE, col="grey",border="black")
lines(density(data_cexio$diff), lwd = 2, col = "blue")
lines(x1, foo_CE$lambda[1]*dgamma(x1, foo_CE$gamma.pars[1,1], 1/foo_CE$gamma.pars[2,1]),
      col="orange", lwd=3)
lines(x1, foo_CE$lambda[2]*dgamma(x1, foo_CE$gamma.pars[1,2], 1/foo_CE$gamma.pars[2,2]),
      col="magenta", lwd=3)
legend("topright", c("0.20 * Gamma ~ (1.70,0.24)","0.80 * Gamma ~ (6.40,0.22)" ), fill=c( "orange","magenta"))

KS检验尝试及结果 #

我针对Gamma混合分布修改了KS检验的实现：

# 双Gamma混合分布的CDF
CDF_gamma <- function(x, shape, scale, p) {
  p[1]*pgamma(x,shape[1],1/scale[1]) + p[2]*pgamma(x,shape[2],1/scale[2])
}
test_CE <- ks.test(data_cexio$diff, CDF_gamma, shape=foo_CE$gamma.pars[1,], scale=foo_CE$gamma.pars[2,], p=foo_CE$lambda)

但检验得到的p值<2.2e-16，根据KS检验结果拟合效果不佳。

卡方拟合优度检验实现方法 #

卡方拟合优度检验需要将数据分组后对比观测频数与理论期望频数，具体R代码实现如下：

# 1. 对数据分组（确保每组期望频数≥5，可调整分组数量）
breaks <- quantile(data_cexio$diff, seq(0, 1, length.out = 10)) # 分成9个区间
groups <- cut(data_cexio$diff, breaks = breaks, include.lowest = TRUE)
obs_freq <- table(groups)

# 2. 计算每个区间的理论概率与期望频数
prob_intervals <- diff(CDF_gamma(breaks, shape=foo_CE$gamma.pars[1,], scale=foo_CE$gamma.pars[2,], p=foo_CE$lambda))
exp_freq <- length(data_cexio$diff) * prob_intervals

# 3. 执行卡方检验（若有组期望频数<5，需合并相邻组）
chi_sq_test <- chisq.test(obs_freq, p = prob_intervals, rescale.p = TRUE)

注意：分组时必须保证每组的期望频数不小于5，否则检验结果会失真；若存在期望频数不足的组，需合并相邻区间调整分组。

其他混合分布适配性验证方法 #

• Q-Q图检验：绘制样本分位数与混合Gamma分布的理论分位数对比图，若点基本贴合直线则说明拟合效果好

  # 生成理论分位数
  theoretical_quantiles <- qgamma(ppoints(length(data_cexio$diff)), 
                                  shape=foo_CE$gamma.pars[1,1], scale=foo_CE$gamma.pars[2,1])*foo_CE$lambda[1] +
                          qgamma(ppoints(length(data_cexio$diff)), 
                                  shape=foo_CE$gamma.pars[1,2], scale=foo_CE$gamma.pars[2,2])*foo_CE$lambda[2]
  qqplot(data_cexio$diff, theoretical_quantiles)
  abline(a=0, b=1, col="red")
  

• 似然比检验：对比混合Gamma模型与单一Gamma模型的对数似然值，判断混合模型是否显著更优
• 残差分析：计算每个样本点的观测密度与模型预测密度的残差，若残差随机分布在0附近则拟合效果理想
• 交叉验证：将数据集拆分为训练集和测试集，用训练集拟合模型后，评估测试集的拟合程度（如计算K-L散度）

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