Como obter a região da elipse a partir de dados distribuídos normais bivariados?

11

Eu tenho dados parecidos com:

Figura

Tentei aplicar a distribuição normal (a estimativa de densidade do kernel funciona melhor, mas não preciso de uma precisão tão grande) nela e funciona muito bem. O gráfico de densidade faz uma elipse.

Eu preciso que a função elipse decida se um ponto está dentro da região da elipse ou não. Como fazer isso?

Código R ou Mathematica são bem-vindos.

matejuh
fonte

Respostas:

18

O Corsario fornece uma boa solução em um comentário: use a função de densidade do kernel para testar a inclusão em um conjunto de níveis.

Outra interpretação da questão é que ela solicita um procedimento para testar a inclusão nas elipses criadas por uma aproximação normal bivariada dos dados. Para começar, vamos gerar alguns dados que se parecem com a ilustração na pergunta:

library(mvtnorm) # References rmvnorm()
set.seed(17)
p <- rmvnorm(1000, c(250000, 20000), matrix(c(100000^2, 22000^2, 22000^2, 6000^2),2,2))

As elipses são determinadas pelo primeiro e segundo momentos dos dados:

center <- apply(p, 2, mean)
sigma <- cov(p)

A fórmula requer inversão da matriz de variância-covariância:

sigma.inv = solve(sigma, matrix(c(1,0,0,1),2,2))

A função "altura" da elipse é negativa do logaritmo da densidade normal bivariada :

ellipse <- function(s,t) {u<-c(s,t)-center; u %*% sigma.inv %*% u / 2}

(Eu ignorei uma constante aditiva igual a .)log(2πdet(Σ))

Para testar isso , vamos desenhar alguns de seus contornos. Isso requer a geração de uma grade de pontos nas direções x e y:

n <- 50
x <- (0:(n-1)) * (500000/(n-1))
y <- (0:(n-1)) * (50000/(n-1))

Calcule a função de altura nesta grade e plote-a:

z <- mapply(ellipse, as.vector(rep(x,n)), as.vector(outer(rep(0,n), y, `+`)))
plot(p, pch=20, xlim=c(0,500000), ylim=c(0,50000), xlab="Packets", ylab="Flows")
contour(x,y,matrix(z,n,n), levels=(0:10), col = terrain.colors(11), add=TRUE)

Gráfico de contorno

Evidentemente, funciona. Portanto, o teste para determinar se um ponto está dentro de um contorno elíptico no nível é(s,t)c

ellipse(s,t) <= c

O Mathematica faz o trabalho da mesma maneira: calcule a matriz de variância-covariância dos dados, inverta isso, construa a ellipsefunção e está tudo pronto.

whuber
fonte
Obrigado a todos, especialmente @whuber. É exatamente disso que eu preciso.
matejuh
Btw. existe alguma solução simples para contornos de estimativa de densidade de kernel? Porque se eu quiser ser mais rigoroso, meus dados se parecem com: github.com/matejuh/doschecker_wiki_images/raw/master/… resp. github.com/matejuh/doschecker_wiki_images/raw/master/…
matejuh
Não consigo encontrar uma solução simples em R. Considere usar a função "SmoothKernelDistribution" do Mathematica 8.
whuber
2
Os níveis correspondem ao nível de confiança? Acho que não. Como posso fazer isso, por favor?
matejuh
Isso precisa de uma nova pergunta, porque você precisa especificar do que busca a confiança e - a julgar pelas suas plotagens - há preocupações sobre se essas elipses são descrições adequadas dos dados em primeiro lugar.
whuber
9

A plotagem é direta com a ellipse()função do mixtoolspacote para R:

library(mixtools)
library(mvtnorm) 
set.seed(17)
p <- rmvnorm(1000, c(250000, 20000), matrix(c(100000^2, 22000^2, 22000^2, 6000^2),2,2))
plot(p, pch=20, xlim=c(0,500000), ylim=c(0,50000), xlab="Packets", ylab="Flows")
ellipse(mu=colMeans(p), sigma=cov(p), alpha = .05, npoints = 250, col="red") 

insira a descrição da imagem aqui

Stéphane Laurent
fonte
5

Primeira abordagem

Você pode tentar essa abordagem no Mathematica.

Vamos gerar alguns dados bivariados:

data = Table[RandomVariate[BinormalDistribution[{50, 50}, {5, 10}, .8]], {1000}];

Então precisamos carregar este pacote:

Needs["MultivariateStatistics`"]

E agora:

ellPar=EllipsoidQuantile[data, {0.9}]

fornece uma saída que define uma elipse de confiança de 90%. Os valores que você obtém dessa saída estão no seguinte formato:

{Ellipsoid[{x1, x2}, {r1, r2}, {{d1, d2}, {d3, d4}}]}

x1 e x2 especificam o ponto no qual a elipse no centro, r1 e r2 especificam os raios do semi-eixo, e d1, d2, d3 e d4 especificam a direção do alinhamento.

Você também pode traçar isso:

Show[{ListPlot[data, PlotRange -> {{0, 100}, {0, 100}}, AspectRatio -> 1],  Graphics[EllipsoidQuantile[data, 0.9]]}]

A forma paramétrica geral da elipse é:

ell[t_, xc_, yc_, a_, b_, angle_] := {xc + a Cos[t] Cos[angle] - b Sin[t] Sin[angle],
    yc + a Cos[t] Sin[angle] + b Sin[t] Cos[angle]}

E você pode plotá-lo desta maneira:

ParametricPlot[
    ell[t, ellPar[[1, 1, 1]], ellPar[[1, 1, 2]], ellPar[[1, 2, 1]], ellPar[[1, 2, 2]],
    ArcTan[ellPar[[1, 3, 1, 2]]/ellPar[[1, 3, 1, 1]]]], {t, 0, 2 \[Pi]},
    PlotRange -> {{0, 100}, {0, 100}}]

Você pode executar uma verificação com base em informações geométricas puras: se a distância euclidiana entre o centro da elipse (ellPar [[1,1]]) e seu ponto de dados for maior que a distância entre o centro da elipse e a borda da elipse a elipse (obviamente, na mesma direção em que seu ponto está localizado), esse ponto de dados está fora da elipse.

Segunda abordagem

Essa abordagem é baseada na distribuição suave do kernel.

Estes são alguns dados distribuídos de maneira semelhante aos seus dados:

data1 = RandomVariate[BinormalDistribution[{.3, .7}, {.2, .3}, .8], 500];
data2 = RandomVariate[BinormalDistribution[{.6, .3}, {.4, .15}, .8], 500];
data = Partition[Flatten[Join[{data1, data2}]], 2];

Obtemos uma distribuição suave do kernel nesses valores de dados:

skd = SmoothKernelDistribution[data];

Obtemos um resultado numérico para cada ponto de dados:

eval = Table[{data[[i]], PDF[skd, data[[i]]]}, {i, Length[data]}];

Fixamos um limite e selecionamos todos os dados que são maiores que esse limite:

threshold = 1.2;
dataIn = Select[eval, #1[[2]] > threshold &][[All, 1]];

Aqui temos os dados que ficam fora da região:

dataOut = Complement[data, dataIn];

E agora podemos plotar todos os dados:

Show[ContourPlot[Evaluate@PDF[skd, {x, y}], {x, 0, 1}, {y, 0, 1}, PlotRange -> {{0, 1}, {0, 1}}, PlotPoints -> 50],
ListPlot[dataIn, PlotStyle -> Darker[Green]],
ListPlot[dataOut, PlotStyle -> Red]]

Os pontos de cor verde são aqueles acima do limite e os pontos de cor vermelha são aqueles abaixo do limite.

insira a descrição da imagem aqui

VLC
fonte
Obrigado, sua segunda abordagem me ajuda muito na distribuição do Kernel. Sou programador, não estatístico e sou novato em Mathmatica e R, por isso agradeço muito sua ajuda. Na sua segunda abordagem, fica claro para mim como testar um ponto em que está. Mas como fazer isso na primeira abordagem? Suponho que tenho que comparar meu argumento com a definição de elipsóide. Tou pode fornecer como? Agora eu tenho a esperança de que existem mesmas definições em R, porque eu preciso para usá-lo em RinRuby ...
matejuh
@matejuh Acabei de adicionar mais algumas linhas sobre a primeira abordagem que pode direcioná-lo para uma solução.
VLC
2

A ellipsefunção no ellipsepacote para R gerará essas elipses (na verdade, um polígono que se aproxima da elipse). Você poderia usar essa elipse.

O que pode ser realmente mais fácil é calcular a altura da densidade no seu ponto e ver se ela é maior (dentro da elipse) ou mais baixa (fora da elipse) do que o valor do contorno na elipse. As ellipsefunções internas usam um valor para criar a elipse, você pode começar por lá para encontrar a altura a ser usada.χ2

Greg Snow
fonte
1

Encontrei a resposta em: /programming/2397097/how-can-a-data-ellipse-be-superimposed-on-a-ggplot2-scatterplot

#bootstrap
set.seed(101)
n <- 1000
x <- rnorm(n, mean=2)
y <- 1.5 + 0.4*x + rnorm(n)
df <- data.frame(x=x, y=y, group="A")
x <- rnorm(n, mean=2)
y <- 1.5*x + 0.4 + rnorm(n)
df <- rbind(df, data.frame(x=x, y=y, group="B"))

#calculating ellipses
library(ellipse)
df_ell <- data.frame()
for(g in levels(df$group)){
df_ell <- rbind(df_ell, cbind(as.data.frame(with(df[df$group==g,], ellipse(cor(x, y), 
                                         scale=c(sd(x),sd(y)), 
                                         centre=c(mean(x),mean(y))))),group=g))
}
#drawing
library(ggplot2)
p <- ggplot(data=df, aes(x=x, y=y,colour=group)) + geom_point(size=1.5, alpha=.6) +
  geom_path(data=df_ell, aes(x=x, y=y,colour=group), size=1, linetype=2)

insira a descrição da imagem aqui

Guy L
fonte