# identify the point of intersection from two distributions

Dear StackExchange community,

I have a problem of identifying where two distributions $$F$$ and $$G$$ intersect (cross each other). In particular, I have an empirical estimation of $$F$$ and $$G$$ from the data I have and I'm looking for a point above which $$G$$ is likely true over $$F$$. In other words, $$F$$ is a reference distribution and $$G$$ is the target. Therefore, we would like to know if an item with value $$x$$ is positive (which means that $$x$$ has a higher probability under $$G$$ than in $$F$$). The problem is that the empirical estimations of $$F$$ and $$G$$ result in multi-modal distribution, and hence it has been a difficult task for me to obtain one most plausible intersection point. Please also see the attached image for an example (a simple scenario). Note that the empirical estimation was obtained using kernel density estimation (density() of R). se let me know if there is a method I can try to obtain the intersection point.

• How are these curves obtained? You say 'empirical distributions' (a phrase which I'd normally assume meant an ECDF, which is not what we're looking at); are these kernel density estimates, or something else? Aug 27, 2022 at 17:09
• The curves are estimated using the density() function of R Aug 27, 2022 at 17:29
• If the desired precision for the intersection point(s) is just to produce a figure, then use the same from and to options of density and linearly interpolate the locations where the differences in density change sign.
– JimB
Aug 27, 2022 at 19:24
• Could you explain why you want to find these intersection point(s)? What do you suppose they might represent? Because it is rare for this procedure to be of any statistical interest, maybe what you really need is a different procedure altogether. Finding them, btw, is an algorithmic matter as sketched by @JimB. The algorithm is simple: subtract one set of density values from the other into an array d and locate indexes i where d[i] * d[i+1] <= 0. Linear interpolation between the abscissae x[i] and x[i+1] does the trick. You should always find at least two intersections.
– whuber
Aug 27, 2022 at 20:07
• Perhaps I am being overly pedantic, but what comes to mind when I hear the "intersection from two distributions" is the set of points that are in the support of both distributions. I think that something like "crossover point of the PDFs" would be clearer. Aug 29, 2022 at 1:47

A tiny bit of statistics is needed here, if only to point out the need to control the bandwidth and study the sensitivity of the solutions to the bandwidth.

Provided a solution is bracketed closely by data points, it will tend to be stable even when the bandwidth is varied substantially. Here is an example involving datasets with 23 points (black density) and 14 points (light blue). Red vertical lines mark the solutions. The data are shown as rug plots at the bottom. The default bandwidth for these data will be around $$1/2,$$ as shown in the middle panel

You can see from this example how one solution (the right hand one in the right panel) persists across all bandwidths. Another solution (the left hand one in the right panel) varies appreciably because data are scarce in its neighborhood. Spurious solutions pop up when using a relatively small bandwidth (left panel).

These examples were created by this R code.

set.seed(17)
x <- rnorm(23)
y <- rnorm(14, 2, 3/2)
bw <- 0.25 # or 0.5, or 1.5, or even "SP", etc: see the help page for density
obj <- intersect(x, y, kernel = "gaussian", n = 512, bw = bw, from = -4, to = 8)


All kernel densities produce a discrete grid of density estimates. The solution implemented by intersect allows you to exploit the default methods of finding endpoints, bandwidths, etc by first computing a density for the combined data. Those defaults are then used to recompute the densities for the data separately. Because both densities are computed on the same grid, it's a simple matter to locate the places where they cross and interpolate linearly on the grid. Linear interpolation is more than precise enough, because it errs less than the mesh of the grid, which presumably is already small enough for your purposes.

#
# Find all points where density $g$ exceeds density $f.$
#
intersect <- function(x, y, bw = "nrd0", from, to,  ...) {
#
# Compute a density for all points combined.
#
largs <- list(x = c(x,y), bw = bw)
if (!missing(from)) largs <- c(largs, from = from)
if (!missing(to)) largs <- c(largs, to = to)
largs <- c(largs, list(...))
obj <- do.call(density, largs) # Compute a common density
#
# Compute densities for the datasets separately.
#
x.0 <- obj$x f.x <- density(x, bw = obj$bw, from = min(x.0), to = max(x.0), ...)
f.y <- density(y, bw = obj$bw, from = min(x.0), to = max(x.0), ...) # # Find the crossings. # d <- zapsmall(f.y$y - f.x\$y)
abscissae <- sapply(which(d[-1] * d[-length(d)] < 0), function(i) {
w <- d[i+1] - d[i]
if (w > 0) (d[i+1] * x.0[i] - d[i] * x.0[i+1]) / w else (x.0[i] + x.0[i+1]) / 2
})
list(Points = abscissae, xlim = range(x.0), f = f.x, g = f.y)
}

• Much appreciation! This is the solution I'm looking for. Aug 28, 2022 at 17:43

Here is a crude approach to find the intersection point(s).

# Generate data
set.seed(12345)
x <- rnorm(100)
y <- rnorm(150, 1, 3)

# Find global minimum and maximum
xymin <- min(x,y)
xymax <- max(x,y)
# Estimate densities
dx <- density(x, n=512, from=xymin, to=xymax)
dy <- density(y, n=512, from=xymin, to=xymax)

# Plot results
plot(dx, xlim=c(xymin, xymax), type="l", lwd=3, xlab="X", ylab="Density", main="")
lines(dy, col="red", lwd=3)

# Differences in densities
dx$$diff <- dx$$y - dy$$y ex <- NULL # Store the interection points ey <- NULL k = 0 for (i in 2:length(dx$$x)) {
# Look for a change in sign of the difference in densities
if (sign(dx$$diff[i-1]) != sign(dx$$diff[i])) {
k = k + 1
# Linearly interpolate
ex[k] <- dx$$x[i-1] + (dx$$x[i]-dx$$x[i-1])*(0-dx$$diff[i-1])/(dx$$diff[i]-dx$$diff[i-1])
ey[k] <- dx$$y[i-1] + (dx$$y[i]-dx$$y[i-1])*(ex[k]-dx$$x[i-1])/(dx$$x[i]-dx$$x[i-1])
lines(c(ex[k],ex[k]), c(0,ey[k]))
points(ex[k], ey[k], pch=16, col="green" )
}
} cbind(ex, ey)

#            ex         ey
#[1,] -1.957378 0.06736659
#[2,]  2.106521 0.12664663

• Much appreciation! I loved this method as well. Aug 28, 2022 at 17:44
• Good. You're question was very clear as to "what" you wanted, but it would be interesting to know "why" you want to do this (as asked in a comment).
– JimB
Aug 28, 2022 at 17:53