How to calculate overlap between empirical probability densities?

I'm looking for a method to calculate the area of overlap between two kernel density estimates in R, as a measure of similarity between two samples. To clarify, in the following example, I would need to quantify the area of the purplish overlapping region:

library(ggplot2)
set.seed(1234)
d <- data.frame(variable=c(rep("a", 50), rep("b", 30)), value=c(rnorm(50), runif(30, 0, 3)))
ggplot(d, aes(value, fill=variable)) + geom_density(alpha=.4, color=NA) A similar question was discussed here, the difference being that I need to do this for arbitrary empirical data rather than predefined normal distributions. The overlap package addresses this question, but apparently only for timestamp data, which doesn't work for me. The Bray-Curtis index (as implemented in vegan package's vegdist(method="bray") function) also seems relevant but again for somewhat different data.

I'm interested in both the theoretical approach and the R functions I might employ to implement it.

• "quantify the purple area" is a problem in estimation, not in hypothesis testing, so you cannot hope to "accomplish this using a standard citable statistical test". You contradict yourself. Please clarify what you actually want. If all you want is an estimate of the area of overlap of two KDEs, that's a simple calculation. – Glen_b May 14 '14 at 0:39
• @Glen_b thanks for the comment, helped to clarify my non-statistician thinking. I believe the area of overlap between KDEs is indeed what I'm looking for -- I've edited the question to reflect that. – mmk May 14 '14 at 5:25
• I would be very concerned about the risk of arbitrariness in this method. Depending on the kernel bandwidth, the computed overlap between any two datasets could be made to equal any chosen value in the interval $(0,1)$. The default bandwidths are not optimized for this purpose and therefore conceivably could give surprising, arbitrary, or inconsistent results. Datasets with natural bounds (such as non-negative data or proportions, etc.) would further introduce unwanted edge effects. What to do instead? Start with the reason for this calculation: what is this "similarity" intended to mean? – whuber Jun 4 '15 at 19:26
• The same question appeared a few months later but referred to intersection points however there were some valid notes which could be taken into consideration. In the referred question is about two empirical distributions. I add the link as this post only answers this via kernel density estimation and for normal distributions. The link below I think extends on the question for pairs of empirical distributions. stats.stackexchange.com/questions/122857/… – Barnaby 7 hours ago – Barnaby Jul 27 '15 at 10:13

The area of overlap of two kernel density estimates may be approximated to any desired degree of accuracy.

1) Since the original KDEs have probably been evaluated over some grid, if the grid is the same for both (or can easily be made the same), the exercise could be as easy as simply taking $\min(K_1(x),K_2(x))$ at each point and then using the trapezoidal rule, or even a midpoint rule.

If the two are on different grids and can't easily be recalculated on the same grid, interpolation could be used.

2) You might find the point (or points) of intersection and integrate the lower of the two KDEs in each interval where each one is lower. In your diagram above you'd integrate the blue curve to the left of the intersection and the pink one to the right by whatever means you like/have available. This can be done essentially exactly by considering the area under each kernel component $\frac{1}{h}K(\frac{x-x_i}{h})$ to the left or right of that cut-off point.

However, whuber's comments above should be clearly borne in mind -- this is not necessarily a very meaningful thing to do.

• How do you calculate the error associated with method one and method 2? – olliepower May 14 '14 at 7:17
• In normal circumstances, both will be miniscule compared with the error in the kernel density estimates, so I wouldn't worry too much. Error bounds can be calculated on trapezoidal methods and other numerical integration of course - such calculations are pretty standard - but it's pointless worrying given that KDEs have large uncertainties. Method 2 will be accurate to accumulated rounding error of the calculations. – Glen_b May 14 '14 at 8:17
• These methodology suggestions makes sense, thanks very much for your answer. I will work on implementing this in R, but as a novice I would be interested in suggestions on how to code this cleanly. – mmk May 14 '14 at 17:21

For the sake of completeness, here's how I ended up doing this in R:

# simulate two samples
a <- rnorm(100)
b <- rnorm(100, 2)

# define limits of a common grid, adding a buffer so that tails aren't cut off
lower <- min(c(a, b)) - 1
upper <- max(c(a, b)) + 1

# generate kernel densities
da <- density(a, from=lower, to=upper)
db <- density(b, from=lower, to=upper)
d <- data.frame(x=da$x, a=da$y, b=db$y) # calculate intersection densities d$w <- pmin(d$a, d$b)

# integrate areas under curves
library(sfsmisc)
total <- integrate.xy(d$x, d$a) + integrate.xy(d$x, d$b)
intersection <- integrate.xy(d$x, d$w)

# compute overlap coefficient
overlap <- 2 * intersection / total


As noted, there is inherent uncertainty and subjectivity involved in the KDE generation and also in the integration.

First, I might be wrong but I think your solution wouldn't work in case where there is multiples points where the Kernel Density Estimates (KDE) intersect. Second, although the overlap package was created for use with timestamp data, you can still use it to estimate the area of overlap of any two KDEs. You simply have to rescale your data so that it range from 0 to 2π.
For exemple :

# simulate two sample
a <- rnorm(100)
b <- rnorm(100, 2)

# To use overplapTrue(){overlap} the scale must be in radian (i.e. 0 to 2pi)
# To keep the *relative* value of a and b the same, combine a and b in the
# same dataframe before rescaling. You'll need to load the ‘scales‘ library.
# But first add a "Source" column to be able to distinguish between a and b
# after they are combined.
a = data.frame( value = a, Source = "a" )
b = data.frame( value = b, Source = "b" )
d = rbind(a, b)
library(scales)
d$value <- rescale( d$value, to = c(0,2*pi) )

# Now you can created the rescaled a and b vectors
a <- d[d$Source == "a", 1] b <- d[d$Source == "b", 1]

# You can then calculate the area of overlap as you did previously.
# It should give almost exactly the same answers.
# Or you can use either the overlapTrue() and overlapEst() function
# provided with the overlap packages.
# Note that with these function the KDE are fitted using von Mises kernel.
library(overlap)
# Using overlapTrue():
# define limits of a common grid, adding a buffer so that tails aren't cut off
lower <- min(d$value)-1 upper <- max(d$value)+1
# generate kernel densities
da <- density(a, from=lower, to=upper, adjust = 1)
db <- density(b, from=lower, to=upper, adjust = 1)
# Compute overlap coefficient
overlapTrue(da$y,db$y)

# Using overlapEst():
overlapEst(a, b, kmax = 3, adjust=c(0.8, 1, 4), n.grid = 500)

# You can also plot the two KDEs and the region of overlap using overlapPlot()
# but sadly I haven't found a way of changing the x scale so that the scale
# range correspond to the initial x value and not the rescaled value.
# You can only change the maximum value of the scale using the xscale argument
# (i.e. it always range from 0 to n, where n is set with xscale = n).
# So if some of your data take negative value, you're probably better off with
# a different plotting method. You can change the x label with the xlab
# argument.
overlapPlot(a, b, xscale = 10, xlab= "x metrics", rug=T)