2010-08-16 12 views
84

Uso frequentemente i grafici della densità del kernel per illustrare le distribuzioni. Questi sono facili e veloci per creare in R in questo modo:Ombreggiamento di un grafico di densità del kernel tra due punti.

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

che mi dà questo piccolo PDF:

enter image description here

mi piacerebbe per ombreggiare l'area sotto il PDF dal 75 ° al 95esimo percentile. E 'facile calcolare i punti utilizzando la funzione di quantile:

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

ma come faccio a ombreggiare l'area tra il q75 e q95?

+0

Potete fornire esempio di ombreggiatura l'esterno della vostra gamma contro l'interno della vostra gamma? Grazie. – Milktrader

risposta

67

Con la funzione polygon(), vedere la sua pagina di aiuto e credo che anche qui abbiamo avuto domande simili.

È necessario trovare l'indice dei valori quantili per ottenere le coppie effettive (x,y).

Edit: Qui si va:

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")) 

uscita (aggiunto dalla JDL)

enter image description here

+3

Non avrei mai potuto farlo funzionare se non avessi fornito la struttura. Grazie! –

+1

È una di quelle cose ... che sono state in 'demo (grafica)' da prima dell'alba del giorno in cui si incontra ogni tanto. Stessa idea per l'ombreggiatura della regressione NBER, ecc. –

+1

ohhhh. SAPEVO CHE L'avevo visto da qualche parte ma non riuscivo a strappare dal mio indice mentale dove l'avevo visto. Sono contento che il tuo indice mentale sia migliore del mio. –

63

Un'altra soluzione:

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) 

Risultato: alt text

+2

hey che è fantastico! e pieno di qualità ggplot! –

19

Una soluzione ampliata:

Se si voleva ombra entrambe le code (copia & incolla di codice di Dirk) e utilizzare valori x noti:

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")) 

Risultato:

2-tailed poly

+0

Ho il file png e l'ho ospitato su freeimagehosting, e potrebbe non caricarsi perché ... Non sono sicuro. – Milktrader

+0

File molto sfocato.Puoi ricreare e * caricare direttamente qui * SO ha il proprio servizio di server per questo? –

+0

Mi dispiace, ma non riesco a vedere come caricarlo direttamente in SO. – Milktrader

17

Questa domanda richiede una risposta lattice. Ecco uno molto semplice, semplicemente adattando il metodo impiegato da Dirk e altri:

#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)) 

enter image description here

Problemi correlati