7
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Have you considered the 'polygon()' function (stat.ethz.ch/R-manual/R-devel/library/graphics/html/…)? $\endgroup$ – ocram Jul 11 '12 at 10:34
  • 7
    $\begingroup$ The command density implements, by default, a Gaussian kernel with bandwidth given by bw.nrd0. You can then easily write this function: kg = function(x) return(mean(dnorm((x-hx5)/bw.nrd0(hx5)))/bw.nrd0(hx5)). Therefore you can evaluate the density at the desired points using e.g. kg(median(hx5)). $\endgroup$ – user10525 Jul 11 '12 at 10:56
  • $\begingroup$ @Procrastinator Your function works very well. Honestly, I thought it was more simple to find the density point of a parameter. Thanks! $\endgroup$ – this.is.not.a.nick Jul 11 '12 at 11:47
  • 2
    $\begingroup$ @Procrastinator - maybe you should repost as an answer? $\endgroup$ – jbowman Jul 11 '12 at 14:52
8
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ (Dollar signs in the code seem to toggle the indentation--must be an interaction with the support for $\TeX$ markup. Inserting extra dollar signs in comments to make an even number per line works around the problem!) $\endgroup$ – whuber Aug 1 '12 at 22:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.