identify the point of intersection from two distributions

Question

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.

Thanks in advance.

Answer 1

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)
}

Answer 2

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