如何用ggplot2绘制难处理的重尾概率分布

背景 #

部分概率分布（如商分布，柯西分布是其典型特例）因重尾特性，难以生成平滑的可视化图像。以下是使用ggplot2实现的最小可复现示例：

require(ggplot2)
a <- rnorm(1e4, 2, 0.5)
b <- rnorm(1e4, 0.2, 0.1)
c <- a/b
ggplot() + geom_density(aes(c)) + coord_cartesian(xlim = c(-300, 500))

解决方案 #

针对重尾分布的特性，可通过以下几种方式实现平滑可视化：

1. 过滤极端值+稳健带宽选择

重尾分布的极端异常值会严重干扰核密度估计，先通过分位数过滤掉极端值，再使用适配重尾分布的带宽选择方法：

require(ggplot2)
a <- rnorm(1e4, 2, 0.5)
b <- rnorm(1e4, 0.2, 0.1)
c <- a/b

# 保留98%的核心数据（过滤首尾1%的极端值）
q_low <- quantile(c, 0.01)
q_high <- quantile(c, 0.99)
c_filtered <- c[c >= q_low & c <= q_high]

ggplot() + 
  geom_density(aes(c_filtered), bw = "SJ", kernel = "epanechnikov") + 
  labs(x = "a/b")

• bw = "SJ"采用Sheather-Jones带宽选择法，比手动调整更适配重尾分布
• epanechnikov核函数比默认高斯核更稳健，减少极端值的影响

2. 对数变换（适用于正取值为主的场景）

若商分布的取值以正数为主，可通过对数变换压缩尾部区间，再进行密度可视化：

# 仅保留正值并做对数变换
c_pos <- c[c > 0]
c_log <- log(c_pos)

ggplot() + 
  geom_density(aes(c_log), bw = "SJ") + 
  labs(x = "log(a/b)")

后续可通过逆变换还原原分布的形态特征。

3. 直方图叠加密度曲线

用直方图的分组特性弱化极端值干扰，同时叠加平滑的密度曲线，聚焦分布主体：

ggplot() + 
  geom_histogram(aes(c, y = after_stat(density)), bins = 50, alpha = 0.3) + 
  geom_density(aes(c), bw = "SJ", color = "red", linewidth = 1) + 
  coord_cartesian(xlim = c(quantile(c, 0.005), quantile(c, 0.995)))

4. ggridges::geom_density_ridges的适配方案

针对ggridges的情况，同样先过滤极端值，再调整密度参数：

require(ggridges)
# 构造分组示例数据
df <- data.frame(
  value = c,
  group = rep("Group1", length(c))
)

# 过滤极端值
df_filtered <- df[df$value >= q_low & df$value <= q_high, ]

ggplot(df_filtered, aes(x = value, y = group)) +
  geom_density_ridges(scale = 1.5, bw = "SJ", kernel = "epanechnikov")

内容的提问来源于stack exchange，提问作者Luka Seamus Wright