Kernel density estimates - Add line for mean-CodePudding

I would like to add a line to my plot that indicates the mean:

X1 <- rnorm(100)
# Kernel density estimates
density_X1 <- density(X1)
# Compute mode
mode_X1 <- density_X1$x[which.max(density_X1$y)]; mode_X1
# Compute mean
mean_X1 <- mean(X1)
# Create plot 
plot(density_X1, main = "Kernel density estimates", xlab = "X1", ylab = "density")
# Add line for mode
segments(x0 = mode_X1, y0 = 0, x1 = mode_X1, y1 = density_X1$y, col = "red", lty = 1, lwd = 1)
# Add line for mean
segments(x0 = mean_X1, y0 = 0, x1 = mean_X1, y1 = , col = "red", lty = 5, lwd = 2)

What value do I need to enter for segments(y1 = ) to ensure that the upper bound is the kernel density plot?

CodePudding user response：

To calculate the y-value at some x-value from a density estimate at discrete points, you will need to do some sort of interpolation. One way is to use linear interpolation with approx:

approx(x = density_X1$x, y = density_X1$y, xout = mean_X1)$y

where, x and y give the points to be interpolated, and xout the (vector of) x-values where the interpolation is done (See ?approx)

So altogether the code is:

set.seed(72405277)
X1 <- rnorm(100)
# Kernel density estimates
density_X1 <- density(X1)
# Compute mode
mdIdx <- which.max(density_X1$y) # store and resuse
mode_X1 <- density_X1$x[mdIdx]; 
# Compute mean
mean_X1 <- mean(X1)
# Create plot 
plot(density_X1, main = "Kernel density estimates", xlab = "X1", ylab = "density")
# Add line for mode
segments(x0 = mode_X1, y0 = 0, x1 = mode_X1, y1 = density_X1$y[mdIdx], col = "red", lty = 1, lwd = 1)
# Add line for mean
yax2 <- approx(x = density_X1$x, y = density_X1$y, xout = mean_X1)$y
segments(x0 = mean_X1, y0 = 0, x1 = mean_X1, y1 = yax2, col = "blue", lty = 5, lwd = 2)

CodePudding user response：

You could also use abline like this:

X1 <- rnorm(100)
# Kernel density estimates
density_X1 <- density(X1)
# Compute mode
mode_X1 <- density_X1$x[which.max(density_X1$y)]; mode_X1
# Compute mean
mean_X1 <- mean(X1)
# Create plot 
plot(density_X1, main = "Kernel density estimates", xlab = "X1", ylab = "density")
abline(v = mode_X1, col = "red", lty = 1, lwd = 1)
abline(v = mean_X1, col = "red", lty = 5, lwd = 2)

Output: