my.densityplot <- function(x, cols="gray90", probs=c(1,3)/4, verticals=FALSE, ...){
    x <- na.omit(x)
    qtl <- quantile(x, probs=probs)
    panel.densityplot(x, ...)
    if(verticals){
        panel.abline(v=qtl)
    }
    dx <- density(x)
    breaks <- sort(c(range(dx$x), qtl))
    fx <- approxfun(dx$x, dx$y)
    do.polygon <- function(x){
        y <- fx(x)
        return(list(x=c(min(x), x, max(x)), y=c(0, y, 0)))
    }
    seqs <- lapply(1:(length(breaks)-1),
                   function(i){
                       x <- seq(breaks[i], breaks[i+1], l=20)
                       do.polygon(x)
                   })
    for(i in 1:length(seqs)){
        seqs[[i]]$col <- cols[i]
    }
    lapply(seqs, function(i) do.call(panel.polygon, i))
    panel.mathdensity(dmath=dnorm, col="black",
                      args=list(mean=mean(x),sd=sd(x)), lty=3)
}