如何在R的ggplot2中修复Beta分布叠加及更改直方图Y轴为频率？

我编写了如下代码，希望绘制直方图并叠加多种密度函数，除Beta分布外其余均正常。问题源于我为计算Beta密度函数将数据归一化到(0,1)，但直方图仍显示原始数据。

library(ggplot2)
library(rlang)
library(MASS)

normalize_btw_zero_and_one <- function(x){(x-min(x)+.00000001)/(max(x)-min(x)+.00000002)}

get_histogram <- function(data_set, column_name, bin_width, attribute_name) {
  colname     <- as_label(enquo(column_name))
  gamma_par   <- MASS :: fitdistr(data_set[[colname]], "gamma")
  weibull_par <- MASS :: fitdistr(data_set[[colname]], "weibull")
  beta_par    <- MASS :: fitdistr(normalize_btw_zero_and_one(data_set[[colname]]), dbeta, start = list(shape1 = 1, shape2 = 1))
  ggplot(data_set, aes(x= {{ column_name }})) +
    geom_histogram(aes(y= after_stat(density)), binwidth = bin_width, fill= "lightblue", colour="black") +
    xlab(paste0(attribute_name)) +
    stat_function(fun = dnorm , args= list(mean= mean(data_set[[colname]]), sd= sd(data_set[[colname]])), 
                  mapping = aes(colour = "Normal"))+
    stat_function(fun = dlnorm, args = list(meanlog= mean(log(data_set[[colname]])), sdlog= sd(log(data_set[[colname]]))),
                  mapping = aes(colour = "LogNormal")) + 
    stat_function(fun = dgamma, args= list(shape= gamma_par$estimate[[1]]  , rate=gamma_par$estimate[[2]] ),
                  mapping = aes(colour= "Gamma"))+
    stat_function(fun = dweibull, args= list(shape= weibull_par$estimate[[1]]  , scale=weibull_par$estimate[[2]] ),
                  mapping = aes(colour= "Weibull"))+
    stat_function(fun = dexp, mapping = aes(colour = "Exponential")) +
    stat_function(fun = dbeta, args = list(shape1= beta_par$estimate[[1]], shape2= beta_par$estimate[[2]] ),
                  mapping = aes(color = "Beta")) +
    scale_colour_manual("Distribution", values = c("red", "blue", "lightgreen","pink", "purple" , "yellow"))
}

set.seed(30333)
test_dt <- rnorm(10000, 30, 1)
df <- data.frame(test_dt)
rm(test_dt)

get_histogram(df, test_dt, .3, "Test ")

问题1：修复Beta分布函数的叠加 #

Beta分布的支持区间是(0,1)，直接用dbeta计算的是归一化后数据的密度，和原始数据区间不匹配，需要做变量变换的密度校正：

1. 在get_histogram函数中先保存原始数据的极值，用于后续转换：

x_min <- min(data_set[[colname]])
x_max <- max(data_set[[colname]])
beta_par <- MASS::fitdistr(normalize_btw_zero_and_one(data_set[[colname]]), dbeta, start = list(shape1 = 1, shape2 = 1))

2. 自定义适配原始数据区间的Beta密度函数（加入雅可比行列式校正，保证密度值和原始数据匹配）：

dbeta_original <- function(x, shape1, shape2, x_min, x_max) {
  z <- (x - x_min + 1e-8) / (x_max - x_min + 2e-8)
  # 变量变换后密度需乘以dx/dz的倒数
  dbeta(z, shape1, shape2) / (x_max - x_min + 2e-8)
}

3. 替换原Beta对应的stat_function：

stat_function(fun = dbeta_original, 
              args = list(shape1= beta_par$estimate[[1]], 
                          shape2= beta_par$estimate[[2]],
                          x_min = x_min,
                          x_max = x_max),
              mapping = aes(color = "Beta"))

问题2：将Y轴改为频率而非密度 #

频率 = 密度 × 样本量 × 组距，需要同步修改直方图和所有密度曲线的计算逻辑：

1. 修改直方图的y轴映射为计数（频率）：

geom_histogram(aes(y= after_stat(count)), binwidth = bin_width, fill= "lightblue", colour="black")

2. 在函数内计算样本量：

n <- nrow(data_set)

3. 给每个stat_function的密度值乘以n * bin_width转换为频率，以正态分布为例：

stat_function(fun = function(x) dnorm(x, mean= mean(data_set[[colname]]), sd= sd(data_set[[colname]])) * n * bin_width, 
              mapping = aes(colour = "Normal"))

Beta分布的频率函数也要同步调整：

dbeta_original_freq <- function(x, shape1, shape2, x_min, x_max, n, bin_width) {
  z <- (x - x_min + 1e-8) / (x_max - x_min + 2e-8)
  dbeta(z, shape1, shape2) / (x_max - x_min + 2e-8) * n * bin_width
}

stat_function(fun = dbeta_original_freq, 
              args = list(shape1= beta_par$estimate[[1]], 
                          shape2= beta_par$estimate[[2]],
                          x_min = x_min,
                          x_max = x_max,
                          n = n,
                          bin_width = bin_width),
              mapping = aes(color = "Beta"))

其他分布的stat_function按照同样逻辑修改即可。

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