¿Cómo obtener la región de la elipse a partir de datos con distribución normal bivariante?

Question

Tengo datos que parecen:

Intenté aplicar la distribución normal (la estimación de la densidad del núcleo funciona mejor, pero no necesito una precisión tan grande) en ella y funciona bastante bien. El gráfico de densidad hace una elipse.

Necesito que esa función de la elipse decida si un punto se encuentra dentro de la región de la elipse o no. ¿Cómo hacerlo?

Se aceptan códigos de R o Mathematica.

Answer 1

Corsario ofrece una buena solución en un comentario: utilizar la función de densidad del núcleo para comprobar la inclusión en un conjunto de niveles.

Otra interpretación de la pregunta es que solicita un procedimiento para comprobar la inclusión dentro de las elipses creadas por un normal bivariante aproximación a los datos. Para empezar, vamos a generar unos datos que se parezcan a la ilustración de la pregunta:

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

Las elipses están determinadas por el primer y segundo momento de los datos:

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

La fórmula requiere la inversión de la matriz de varianza-covarianza:

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

La función "altura" de la elipse es el negativo del logaritmo de la densidad normal bivariada :

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

(He ignorado una constante aditiva igual a $\log(2\pi\sqrt{\det(\Sigma)})$ .)

Para probar esto dibujemos algunos de sus contornos. Para ello es necesario generar una cuadrícula de puntos en las direcciones x e y:

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

Calcule la función de altura en esta cuadrícula y represéntela:

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)

Evidentemente, funciona. Por lo tanto, la prueba para determinar si un punto $(s,t)$ se encuentra dentro de un contorno elíptico en el nivel $c$ es

ellipse(s,t) <= c

Mathematica hace el trabajo de la misma manera: calcula la matriz de varianza-covarianza de los datos, la invierte, construye la ellipse y ya está todo listo.

Answer 2

La trama es sencilla con el ellipse() función del mixtools 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") 

Answer 3

Primera aproximación

Puedes probar este enfoque en Mathematica.

Vamos a generar algunos datos bivariados:

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

Entonces tenemos que cargar este paquete:

Needs["MultivariateStatistics`"]

Y, ahora:

ellPar=EllipsoidQuantile[data, {0.9}]

da una salida que define una elipse de confianza del 90%. Los valores que se obtienen de esta salida tienen el siguiente formato:

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

x1 y x2 especifican el punto en el que se centra la elipse, r1 y r2 especifican los radios del semieje, y d1, d2, d3 y d4 especifican la dirección de alineación.

También puedes trazar esto:

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

La forma paramétrica general de la elipse es:

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]}

Y se puede trazar de esta manera:

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}}]

Podría realizar una comprobación basada en información puramente geométrica: si la distancia euclidiana entre el centro de la elipse (ellPar[[1,1]]) y su punto de datos es mayor que la distancia entre el centro de la elipse y el borde de la misma (obviamente, en la misma dirección en la que se encuentra su punto), entonces ese punto de datos está fuera de la elipse.

Segundo enfoque

Este enfoque se basa en la distribución kernel suave.

Estos son algunos datos distribuidos de forma similar a los suyos:

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];

Obtenemos una distribución kernel suave sobre estos valores de datos:

skd = SmoothKernelDistribution[data];

Obtenemos un resultado numérico para cada punto de datos:

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

Fijamos un umbral y seleccionamos todos los datos que sean superiores a este umbral:

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

Aquí obtenemos los datos que caen fuera de la región:

dataOut = Complement[data, dataIn];

Y ahora podemos trazar todos los datos:

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]]

Los puntos de color verde son los que están por encima del umbral y los puntos de color rojo son los que están por debajo del umbral.

Answer 4

El ellipse en la función ellipse para R generará estas elipses (en realidad un polígono que se aproxima a la elipse). Puedes usar esa elipse.

Lo que en realidad podría ser más fácil es calcular la altura de la densidad en su punto y ver si es mayor (dentro de la elipse) o menor (fuera de la elipse) que el valor del contorno en la elipse. El ellipse utilizan una función $\chi^2$ para crear la elipse, podrías empezar por ahí para encontrar la altura a utilizar.

Answer 5

Encontré la respuesta en: https://stackoverflow.com/questions/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)