寻找具有最低AIC值的fitdist类变量名称的快捷方法

我明白你遇到的痛点了——当数据里零值占比太高时，fitdistrplus包的gofstat()会直接罢工，没法依赖它来筛选最优分布。那咱们换个思路：手动遍历候选分布，尝试拟合后计算AIC值，挑出AIC最小的有效拟合结果，再提取对应的分布名称、均值、标准差等统计量。

下面是完整的实现函数，我已经考虑了拟合失败的情况（比如某些分布不能处理零值），用tryCatch跳过报错的拟合：

# 先确保加载必要的包
if (!requireNamespace("fitdistrplus", quietly = TRUE)) {
  stop("需要安装fitdistrplus包，请运行 install.packages('fitdistrplus')")
}
library(fitdistrplus)

fit_best_distribution <- function(data, candidate_dists = c("norm", "pois", "lnorm", "weibull", "nbinom")) {
  # 存储有效拟合结果的列表
  fit_results <- list()
  
  # 遍历每个候选分布，尝试拟合并计算AIC
  for (dist in candidate_dists) {
    # 用tryCatch捕获拟合失败的情况，避免中断循环
    fit <- tryCatch(
      fitdist(data, distr = dist),
      error = function(e) NULL
    )
    
    # 如果拟合成功，计算AIC并保存
    if (!is.null(fit)) {
      aic <- AIC(fit)
      fit_results[[dist]] <- list(fit_obj = fit, aic = aic)
    }
  }
  
  # 检查是否有有效拟合结果
  if (length(fit_results) == 0) {
    stop("所有候选分布均拟合失败，请检查数据或更换候选分布")
  }
  
  # 找出AIC最小的分布
  aic_values <- sapply(fit_results, function(x) x$aic)
  best_dist_name <- names(which.min(aic_values))
  best_fit <- fit_results[[best_dist_name]]$fit_obj
  
  # 根据分布类型计算均值和标准差
  calculate_stats <- function(fit_obj, dist_name) {
    params <- fit_obj$estimate
    switch(dist_name,
           "norm" = {
             mean_val <- params["mean"]
             sd_val <- params["sd"]
           },
           "pois" = {
             mean_val <- params["lambda"]
             sd_val <- sqrt(params["lambda"])
           },
           "lnorm" = {
             mean_val <- exp(params["meanlog"] + params["sdlog"]^2 / 2)
             sd_val <- mean_val * sqrt(exp(params["sdlog"]^2) - 1)
           },
           "weibull" = {
             shape <- params["shape"]
             scale <- params["scale"]
             mean_val <- shape * gamma(1 + 1/shape)
             var_val <- shape^2 * (gamma(1 + 2/shape) - (gamma(1 + 1/shape))^2)
             sd_val <- sqrt(var_val)
           },
           "nbinom" = {
             size <- params["size"]
             mu <- params["mu"]
             mean_val <- mu
             sd_val <- sqrt(mu + mu^2 / size)
           },
           # 如果是其他分布，可以在这里扩展
           {
             warning("暂不支持该分布的均值和标准差计算，返回NA")
             mean_val <- NA
             sd_val <- NA
           }
    )
    return(list(mean = mean_val, sd = sd_val))
  }
  
  stats <- calculate_stats(best_fit, best_dist_name)
  
  # 返回结果向量
  result <- c(
    distribution = best_dist_name,
    mean = unname(stats$mean),
    sd = unname(stats$sd),
    aic = best_fit$aic
  )
  
  return(result)
}

代码说明： #

• 候选分布选择：默认选了几个常用且对零值兼容性较好的分布（比如泊松、负二项分布天生支持零值；正态、对数正态、威布尔如果有零值可能需要调整，不过tryCatch会自动跳过拟合失败的情况），你可以根据自己的数据类型修改candidate_dists参数。
• 拟合容错处理：用tryCatch捕获拟合过程中的错误，比如gamma分布遇到零值会报错，这时候会直接跳过该分布，不影响其他分布的拟合。
• AIC筛选逻辑：遍历所有成功拟合的分布，找出AIC值最小的那个（AIC越小，拟合效果越好）。
• 统计量计算：针对不同分布的参数，手动计算均值和标准差（因为不同分布的参数和均值/标准差的对应关系不一样，比如威布尔分布需要用伽马函数转换）。

使用示例： #

# 生成带大量零值的测试数据
set.seed(123)
test_data <- c(rpois(50, 0.5), rep(0, 100))

# 运行函数
best_fit_result <- fit_best_distribution(test_data)
print(best_fit_result)

这样就能得到最优分布的名称、均值、标准差和对应的AIC值啦，完美避开gofstat()失效的问题。

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