Solved – How to draw mean, median and mode lines in R that end at density

data visualizationmeanmedianmoder

I have drawn a skewed distribution, using these commands:

x <- seq(-2.5, 10, length=1000000)
hx5 <- rnorm(x,0,1) + rexp(x,1/5) # tau=5 (rate = 1/tau)
plot(density(hx5), xlim=c(-2.5,10), type="l", col="green",
     xlab="x", main="ExGaussian  curve",lwd=2)

Now I want to draw three lines, for the mean, the mode and the median of the distribution. If I simply write, for example:

abline(v=median(hx5))

the line go out of the curve, but I want to end the line with the density point of the parameter. So, my problem is:

how can I find the values of the density at the mean, mode and median of my observations in order to set the correct coordinates for drawing?

Best Answer

The density is represented as a polyline, which is a pair of parallel arrays, one for $x$, one for $y$, forming vertices along the graph of the density (with equal spacings in the $x$ direction). As such it is a discrete approximation to the idealized continuous density and we can use discrete versions of the relevant integrals to compute statistics. Because the spacing is typically so close, there's probably little need to interpolate between successive points: we can use simple algorithms.

Whence,

x <- seq(-2.5, 10, length=1000000)
hx5 <- rnorm(x,0,1) + rexp(x,1/5) # tau=5 (rate = 1/tau)
#
# Compute the density.
#
dens <- density(hx5)
#
# Compute some measures of location.
#
n <- length(dens$y)                       #$
dx <- mean(diff(dens$x))                  # Typical spacing in x $
y.unit <- sum(dens$y) * dx                # Check: this should integrate to 1 $
dx <- dx / y.unit                         # Make a minor adjustment
x.mean <- sum(dens$y * dens$x) * dx
y.mean <- dens$y[length(dens$x[dens$x < x.mean])] #$
x.mode <- dens$x[i.mode <- which.max(dens$y)]
y.mode <- dens$y[i.mode]                  #$
y.cs <- cumsum(dens$y)                    #$
x.med <- dens$x[i.med <- length(y.cs[2*y.cs <= y.cs[n]])] #$
y.med <- dens$y[i.med]                                    #$
#
# Plot the density and the statistics.
#
plot(dens, xlim=c(-2.5,10), type="l", col="green",
     xlab="x", main="ExGaussian curve",lwd=2)
temp <- mapply(function(x,y,c) lines(c(x,x), c(0,y), lwd=2, col=c), 
               c(x.mean, x.med, x.mode), 
               c(y.mean, y.med, y.mode), 
               c("Blue", "Gray", "Red"))

Plot