¿Cómo obtener la región de elipse a partir de datos distribuidos normales bivariados?

Tengo datos que se parecen a:

Traté de aplicar una distribución normal (la estimación de la densidad del núcleo funciona mejor, pero no necesito tanta precisión) y funciona bastante bien. La gráfica de densidad hace una elipse.

Necesito obtener esa función de elipse para decidir si un punto se encuentra dentro de la región de la elipse o no. ¿Como hacer eso?

R o código de Mathematica son bienvenidos.

Corsario proporciona una buena solución en un comentario: use la función de densidad del núcleo para probar la inclusión dentro de un conjunto de niveles.

Otra interpretación de la pregunta es que solicita un procedimiento para evaluar la inclusión dentro de las elipses creadas por una aproximación normal bivariada a los datos. Para comenzar, generemos algunos datos que se parecen 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 de "altura" de 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. Eso requiere 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 tracela:

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 encuentra dentro de un contorno elíptico en el nivel es$(s,t)$$c$

ellipse(s,t) <= c

Mathematica hace el trabajo de la misma manera: calcula la matriz de varianza-covarianza de los datos, invierte eso, construye la ellipsefunción y estás listo.

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

Primer enfoque

Puede probar este enfoque en Mathematica.

Generemos algunos datos bivariados:

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

Entonces necesitamos 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 obtiene 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 la elipse en el centro, 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 puedes trazarlo 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}}]

Puede realizar una verificación basada en información geométrica pura: 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 elipse (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 fluida del núcleo.

Estos son algunos datos distribuidos de manera similar a sus datos:

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 de kernel suave en 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]}];

Arreglamos un umbral y seleccionamos todos los datos que son más altos que este umbral:

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

Aquí obtenemos los datos que quedan 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.

La ellipsefunción en el ellipsepaquete para R generará estas elipses (en realidad un polígono que se aproxima a la elipse). Podrías 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 más alta (dentro de la elipse) o más baja (fuera de la elipse) que el valor del contorno en la elipse. Las ellipsefunciones internas usan un valor para crear la elipse, puede comenzar allí para encontrar la altura a usar.$\chi^2$

Encontré la respuesta en: /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)