sm.density.compare(): отображение нескольких оценок плотности на одном графике

Question

Я пытаюсь наложить три графика различной плотности в R, чтобы создать один график, отображающий все три линии (наложение). у меня есть sm пакет установлен / загружен, но я попытался использовать его с моими данными, но безрезультатно. Я создал три отдельных графика данных, просто используя 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")

Я получаю сообщение об ошибке:

Ошибка в sm.density.compare (dx, dy, model = "equal"):
sm.density.compare может обрабатывать только 1-ю трассировку данных:

Кто-нибудь знает, как я могу это исправить? Я много исследовал, но безуспешно. Я довольно плохо знаком с 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 представляет собой матрицу значений оценки плотности. (У Адриана здесь есть опечатка, это должна быть "оценка", а не "оценка", LOL).

Все функции от 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)

В этой ситуации оценка плотности для каждой группы выполняется независимо. Каждый объект "плотность" представляет собой список, например:

> 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 разные, из-за независимой оценки. Тем не менее, вы можете вручную указать универсальный bw при звонке density() используя аргумент 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)