Prozentsatz der überlappenden Bereiche zweier Normalverteilungen

46

Ich wunderte mich, zwei Normalverteilungen mit undσ1, μ1σ2, μ2

  • Wie kann ich den Prozentsatz überlappender Bereiche zweier Verteilungen berechnen?
  • Ich nehme an, dieses Problem hat einen bestimmten Namen. Kennen Sie einen bestimmten Namen, der dieses Problem beschreibt?
  • Ist Ihnen eine Implementierung davon bekannt (z. B. Java-Code)?
Ali Salehi
quelle
2
Was meinst du mit überlappender Region? Meinen Sie den Bereich, der unter beiden Dichtekurven liegt?
Nick Sabbe
Ich meine den Schnittpunkt zweier Gebiete
Ali Salehi
4
Kurz gesagt, wenn Sie die beiden pdfs als und schreiben , möchten Sie wirklich berechnen ? Können Sie uns erklären, in welchem ​​Kontext dies geschieht und wie es interpretiert wird? fGMindest(f(X),G(X))dX
Whuber

Antworten:

41

Dies wird auch oft als "Überlappungskoeffizient" (OVL) bezeichnet. Wenn du dafür googelst, erhältst du viele Treffer. Ein Nomogramm für den Bi-Normalfall finden Sie hier . Ein nützliches Papier kann sein:

  • Henry F. Inman; Edwin L. Bradley Jr. (1989). Der Überlappungskoeffizient als Maß für die Übereinstimmung zwischen Wahrscheinlichkeitsverteilungen und der Punktschätzung der Überlappung zweier normaler Dichten. Kommunikationen in der Statistik - Theorie und Methoden, 18 (10), 3851-3874. ( Link )

Bearbeiten

Jetzt hast du mich mehr dafür interessiert, also habe ich R-Code erstellt, um dies zu berechnen (es ist eine einfache Integration). Ich habe eine Darstellung der beiden Verteilungen, einschließlich der Schattierung des überlappenden Bereichs, eingefügt:

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)

Für dieses Beispiel lautet das Ergebnis: 0.6099324mit absolutem Fehler < 1e-04. Abbildung unten.

Beispiel

Wolfgang
quelle
10
(+1) Beim Googeln werden mindestens drei unterschiedliche Definitionen angezeigt (Matsushita, Morisita und Weitzman). Ihre Implementierung ist die von Weitzman.
whuber
1
0,60993 24 ist eine Näherung für 0,60993 43398 78944 33895 ....
Whuber
10

Dies ist durch den Bhattacharyya-Koeffizienten gegeben . Für andere Distributionen siehe auch die verallgemeinerte Version, der Hellinger-Abstand zwischen zwei Distributionen.

Ich kenne keine Bibliotheken, um dies zu berechnen, aber angesichts der expliziten Formulierung in Bezug auf Mahalanobis-Abstände und Determinante von Varianzmatrizen sollte die Implementierung kein Problem sein.

user603
quelle
3
Der Bhattacharyya-Koeffizient ist ein Maß für die Überlappung, aber es ist nicht dasselbe, oder?
Stéphane Laurent
7

Ich weiß nicht, ob es dafür einen offensichtlichen Standard gibt, aber:

Zunächst finden Sie die Schnittpunkte zwischen den beiden Dichten. Dies kann leicht erreicht werden, indem beide Dichten gleichgesetzt werden, was für die Normalverteilung zu einer quadratischen Gleichung für x führen sollte.

(X-μ2)22σ22-(X-μ1)22σ12=Logσ1σ2

Dies kann mit der Grundrechnung gelöst werden.

Sie haben also entweder Null, einen oder zwei Schnittpunkte. Diese Schnittpunkte teilen nun die reale Linie in 1, 2 oder drei Teile, wobei eine der beiden Dichten die niedrigste ist. Wenn Ihnen nichts Mathematischeres einfällt, versuchen Sie einfach einen beliebigen Punkt in einem der Teile, um herauszufinden, welcher der niedrigste ist.

Ihr interessierender Wert ist jetzt die Summe der Bereiche unter der Kurve mit der niedrigsten Dichte in jedem Teil. Dieser Bereich kann nun über die kumulative Verteilungsfunktion ermittelt werden (subtrahieren Sie einfach den Wert an beiden Kanten des 'Teils'.

Nick Sabbe
quelle
4
σ1σ2μ1μ2σ1=σ2
2
@whuber Könntest du daraus eine vollständige Antwort machen? Oder vielleicht kann Nick seine bearbeiten.
Aleksandr Dubinsky
σ1σ2μ1μ2
@Stéphane Ich denke, Sie haben Recht, dass die SDs die Reihenfolge bestimmen: Die Dichte mit kleinerem SD wird eventuell kleinere Schwänze sowohl in der positiven als auch in der negativen Richtung haben und daher die größeren Werte zwischen den Nullen und die kleineren Werte an anderer Stelle haben.
whuber
@whuber Ja, und in der Tat ist leicht zu erkennen, dass die Reihenfolge der SDs das Vorzeichen des Koeffizienten 2. Ordnung des von Nick abgeleiteten Polynoms bestimmt.
Stéphane Laurent
1

Für die Nachwelt hat die Lösung von wolfgang bei mir nicht funktioniert - ich bin auf Fehler in der integrateFunktion gestoßen. Also habe ich es mit der Antwort von Nick Staubbe kombiniert, um die folgende kleine Funktion zu entwickeln. Sollte schneller und weniger fehleranfällig sein als die numerische 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
}
generic_user
quelle
sollte es nicht zurückkehren (l+r)/2?
RSHAP
0

Hier ist die 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;

}
Vithun Venugopalan
quelle
0

Ich denke, so etwas könnte die Lösung in MATLAB sein:

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

Zumindest konnte ich den Wert 0,8026, der unten in Abb.1 angegeben ist, in diesem PDF reproduzieren .

Sie müssen nur Start- und End- sowie Intervallwerte anpassen, um genau zu sein, da dies nur eine numerische Lösung ist.

Danny K
quelle