I was wondering, given two normal distributions with $\sigma_1,\ \mu_1$ and $\sigma_2, \ \mu_2$

  • how can I calculate the percentage of overlapping regions of two distributions?
  • I suppose this problem has a specific name, are you aware of any particular name describing this problem?
  • Are you aware of any implementation of this (e.g., Java code)?
  • 2
    What do you mean with overlapping region? Do you mean the area that is below both density curves? – Nick Sabbe Jun 22 '11 at 8:11
  • I mean the intersection of two areas – Ali Salehi Jun 22 '11 at 8:15
  • 4
    In short, writing the two pdfs as $f$ and $g$, do you really want to compute $\int \min(f(x),g(x))dx$ ? Could you enlighten us about the context in which this arises and how it would be interpreted? – whuber Jun 22 '11 at 14:38
  • See also: stats.stackexchange.com/questions/103800/… – wolfies Jan 15 '17 at 16:39
up vote 40 down vote accepted

This is also often called the "overlapping coefficient" (OVL). Googling for this will give you lots of hits. You can find a nomogram for the bi-normal case here. A useful paper may be:

  • Henry F. Inman; Edwin L. Bradley Jr (1989). The overlapping coefficient as a measure of agreement between probability distributions and point estimation of the overlap of two normal densities. Communications in Statistics - Theory and Methods, 18(10), 3851-3874. (Link)

Edit

Now you got me interested in this more, so I went ahead and created R code to compute this (it's a simple integration). I threw in a plot of the two distributions, including the shading of the overlapping region:

min.f1f2 <- function(x, mu1, mu2, sd1, sd2) {
    f1 <- dnorm(x, mean=mu1, sd=sd1)
    f2 <- dnorm(x, mean=mu2, sd=sd2)
    pmin(f1, f2)
}

mu1 <- 2;    sd1 <- 2
mu2 <- 1;    sd2 <- 1

xs <- seq(min(mu1 - 3*sd1, mu2 - 3*sd2), max(mu1 + 3*sd1, mu2 + 3*sd2), .01)
f1 <- dnorm(xs, mean=mu1, sd=sd1)
f2 <- dnorm(xs, mean=mu2, sd=sd2)

plot(xs, f1, type="l", ylim=c(0, max(f1,f2)), ylab="density")
lines(xs, f2, lty="dotted")
ys <- min.f1f2(xs, mu1=mu1, mu2=mu2, sd1=sd1, sd2=sd2)
xs <- c(xs, xs[1])
ys <- c(ys, ys[1])
polygon(xs, ys, col="gray")

### only works for sd1 = sd2
SMD <- (mu1-mu2)/sd1
2 * pnorm(-abs(SMD)/2)

### this works in general
integrate(min.f1f2, -Inf, Inf, mu1=mu1, mu2=mu2, sd1=sd1, sd2=sd2)

For this example, the result is: 0.6099324 with absolute error < 1e-04. Figure below.

Example

  • 10
    (+1) Googling turns up at least three distinct definitions (Matsushita, Morisita, and Weitzman). Your implementation is Weitzman's. – whuber Jun 22 '11 at 15:26
  • 1
    0.60993 24 is an approximation for 0.60993 43398 78944 33895... . – whuber Jun 22 '11 at 22:26

This is given by the Bhattacharyya coefficient. For other distributions, see also the generalised version, the Hellinger distance between two distributions.

I don't know of any libraries to compute this, but given the explicit formulation in terms of Mahalanobis distances and determinant of variance matrices, implementation should not be an issue.

  • 3
    The Bhattacharyya coefficient is a measure of overlap but it is not the same, is it ? – Stéphane Laurent Sep 14 '13 at 14:00

I don't know if there is an obvious standard way of doing this, but:

First, you find the intersection points between the two densities. This can be easily achieved by equating both densities, which, for the normal distribution, should result in a quadratic equation for x.

Something close to: $$ \frac{(x-\mu_2)^2}{2\sigma_2^2} - \frac{(x-\mu_1)^2}{2\sigma_1^2} = \log{\frac{\sigma_1}{\sigma_2}} $$

This can be solved with basic calculus.

Thus you have either zero, one or two intersection points. Now, these intersection points divide the real line into 1, 2 or three parts, where either of the two densities is the lowest one. If nothing more mathematical comes to mind, just try any point within one of the parts to find which one is the lowest.

Your value of interest is now the sum of the areas under the lowest density curve in each part. This area can now be found from the cumulative distribution function (just subtract the value in both edges of the 'part'.

  • 4
    (+1) Actually, when $\sigma_1 \ne \sigma_2$, the equation can be solved with the quadratic formula: no need for calculus. If we arrange (wlg) for $\mu_1 \ge \mu_2$, then the second density is smallest between the two zeros and otherwise the first density is smallest. This reduces the calculation to four evaluations of a Normal CDF. The situation with $\sigma_1 = \sigma_2$ is even simpler, requiring solution of a linear equation and only two evaluations of a CDF. – whuber Jun 22 '11 at 15:30
  • 1
    @whuber Could you turn this into a full answer? Or maybe Nick can edit his. – Aleksandr Dubinsky Feb 10 '13 at 9:02
  • @whuber Didn't you mean $\sigma_1 \geq \sigma_2$ rather than $\mu_1 \geq \mu_2$ ? – Stéphane Laurent Sep 14 '13 at 14:23
  • @Stéphane I think you're correct that the SDs determine the order: the density with smaller SD will eventually have smaller tails in both the positive and negative directions and therefore will have the larger values between the zeros and smaller values elsewhere. – whuber Sep 14 '13 at 19:21
  • @whuber Yes, and indeed it is easy to see that the order of the SDs determine the sign of the 2nd order coefficient of the polynom derived by Nick. – Stéphane Laurent Sep 14 '13 at 23:51

For posterity, wolfgang's solution didn't work for me—I ran into bugs in the integrate function. So I combined it with Nick Staubbe's answer to develop the following little function. Should be quicker and less buggy than using numerical integration:

get_overlap_coef <- function(mu1, mu2, sd1, sd2){
  xs  <- seq(min(mu1 - 4*sd1, mu2 - 4*sd2), 
             max(mu1 + 4*sd1, mu2 + 4*sd2), 
             length.out = 500)
  f1  <- dnorm(xs, mean=mu1, sd=sd1)
  f2  <- dnorm(xs, mean=mu2, sd=sd2)
  int <- xs[which.max(pmin(f1, f2))]
  l   <- pnorm(int, mu1, sd1, lower.tail = mu1>mu2)
  r   <- pnorm(int, mu2, sd2, lower.tail = mu1<mu2)
  l+r
}

Here is the Java version, Apache Commons Mathematics Library:

import org.apache.commons.math3.distribution.NormalDistribution;

public static double overlapArea(double mean1, double sd1, double mean2, double sd2) {

    NormalDistribution normalDistribution1 = new NormalDistribution(mean1, sd1);
    NormalDistribution normalDistribution2 = new NormalDistribution(mean2, sd2);

    double min = Math.min(mean1 - 6 * sd1, mean2 - 6 * sd2);
    double max = Math.max(mean1 + 6 * sd1, mean2 + 6 * sd2);
    double range = max - min;

    int resolution = (int) (range/Math.min(sd1, sd2));

    double partwidth = range / resolution;

    double intersectionArea = 0;

    int begin = (int)((Math.max(mean1 - 6 * sd1, mean2 - 6 * sd2)-min)/partwidth);
    int end = (int)((Math.min(mean1 + 6 * sd1, mean2 + 6 * sd2)-min)/partwidth);

    /// Divide the range into N partitions
    for (int ii = begin; ii < end; ii++) {

        double partMin = partwidth * ii;
        double partMax = partwidth * (ii + 1);

        double areaOfDist1 = normalDistribution1.probability(partMin, partMax);
        double areaOfDist2 = normalDistribution2.probability(partMin, partMax);

        intersectionArea += Math.min(areaOfDist1, areaOfDist2);
    }

    return intersectionArea;

}

I think something like this could be the solution in MATLAB:

[overlap] = calc_overlap_twonormal(2,2,0,1,-20,20,0.01)

% numerical integral of the overlapping area of two normal distributions:
% s1,s2...sigma of the normal distributions 1 and 2
% mu1,mu2...center of the normal distributions 1 and 2
% xstart,xend,xinterval...defines start, end and interval width
% example: [overlap] = calc_overlap_twonormal(2,2,0,1,-10,10,0.01)

function [overlap2] = calc_overlap_twonormal(s1,s2,mu1,mu2,xstart,xend,xinterval)

clf
x_range=xstart:xinterval:xend;
plot(x_range,[normpdf(x_range,mu1,s1)' normpdf(x_range,mu2,s2)']);
hold on
area(x_range,min([normpdf(x_range,mu1,s1)' normpdf(x_range,mu2,s2)']'));
overlap=cumtrapz(x_range,min([normpdf(x_range,mu1,s1)' normpdf(x_range,mu2,s2)']'));
overlap2 = overlap(end);

[overlap] = calc_overlap_twonormal(2,2,0,1,-10,10,0.01) 

At least I could reproduce the value 0.8026 given below Fig.1 in this pdf.

You just need to adapt the start and end and interval values to be precise as this is only a numerical solution.

protected by kjetil b halvorsen Nov 1 '17 at 19:14

Thank you for your interest in this question. Because it has attracted low-quality or spam answers that had to be removed, posting an answer now requires 10 reputation on this site (the association bonus does not count).

Would you like to answer one of these unanswered questions instead?

Not the answer you're looking for? Browse other questions tagged or ask your own question.