First, calculate both densities in their united range u.
u <- range(c(x, y))
dx <- density(x, from=u[1], to=u[2])
dy <- density(y, from=u[1], to=u[2])
Second subtract the estimations yfrom each other.
dd_xy <- dx$y - dy$y
The x sould be the same.
stopifnot(all.equal(dx$x, dy$x))
Plot
Then plot one density and use lines to add the others.
plot(dx, col=4, ylim=c(-.25, .45), main='Density distributions', xlab='')
abline(h=0, lty=3, col=8)
lines(dy, col=3)
lines(dx$x, dd_xy, col=2, lty=2, lwd=2) ## <---------------- difference
mtext(sprintf('N(x) = %s Bandwidth(x) = %s', dx$n, signif(dx$bw, 3)), 1, 2)
mtext(sprintf('N(y) = %s Bandwidth(y) = %s', dy$n, signif(dy$bw, 3)), 1, 3)
legend('topleft', legend=c('x', 'y', 'x - y'), col=4:2,
lty=c(1, 1, 2), lwd=c(1, 1, 2), title='density')
![enter image description here]()
Calculations
sapply(c('sum', 'mean', 'sd', 'min', 'max'), \(x) do.call(x, list(dd_xy))) |>
signif(3)
# sum mean sd min max
# -0.049700 -0.000097 0.088000 -0.249000 0.122000
Data:
set.seed(42)
x <- rnorm(1000, 0, 2)
y <- rnorm(500, 1, 1)