Comment obtenir une région d'ellipse à partir de données distribuées normales bivariées?

11

J'ai des données qui ressemblent à:

Figure

J'ai essayé d'appliquer une distribution normale (l'estimation de la densité du noyau fonctionne mieux, mais je n'ai pas besoin d'une telle précision) et cela fonctionne assez bien. Le tracé de densité fait une ellipse.

J'ai besoin d'obtenir cette fonction d'ellipse pour décider si un point se trouve dans la région de l'ellipse ou non. Comment faire ça?

Les codes R ou Mathematica sont les bienvenus.

matejuh
la source

Réponses:

18

Corsario fournit une bonne solution dans un commentaire: utilisez la fonction de densité du noyau pour tester l'inclusion dans un ensemble de niveaux.

Une autre interprétation de la question est qu'elle demande une procédure pour tester l'inclusion dans les ellipses créées par une approximation normale bivariée des données. Pour commencer, générons des données qui ressemblent à l'illustration de la question:

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

Les ellipses sont déterminées par les premier et deuxième moments des données:

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

La formule nécessite l'inversion de la matrice variance-covariance:

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

La fonction "hauteur" de l'ellipse est le négatif du logarithme de la densité normale bivariée :

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

(J'ai ignoré une constante additive égale à .)log(2πdet(Σ))

Pour tester cela , dessinons certains de ses contours. Cela nécessite de générer une grille de points dans les directions x et y:

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

Calculez la fonction de hauteur sur cette grille et tracez-la:

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)

Tracé de contour

Evidemment ça marche. Par conséquent, le test pour déterminer si un point se trouve à l'intérieur d'un contour elliptique au niveau est(s,t)c

ellipse(s,t) <= c

Mathematica fait le travail de la même manière: calculez la matrice de variance-covariance des données, inversez-la, construisez la ellipsefonction, et vous êtes tous ensemble.

whuber
la source
Merci à tous, en particulier @whuber. C'est exactement ce dont j'ai besoin.
matejuh
Btw. existe-t-il une solution simple pour les contours d'estimation de la densité du noyau? Parce que si je veux être plus strict, mes données ressemblent à: github.com/matejuh/doschecker_wiki_images/raw/master/… resp. github.com/matejuh/doschecker_wiki_images/raw/master/…
matejuh
Je ne trouve pas de solution simple dans R. Envisagez d'utiliser la fonction "SmoothKernelDistribution" de Mathematica 8.
whuber
2
Les niveaux correspondent-ils au niveau de confiance? Je ne pense pas. Comment puis-je faire ça s'il vous plait?
matejuh
Cela nécessite une nouvelle question, car vous devez spécifier ce que vous recherchez et, à en juger par vos graphiques, vous vous demandez si ces ellipses sont des descriptions adéquates des données en premier lieu.
whuber
9

L'intrigue est simple avec la ellipse()fonction du mixtoolspackage pour 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") 

entrez la description de l'image ici

Stéphane Laurent
la source
5

Première approche

Vous pouvez essayer cette approche dans Mathematica.

Générons quelques données bivariées:

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

Ensuite, nous devons charger ce package:

Needs["MultivariateStatistics`"]

Et maintenant:

ellPar=EllipsoidQuantile[data, {0.9}]

donne une sortie qui définit une ellipse de confiance à 90%. Les valeurs que vous obtenez à partir de cette sortie sont au format suivant:

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

x1 et x2 spécifient le point auquel l'ellipse centrée, r1 et r2 spécifient les rayons semi-axes, et d1, d2, d3 et d4 spécifient la direction d'alignement.

Vous pouvez également tracer ceci:

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

La forme paramétrique générale de l'ellipse est:

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

Et vous pouvez le tracer de cette façon:

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

Vous pouvez effectuer une vérification basée sur des informations géométriques pures: si la distance euclidienne entre le centre de l'ellipse (ellPar [[1,1]]) et votre point de données est supérieure à la distance entre le centre de l'ellipse et la frontière de l'ellipse (évidemment, dans la même direction où se trouve votre point), alors ce point de données est en dehors de l'ellipse.

Deuxième approche

Cette approche est basée sur la distribution fluide du noyau.

Voici quelques données distribuées de manière similaire à vos données:

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

Nous obtenons une distribution fluide du noyau sur ces valeurs de données:

skd = SmoothKernelDistribution[data];

Nous obtenons un résultat numérique pour chaque point de données:

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

Nous fixons un seuil et nous sélectionnons toutes les données supérieures à ce seuil:

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

Ici, nous obtenons les données qui se trouvent en dehors de la région:

dataOut = Complement[data, dataIn];

Et maintenant, nous pouvons tracer toutes les données:

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

Les points de couleur verte sont ceux au-dessus du seuil et les points de couleur rouge sont ceux au-dessous du seuil.

entrez la description de l'image ici

VLC
la source
Merci, votre deuxième approche m'aide beaucoup avec la distribution du noyau. Je suis programmeur, pas statistique et je suis novice en Mathmatica et R, donc j'apprécie beaucoup votre aide. Dans votre deuxième approche, il est clair pour moi comment tester un point où il se trouve. Mais comment faire cela en première approche? Je suppose que je dois comparer mon point avec la définition des ellipsoïdes. Pouvez-vous s'il vous plaît fournir comment? Maintenant, je dois espérer qu'il y a les mêmes définitions dans R, car je dois l'utiliser dans RinRuby ...
matejuh
@matejuh Je viens d'ajouter quelques lignes supplémentaires sur la première approche qui pourraient vous orienter vers une solution.
VLC
2

La ellipsefonction dans le ellipsepackage pour R générera ces ellipses (en fait un polygone se rapprochant de l'ellipse). Vous pouvez utiliser cette ellipse.

Ce qui pourrait en fait être plus facile, c'est de calculer la hauteur de la densité à votre point et de voir si elle est plus élevée (à l'intérieur de l'ellipse) ou plus basse (à l'extérieur de l'ellipse) que la valeur du contour à l'ellipse. Les ellipsefonctions internes utilisent une valeur pour créer l'ellipse, vous pouvez commencer par là pour trouver la hauteur à utiliser.χ2

Greg Snow
la source
1

J'ai trouvé la réponse sur: /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)

entrez la description de l'image ici

Guy L
la source