50
$\begingroup$

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)?
Improve this question
$\endgroup$
4
  • 2
    $\begingroup$ What do you mean with overlapping region? Do you mean the area that is below both density curves? $\endgroup$
    – Nick Sabbe
    Jun 22, 2011 at 8:11
  • $\begingroup$ I mean the intersection of two areas $\endgroup$
    – Ali Salehi
    Jun 22, 2011 at 8:15
  • 5
    $\begingroup$ 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? $\endgroup$
    – whuber ♦
    Jun 22, 2011 at 14:38
  • $\begingroup$ See also: stats.stackexchange.com/questions/103800/… $\endgroup$
    – wolfies
    Jan 15, 2017 at 16:39

6 Answers 6

Reset to default
44
$\begingroup$

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

Improve this answer
$\endgroup$
2
  • 12
    $\begingroup$ (+1) Googling turns up at least three distinct definitions (Matsushita, Morisita, and Weitzman). Your implementation is Weitzman's. $\endgroup$
    – whuber ♦
    Jun 22, 2011 at 15:26
  • 1
    $\begingroup$ 0.60993 24 is an approximation for 0.60993 43398 78944 33895... . $\endgroup$
    – whuber ♦
    Jun 22, 2011 at 22:26
11
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • 3
    $\begingroup$ The Bhattacharyya coefficient is a measure of overlap but it is not the same, is it ? $\endgroup$
    – Stéphane Laurent
    Sep 14, 2013 at 14:00
7
$\begingroup$

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'.

Improve this answer
$\endgroup$
5
  • 4
    $\begingroup$ (+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. $\endgroup$
    – whuber ♦
    Jun 22, 2011 at 15:30
  • 2
    $\begingroup$ @whuber Could you turn this into a full answer? Or maybe Nick can edit his. $\endgroup$
    – Aleksandr Dubinsky
    Feb 10, 2013 at 9:02
  • $\begingroup$ @whuber Didn't you mean $\sigma_1 \geq \sigma_2$ rather than $\mu_1 \geq \mu_2$ ? $\endgroup$
    – Stéphane Laurent
    Sep 14, 2013 at 14:23
  • $\begingroup$ @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. $\endgroup$
    – whuber ♦
    Sep 14, 2013 at 19:21
  • $\begingroup$ @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. $\endgroup$
    – Stéphane Laurent
    Sep 14, 2013 at 23:51
1
$\begingroup$

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
}
Improve this answer
$\endgroup$
1
  • $\begingroup$ shouldn't it return (l+r)/2? $\endgroup$
    – RSHAP
    Apr 3, 2019 at 16:00
0
$\begingroup$

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;

}
Improve this answer
$\endgroup$
0
$\begingroup$

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.

Improve this answer
$\endgroup$

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