sm.density.compare()：在单个图中显示多个密度估计

Question

我正在尝试在 R 中叠加三个不同的密度图，以创建一个显示所有三条线的图（叠加）。我有 sm 包 installed/loaded 但我尝试将它与我的数据一起使用但无济于事。我简单地使用 density() 创建了三个单独的数据图并绘制了值。我的代码如下所示：

library(sm)

set.seed(0)
x <- rnorm(100, 0, 1)
y <- rnorm(126, 0.3, 1.2)
z <- rnorm(93, -0.5, 0.7)
dx <- density(x)
dy <- density(y)
dz <- density(z)

plot(dx)
plot(dy)
plot(dz)

但是当我尝试使用 sm.density.compare() 叠加图表时：

sm.density.compare(dx,dy,model="equal")

我收到一条错误消息：

Error in sm.density.compare(dx, dy, model = "equal") :
sm.density.compare can handle only 1-d data trace:

有人知道我该如何解决这个问题吗？我研究了很多但没有成功。我是 R 的新手，真的可以使用帮助。

Answer 1

如果要使用sm.density.compare()，则不要使用density()。

sm.density.compare()本身就是在做密度估计。具体来说就是对分组数据做密度估计，这样就可以在同一张图上绘制不同组的密度。

这是您真正需要做的事情：

## three groups, each of length: length(x), length(y), length(z)
group.index <- rep(1:3, c(length(x), length(y), length(z)))
## collect data together and use sm.density.compare()
den <- sm.density.compare(c(x,y,z), group = group.index, model = "equal")
## plot will be generated automatically

当使用model = "equal"时，sm.density.compare()有返回值。看看 str(den):

List of 4
 $ p          : num 0
 $ estimaate  : num [1:3, 1:50] 2.37e-07 3.81e-06 6.06e-10 2.17e-06 2.26e-05 ...
 $ eval.points: num [1:50] -4.12 -3.94 -3.76 -3.58 -3.4 ...
 $ h          : num 0.376

h包含用于所有密度估计的带宽，eval.points包含估计点，而estimaate是密度估计值的矩阵。（Adrian这里打错了，应该是"estimate"，不是"estimaate"，哈哈）。

sm 包中的所有函数，以前缀 sm. 开头，接受可选参数 ...，传递给 sm.options。阅读 ?sm.options，您会发现您可以完全控制颜色显示、线型和线宽、带宽选择方法等。

参考波段只会添加到两组数据中。即，对于成对比较，sm.density.compare() 可以做得更多。例如：

den2 <- sm.density.compare(c(x,y), group = rep(1:2, c(length(x), length(y))),
                           model = "equal")

> str(den2)
List of 6
 $ p          : num 0.22
 $ estimate   : num [1:2, 1:50] 4.92e-06 2.70e-05 2.51e-05 1.00e-04 1.09e-04 ...
 $ eval.points: num [1:50] -4.12 -3.94 -3.76 -3.58 -3.4 ...
 $ upper      : num [1:50] 0.00328 0.00373 0.00459 0.00614 0.00886 ...
 $ lower      : num [1:50] 0 0 0 0 0 ...
 $ h          : num 0.44

其中 lower 和 upper 给出了参考带/置信区域的界限。

如果使用density()，则不要使用sm.density.compare()

## set universal estimation range
xlim <- range(x, y, z)
dx <- density(x, from = xlim[1], to = xlim[2], n = 200)
dy <- density(y, from = xlim[1], to = xlim[2], n = 200)
dz <- density(z, from = xlim[1], to = xlim[2], n = 200)

在这种情况下，每个组的密度估计都是独立完成的。每个"density"对象都是一个列表，例如：

> str(dx)
List of 7
 $ x        : num [1:200] -2.64 -2.61 -2.58 -2.55 -2.52 ...
 $ y        : num [1:200] 0.023 0.026 0.0291 0.0323 0.0356 ...
 $ bw       : num 0.31
 $ n        : int 100
 $ call     : language density.default(x = x, n = 200, from = xlim[1], to = xlim[2])
 $ data.name: chr "x"
 $ has.na   : logi FALSE
 - attr(*, "class")= chr "density"

x是评估点，y是估计密度，bw是使用的带宽。您会看到，由于独立估计，dx$bw、dy$bw 和 dz$bw 是不同的。但是，您可以在调用 density() 时通过使用参数 bw 手动指定通用 bw。参见?density，这里就不举例了。

现在，要叠加这些密度图，您需要自己动手。

## set global plotting range
ylim <- range(dx$y, dy$y, dz$y)
## make plot
plot(dx$x, dx$y, col = 1, lwd = 2, type = "l", xlim = xlim, ylim = ylim)
lines(dy$x, dy$y, col = 2, lwd = 2)
lines(dz$x, dz$y, col = 3, lwd = 2)

sm.density.compare()：在单个图中显示多个密度估计

sm.density.compare(): displaying multiple density estimations in a single plot

r

density-plot