Procent nakładających się regionów dwóch rozkładów normalnych

Zastanawiałem się, biorąc pod uwagę dwie normalne dystrybucje z iσ 2 , μ $\sigma_1,\ \mu_1$$\sigma_2, \ \mu_2$

• jak mogę obliczyć procent nakładających się regionów dwóch rozkładów?
• Podejrzewam, że ten problem ma konkretną nazwę. Czy znasz jakieś konkretne nazwy opisujące ten problem?
• Czy znasz jakieś implementacje tego (np. Kod Java)?

Jest to często nazywane „współczynnikiem nakładania się” (OVL). Googlowanie za to da ci wiele trafień. Można znaleźć nomogramu dla bi-normalnym przypadku tutaj . Przydatnym papierem może być:

• Henry F. Inman; Edwin L. Bradley Jr (1989). Współczynnik nakładania się jako miara zgodności między rozkładami prawdopodobieństwa i oszacowanie punktowe nakładania się dwóch normalnych gęstości. Communications in Statistics - Theory and Methods, 18 (10), 3851-3874. ( Link )

Edytować

Teraz zainteresowałeś mnie tym bardziej, więc stworzyłem kod R, aby to obliczyć (jest to prosta integracja). Wrzuciłem wykres dwóch rozkładów, w tym cieniowanie nakładającego się regionu:

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)

W tym przykładzie wynik jest następujący: 0.6099324z błędem bezwzględnym < 1e-04. Niżej wymienione.

Daje to współczynnik Bhattacharyya . W przypadku innych dystrybucji zobacz także wersję uogólnioną, odległość Hellingera między dwiema dystrybucjami.

Nie znam żadnych bibliotek do obliczenia tego, ale biorąc pod uwagę wyraźne sformułowanie w odniesieniu do odległości Mahalanobisa i wyznacznika macierzy wariancji, wdrożenie nie powinno stanowić problemu.

Nie wiem, czy istnieje oczywisty standardowy sposób na zrobienie tego, ale:

Najpierw znajdź punkty przecięcia dwóch gęstości. Można to łatwo osiągnąć przez zrównanie obu gęstości, co przy rozkładzie normalnym powinno skutkować równaniem kwadratowym dla x.

Coś w pobliżu:

Można to rozwiązać za pomocą rachunku różniczkowego.

Zatem masz zero, jeden lub dwa punkty przecięcia. Teraz te punkty przecięcia dzielą rzeczywistą linię na 1, 2 lub trzy części, gdzie jedna z dwóch gęstości jest najniższa. Jeśli nie przychodzi ci na myśl nic więcej matematycznego, po prostu wypróbuj dowolny punkt w jednej z części, aby dowiedzieć się, który z nich jest najniższy.

Twoja wartość zainteresowania jest teraz sumą obszarów pod krzywą najniższej gęstości w każdej części. Obszar ten można teraz znaleźć w funkcji skumulowanego rozkładu (wystarczy odjąć wartość na obu krawędziach „części”).

Dla potomnych rozwiązanie Wolfganga nie zadziałało - natknąłem się na błędy w integratefunkcji. Połączyłem to z odpowiedzią Nicka Staubbe, aby opracować następującą małą funkcję. Powinny być szybsze i mniej obciążające niż przy użyciu integracji numerycznej:

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
}

Oto wersja Java, 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;

}

Myślę, że coś takiego może być rozwiązaniem w 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) 

Przynajmniej mógłbym odtworzyć wartość 0,8026 podaną poniżej na ryc. 1 w tym pliku pdf .

Musisz tylko precyzyjnie dostosować wartości początkową i końcową oraz wartości przedziałów, ponieważ jest to tylko rozwiązanie numeryczne.