小样本下R中Beta分布叠加直方图失效的原因及代码疑问

问题背景 #

我有多个30至200行的小数据集，需要确定最优拟合分布，计划用可视化方法对比Normal和Beta分布。Beta分布要求数据落在[0,1]区间，所以需要对原始数据做缩放转换。JMP软件对任意样本量都能正常输出拟合图形，但在R中复现时，样本量小于80的小数据集会出现无意义结果（比如y轴密度值高达400），只有当数据量达到5000条时，R和JMP的结果才一致。

核心函数代码 #

get_histogram <- function(data_set, column_name, bin_width) {
  
  colname <- as_label(enquo(column_name))
  x <- data_set[[colname]]
  shift <- min(x)
  scale <- diff(range(x))
  norm_it <- function(x) (x - shift + 1e-8) / (diff(range(x)) + 2e-8)
  
  beta_par    <- fitdistr(norm_it(x), dbeta, 
                          start = list(shape1 = 1, shape2 = 1))$estimate
  
  ggplot(data_set, aes(x = {{ column_name }})) +
    geom_histogram(aes(y = after_stat(density)), binwidth = bin_width, 
                   fill = "lightblue", colour="black") +
    stat_function(fun = dnorm, args = list(mean = mean(x), sd = sd(x)), 
                  mapping = aes(colour = "Normal")) +
    stat_function(fun = function(x) dbeta(norm_it(x), beta_par[1], beta_par[2])/(scale),
                  mapping = aes(color = "Beta")) +
    scale_colour_manual("Distribution", 
                        values = c("red", "blue"))
}

测试代码 #

set.seed(2122)
test_dt <- rnorm(50, 30 , 2)
new_col <- sample(size = 50, x= c("a", "b", "c"), replace = TRUE)
df <- data.frame(test_dt, new_col)
rm(test_dt)
rm(new_col)

get_histogram(df, test_dt, .05)

疑问 #

1. 在stat_function()中叠加Beta分布时，分母应该用/ scale还是/ (scale + shift)？
2. 为什么该代码在小数据集（即使样本量大于30）下会失效？需要R相关的技术解决方案。

解答 #

问题1：Beta分布密度转换的分母选择 #

应该用/ scale，理由如下：
当你把原始数据x通过norm_it(x)转换到[0,1]区间时，核心是线性变换：
$$z = \frac{x - \text{shift}}{\text{scale}}$$
（代码中添加的1e-8是为了避免边界值等于0或1，不影响变换的核心逻辑）

根据概率密度的变量替换公式，若$z = g(x)$是单调可导变换，原始变量$x$的密度$f_X(x)$与变换后变量$z$的密度$f_Z(z)$满足：
$$f_X(x) = f_Z(g(x)) \times \left|\frac{dz}{dx}\right|^{-1}$$

这里$dz/dx = 1/\text{scale}$，逆变换的导数就是$\text{scale}$，因此$f_X(x) = f_Z(z) / \text{scale}$。shift是平移量，不影响密度的缩放比例，无需加入分母。

问题2：小数据集下失效的原因及解决方案 #

核心原因

1. Beta参数拟合不稳定：fitdistr对小样本的Beta分布参数估计容易发散，尤其是初始值shape1=1, shape2=1（对应均匀分布）与真实数据偏差较大时，极易得到极端形状参数（如趋近于0或无穷大），导致密度值异常飙升。
2. 直方图bin宽度不合理：测试代码中bin_width=.05与原始数据（均值30、标准差2，范围约25-35）不匹配，过小的bin宽度会生成大量空bin或极窄bin，而密度计算为频数/(样本量*bin宽度)，bin宽度太小会直接放大密度值。

解决方案

1. 优化Beta参数拟合稳定性
   • 改用fitdistrplus包的fitdist函数，支持更多优化算法与收敛控制：

     # 先安装包
     # install.packages("fitdistrplus")
     library(fitdistrplus)
     
     # 替换原代码中的beta_par计算部分
     beta_par <- fitdist(norm_it(x), "beta", start = list(shape1 = 1, shape2 = 1), 
                         lower = c(1e-8, 1e-8))$estimate
     

   • 添加参数下界约束，避免形状参数趋近于0。
2. 动态调整直方图bin宽度
   放弃手动设置过小的bin宽度，改用自适应方法：

   # 替换geom_histogram部分
   geom_histogram(aes(y = after_stat(density)), bins = nclass.Sturges(x), 
                  fill = "lightblue", colour="black")
   

   nclass.Sturges()会根据样本量自动计算合理的bin数量，避免密度值异常。
3. 增强边界处理鲁棒性
   优化norm_it函数，确保转换后数据严格落在(0,1)区间：

   norm_it <- function(x) {
     z <- (x - shift) / scale
     z[z <= 0] <- 1e-8
     z[z >= 1] <- 1 - 1e-8
     z
   }
   

4. 添加拟合结果校验
   在拟合后检查参数合理性，避免绘制无意义曲线：

   if(any(beta_par > 100 | beta_par < 0.01)) {
     warning("Beta分布参数估计异常，可能因样本量过小导致")
     beta_par <- c(1,1) # 回退到均匀分布参数
   }
   

内容的提问来源于stack exchange，提问作者Joe the Second