Density plot with shadow colored

Category:R

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
# transabi.R
# statistics on abi expression data set
#
 
library("Hmisc")
library("lattice")
library("cluster")
library("gdata") #nobs
library("RColorBrewer")
 
polysection <- function(a,b,dist=dnorm,col="blue",n=11){
    dx <- seq(a,b,length.out=n)
    polygon(c(a,dx,b),c(0,dist(dx),0),border=NA,col=col)
}
 
polysection.funct <- function(a, b, x, y, col="blue",n=5) {
    sx <- seq(a, b, 5);
    polygon(c(x[a],x[sx],x[b]),c(0,y[sx],0),
                border=NA,col=col);
         
}
 
 
nice.dnorm <- function() {
 
x <- seq(-4,4,.1)
plot(x,dnorm(x),type="n",xaxs="i",yaxs="i",ylim=c(0,.4),bty="l",xaxt="n")
library(RColorBrewer)
cols<-rev(brewer.pal(4,"Blues"))
for(i in 0:3){
   polysection(i,i+1,col=cols[i+1])
   polysection(-i-1,-i,col=cols[i+1])
}
 
sx <- -3:3
segments(sx,0,sx,dnorm(sx),col="white")
lines(x,dnorm(x))
axis(1,sx,sx)
 
}
 
density.integrand <- function(x, dens){
    return(dens$y[x]);
}
 
integrate.density <- function(a, b, dens) {
    loc_a <- which(dens$x <= a);
    if (length(loc_a) == 0) loc_a <- 1;
    loc_b <- which(dens$x <= b);
    if (length(loc_b) == 0) loc_b <- 1;
    p <- loc_a[length(loc_a)];
    q <- loc_b[length(loc_b)];
 
    int<-integrate(density.integrand, lower=p, upper=q, dens=dens,
        subdivisions=1000);
     
}
 
nice.density <-function(df, ...) {
 
    m <- mean(df);
    d <- density(as.double(as.matrix(df)));
    plot(d, ...);
 
    w <- which(d$x <= m);
    loc <- w[length(w)];
 
    sx <- seq(loc, length(d$x), 50);
    sy <- seq(loc, 0, -50);
 
 
     
    segments(d$x[sx], 0, d$x[sx], d$y[sx], col="lightblue");
    segments(d$x[sy], 0, d$x[sy], d$y[sy], col="lightblue");
 
    library(RColorBrewer)
    cols<-rev(brewer.pal(length(sx),"Blues"))
    int <- integrate(density.integrand, lower=1, upper=length(d$x),
        dens=d, subdivisions=10000);
     
     
    for (i in 1:(length(sx)-1)) {
        sub<-integrate(density.integrand, lower=sx[i], upper=sx[i+1],
        dens=d, subdivisions=50000);
     
        polysection.funct(sx[i],sx[i+1], d$x, d$y, col=cols[i]);
        val<-as.character(round(sub$value/int$value, digits=3));
        text((d$x[sx[i]]+d$x[sx[i+1]])/2, 0.2, val);
    }
 
    for (i in 1:(length(sy)-1)) {
        sub<-integrate(density.integrand, lower=sy[i+1], upper=sy[i],
        dens=d, subdivisions=50000);
     
        polysection.funct(sy[i+1],sy[i], d$x, d$y, col=cols[i]);
        val<-as.character(round(sub$value/int$value, digits=3));
        text((d$x[sy[i]]+d$x[sy[i+1]])/2, 0.2, val);
    }
    text(d$x[loc], 0.5, paste("mean=", round(m,digits=3)));
}
 
### example
nice.dnorm()
nice.density(runif(1000))
Homepage
Comments

Hide Comments