Mengarsir plot kepadatan kernel antara dua titik.

95

Saya sering menggunakan plot kepadatan kernel untuk menggambarkan distribusi. Ini mudah dan cepat untuk dibuat di R seperti:

set.seed(1)
draws <- rnorm(100)^2
dens <- density(draws)
plot(dens)
#or in one line like this: plot(density(rnorm(100)^2))

Yang memberi saya PDF kecil yang bagus ini:

masukkan deskripsi gambar di sini

Saya ingin menaungi area di bawah PDF dari persentil ke-75 hingga ke-95. Mudah untuk menghitung poin menggunakan quantilefungsi:

q75 <- quantile(draws, .75)
q95 <- quantile(draws, .95)

Tapi bagaimana cara menaungi area antara q75dan q95?

JD Long
sumber
Dapatkah Anda memberikan contoh bayangan bagian luar jangkauan Anda versus bagian dalam jangkauan Anda? Terima kasih.
Milktrader

Jawaban:

75

Dengan polygon()fungsinya, lihat halaman bantuannya dan saya yakin kami juga memiliki pertanyaan serupa di sini.

Anda perlu mencari indeks dari nilai-nilai kuantil untuk mendapatkan (x,y)pasangan yang sebenarnya .

Edit: Ini dia:

x1 <- min(which(dens$x >= q75))  
x2 <- max(which(dens$x <  q95))
with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="gray"))

Output (ditambahkan oleh JDL)

masukkan deskripsi gambar di sini

Dirk Eddelbuettel
sumber
3
Saya tidak akan pernah berhasil jika Anda tidak memberikan strukturnya. Terima kasih!
JD Long
2
Itu salah satu dari hal-hal itu ... yang telah ada demo(graphics)sejak sebelum fajar tepat waktu sehingga orang datang sesekali. Ide yang sama untuk bayangan regresi NBER dll.
Dirk Eddelbuettel
1
ohhhh. SAYA TAHU saya telah melihatnya di suatu tempat tetapi tidak dapat menarik dari indeks mental saya di mana saya telah melihatnya. Saya senang indeks mental Anda lebih baik dari saya.
JD Long
70

Solusi lain:

dd <- with(dens,data.frame(x,y))

library(ggplot2)

qplot(x,y,data=dd,geom="line")+
  geom_ribbon(data=subset(dd,x>q75 & x<q95),aes(ymax=y),ymin=0,
              fill="red",colour=NA,alpha=0.5)

Hasil:

teks alt

Ben Bolker
sumber
21

Solusi yang diperluas:

Jika Anda ingin membuat bayangan pada kedua ekor (salin & tempel kode Dirk) dan gunakan nilai x yang diketahui:

set.seed(1)
draws <- rnorm(100)^2
dens <- density(draws)
plot(dens)

q2     <- 2
q65    <- 6.5
qn08   <- -0.8
qn02   <- -0.2

x1 <- min(which(dens$x >= q2))  
x2 <- max(which(dens$x <  q65))
x3 <- min(which(dens$x >= qn08))  
x4 <- max(which(dens$x <  qn02))

with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="gray"))
with(dens, polygon(x=c(x[c(x3,x3:x4,x4)]), y= c(0, y[x3:x4], 0), col="gray"))

Hasil:

Poli berekor 2

Milktrader
sumber
Saya memiliki file png dan menyimpannya di freeimagehosting, dan file tersebut mungkin tidak dimuat karena ... Saya tidak yakin.
Milktrader
File sangat buram. Bisakah Anda membuatnya kembali dan mengunggahnya di sini secara langsung. SO memiliki layanan server sendiri untuk ini?
Dirk Eddelbuettel
Maaf, tapi saya tidak bisa melihat cara menguploadnya ke SO secara langsung.
Milktrader
18

Pertanyaan ini butuh latticejawaban. Ini yang paling mendasar, cukup mengadaptasi metode yang digunakan oleh Dirk dan lainnya:

#Set up the data
set.seed(1)
draws <- rnorm(100)^2
dens <- density(draws)

#Put in a simple data frame   
d <- data.frame(x = dens$x, y = dens$y)

#Define a custom panel function;
# Options like color don't need to be hard coded    
shadePanel <- function(x,y,shadeLims){
    panel.lines(x,y)
    m1 <- min(which(x >= shadeLims[1]))
    m2 <- max(which(x <= shadeLims[2]))
    tmp <- data.frame(x1 = x[c(m1,m1:m2,m2)], y1 = c(0,y[m1:m2],0))
    panel.polygon(tmp$x1,tmp$y1,col = "blue")
}

#Plot
xyplot(y~x,data = d, panel = shadePanel, shadeLims = c(1,3))

masukkan deskripsi gambar di sini

joran
sumber
3

Berikut ggplot2varian lain berdasarkan fungsi yang mendekati kepadatan kernel pada nilai data asli:

approxdens <- function(x) {
    dens <- density(x)
    f <- with(dens, approxfun(x, y))
    f(x)
}

Menggunakan data asli (daripada menghasilkan kerangka data baru dengan nilai x dan y perkiraan kepadatan) memiliki manfaat juga bekerja di plot segi di mana nilai-nilai kuantitatif bergantung pada variabel yang mengelompokkan data:

Kode digunakan

library(tidyverse)
library(RColorBrewer)

# dummy data
set.seed(1)
n <- 1e2
dt <- tibble(value = rnorm(n)^2)

# function that approximates the density at the provided values
approxdens <- function(x) {
    dens <- density(x)
    f <- with(dens, approxfun(x, y))
    f(x)
}

probs <- c(0.75, 0.95)

dt <- dt %>%
    mutate(dy = approxdens(value),                         # calculate density
           p = percent_rank(value),                        # percentile rank 
           pcat = as.factor(cut(p, breaks = probs,         # percentile category based on probs
                                include.lowest = TRUE)))

ggplot(dt, aes(value, dy)) +
    geom_ribbon(aes(ymin = 0, ymax = dy, fill = pcat)) +
    geom_line() +
    scale_fill_brewer(guide = "none") +
    theme_bw()



# dummy data with 2 groups
dt2 <- tibble(category = c(rep("A", n), rep("B", n)),
              value = c(rnorm(n)^2, rnorm(n, mean = 2)))

dt2 <- dt2 %>%
    group_by(category) %>% 
    mutate(dy = approxdens(value),    
           p = percent_rank(value),
           pcat = as.factor(cut(p, breaks = probs,
                                include.lowest = TRUE)))

# faceted plot
ggplot(dt2, aes(value, dy)) +
    geom_ribbon(aes(ymin = 0, ymax = dy, fill = pcat)) +
    geom_line() +
    facet_wrap(~ category, nrow = 2, scales = "fixed") +
    scale_fill_brewer(guide = "none") +
    theme_bw()

Dibuat pada 13-07-2018 oleh paket reprex (v0.2.0).

Matt Flor
sumber