Porcentagem de regiões sobrepostas de duas distribuições normais

46

Fiquei me perguntando, dadas duas distribuições normais com eσ1, μ1σ2, μ2

  • como posso calcular a porcentagem de regiões sobrepostas de duas distribuições?
  • Suponho que esse problema tenha um nome específico. Você conhece algum nome específico que descreva esse problema?
  • Você está ciente de alguma implementação disso (por exemplo, código Java)?
Ali Salehi
fonte
2
O que você quer dizer com região sobreposta? Você quer dizer a área abaixo das duas curvas de densidade?
Nick Sabbe
Quero dizer a interseção de duas áreas
Ali Salehi
4
Em resumo, escrevendo os dois pdfs como e , você realmente deseja calcular ? Você poderia nos esclarecer sobre o contexto em que isso ocorre e como seria interpretado? fgmin(f(x),g(x))dx
whuber

Respostas:

41

Isso também é chamado de "coeficiente de sobreposição" (OVL). Ao pesquisar no Google, você obterá muitos hits. Você pode encontrar um nomograma para o caso bi-normal aqui . Um artigo útil pode ser:

  • Henry F. Inman; Edwin L. Bradley Jr (1989). O coeficiente de sobreposição como uma medida de concordância entre distribuições de probabilidade e estimativa pontual da sobreposição de duas densidades normais. Comunicações em Estatística - Teoria e Métodos, 18 (10), 3851-3874. ( Link )

Editar

Agora você me interessou mais por isso, então fui em frente e criei o código R para calcular isso (é uma integração simples). Joguei um gráfico das duas distribuições, incluindo o sombreamento da região sobreposta:

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)

Para este exemplo, o resultado é: 0.6099324com erro absoluto < 1e-04. Figura abaixo.

Exemplo

Wolfgang
fonte
10
(+1) O Google exibe pelo menos três definições distintas (Matsushita, Morisita e Weitzman). Sua implementação é da Weitzman.
whuber
1
0,60993 24 é uma aproximação para 0,60993 43398 78944 33895 ....
whuber
10

Isto é dado pelo coeficiente de Bhattacharyya . Para outras distribuições, consulte também a versão generalizada, a distância de Hellinger entre duas distribuições.

Não conheço nenhuma biblioteca para calcular isso, mas, dada a formulação explícita em termos das distâncias de Mahalanobis e das matrizes determinantes das variações, a implementação não deve ser um problema.

user603
fonte
3
O coeficiente de Bhattacharyya é uma medida de sobreposição, mas não é o mesmo, é?
Stéphane Laurent
7

Não sei se existe uma maneira óbvia de fazer isso, mas:

Primeiro, você encontra os pontos de interseção entre as duas densidades. Isso pode ser facilmente alcançado equacionando as duas densidades, que, para a distribuição normal, devem resultar em uma equação quadrática para x.

Algo próximo a:

(x-μ2)22σ22-(x-μ1)22σ12=registroσ1σ2

Isso pode ser resolvido com cálculo básico.

Assim, você tem zero, um ou dois pontos de interseção. Agora, esses pontos de interseção dividem a linha real em 1, 2 ou três partes, onde uma das duas densidades é a mais baixa. Se nada mais matemático vier à mente, tente qualquer ponto dentro de uma das partes para descobrir qual é a mais baixa.

Seu valor de interesse agora é a soma das áreas sob a curva de densidade mais baixa em cada parte. Agora, essa área pode ser encontrada na função de distribuição cumulativa (basta subtrair o valor nas duas arestas da 'peça'.

Nick Sabbe
fonte
4
(+1) Na verdade, quando , a equação pode ser resolvida com a fórmula quadrática: sem necessidade de cálculo. Se organizarmos (wlg) para , a segunda densidade será menor entre os dois zeros e, caso contrário, a primeira densidade será menor. Isso reduz o cálculo para quatro avaliações de um CDF normal. A situação com é ainda mais simples, exigindo solução de uma equação linear e apenas duas avaliações de um CDF. μ 1μ 2 σ 1 = σ 2σ1σ2μ1μ2σ1=σ2
whuber
2
@whuber Você poderia transformar isso em uma resposta completa? Ou talvez Nick possa editar o dele.
Aleksandr Dubinsky
@whuber Você não quis dizer vez de ? μ 1μ 2σ1σ2μ1μ2
Stéphane Laurent
@ Stéphane Eu acho que você está certo de que os SDs determinam a ordem: a densidade com SD menor acabará por ter caudas menores nas direções positiva e negativa e, portanto, terá os valores maiores entre os zeros e os valores menores em outros lugares.
whuber
@whuber Sim, e de fato é fácil ver que a ordem dos SDs determina o sinal do coeficiente de 2ª ordem do polinômio derivado de Nick.
Stéphane Laurent
1

Para a posteridade, a solução da wolfgang não funcionou para mim - encontrei bugs na integratefunção. Então eu combinei com a resposta de Nick Staubbe para desenvolver a pequena função a seguir. Deve ser mais rápido e com menos bugs do que usar integração numérica:

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
fonte
não deveria voltar (l+r)/2?
RSHAP
0

Aqui está a versão 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;

}
Vithun Venugopalan
fonte
0

Eu acho que algo assim poderia ser a solução no 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) 

Pelo menos eu poderia reproduzir o valor 0,8026 dado abaixo na Fig.1 neste pdf .

Você só precisa adaptar os valores inicial, final e de intervalo para ser preciso, pois essa é apenas uma solução numérica.

Danny K
fonte