R密度图边缘效应与Y轴计数显示问题求助

一、基础R解决方案（优先方案） #

核心修改点 #

• 消除边缘效应：在density()函数中加入from和to参数，限制密度估计范围与绘图x轴完全一致，避免核函数在数据范围外的区域出现不必要衰减
• 密度转计数：将密度值乘以总样本权重（此处为28），由于密度曲线下面积为1，转换后曲线下面积等于总计数，Y轴数值直接对应计数密度
• 优化循环逻辑：简化data.table的模拟赋值语法，提升运行效率

修改后的完整代码 #

bw = 500
total_weight = sum(burialsToRun$Weight)
x_range = c(-9500, -3500)

# 初始密度计算并转换为计数
y <- density(check$sim, bw = bw, n=512, kernel = "gaussian", 
             weights = burialsToRun$Weight, from = x_range[1], to = x_range[2])
y$y <- y$y * total_weight # 密度转计数

# 自动计算Y轴最大值
ymax = max(y$y) * 1.1

par(mar=c(4.1,2.1,2.1,6.1))
plot(y, xlim = x_range, xaxs = "i", ylim = c(0, ymax), main="", col = rgb(0,0,0,0.05), axes = FALSE, ylab = " ", xlab = " ")
axis(1, at=seq(x_range[1], x_range[2], by = 500), las = 1)
axis(4, at=seq(0, ymax, by = ceiling(ymax/10)), las = 1)
mtext(text="Count", side=4, line=4, las=3) 
mtext(text="Calibrated date (cal. BC)", side=1, line=2.5, las=1)
box()

# 3000次循环绘图
for(i in 1:3000) {
  burialsToRun[, sim := runif(.N)*(End - Start) + Start]
  y_temp <- density(burialsToRun$sim, bw = bw, n=512, kernel = "gaussian", 
                    weights = burialsToRun$Weight, from = x_range[1], to = x_range[2])
  y_temp$y <- y_temp$y * total_weight
  lines(y_temp, col= rgb(0,0,0,0.05))
}

二、ggplot2解决方案及现有代码错误分析 #

现有代码的错误 #

1. 覆盖全局映射/数据：stat_density中设置mapping = NULL和data = NULL，导致无法读取check数据集的内容
2. 计数转换逻辑错误：直接使用after_stat(count)会得到「密度×样本量×带宽」的结果，和预期的计数密度逻辑不符
3. 未实现循环逻辑：ggplot基于数据框绘图，需要先生成所有模拟数据，再通过分组实现多次循环的线图

正确的ggplot2代码 #

library(ggplot2)
library(data.table)

bw = 500
total_weight = sum(burialsToRun$Weight)
x_range = c(-9500, -3500)

# 生成3000次模拟的长格式数据集
sim_data <- rbindlist(lapply(1:3000, function(i) {
  burialsToRun[, sim := runif(.N)*(End - Start) + Start]
  burialsToRun[, .(sim, Weight, iteration = i)]
}))

# 绘制所有模拟的密度线+初始数据的参考线
ggplot(sim_data, aes(x = sim)) +
  stat_density(aes(y = after_stat(density * total_weight), group = iteration),
               bw = bw, kernel = "gaussian", 
               from = x_range[1], to = x_range[2],
               color = rgb(0,0,0,0.05), geom = "line") +
  geom_density(data = check, aes(y = after_stat(density * total_weight)),
               bw = bw, kernel = "gaussian",
               from = x_range[1], to = x_range[2],
               color = "black", linewidth = 1) +
  scale_x_continuous(limits = x_range, breaks = seq(x_range[1], x_range[2], by = 500),
                     name = "Calibrated date (cal. BC)") +
  scale_y_continuous(position = "right", name = "Count") +
  theme_classic() +
  theme(plot.margin = unit(c(1, 1, 1, 1), "lines"))

代码说明 #

• 用rbindlist()生成包含3000次模拟的长数据，通过iteration列区分不同循环
• group = iteration让ggplot为每次模拟单独绘制密度线
• y = after_stat(density * total_weight)实现和基础R一致的密度转计数逻辑
• from和to参数消除边缘效应，可选叠加初始数据的黑色密度线作为参考

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