Wie erhalte ich eine Ellipsenregion aus bivariaten normalverteilten Daten?

11

Ich habe Daten, die aussehen wie:

Zahl

Ich habe versucht, eine Normalverteilung anzuwenden (die Schätzung der Kerneldichte funktioniert besser, aber ich brauche keine so große Präzision), und sie funktioniert recht gut. Dichtediagramm macht eine Ellipse.

Ich brauche diese Ellipsenfunktion, um zu entscheiden, ob ein Punkt innerhalb der Ellipsenregion liegt oder nicht. Wie geht das?

R- oder Mathematica-Code sind willkommen.

matejuh
quelle

Antworten:

17

Corsario bietet eine gute Lösung in einem Kommentar: Verwenden Sie die Kerneldichtefunktion, um die Aufnahme in einen Levelsatz zu testen.

Eine andere Interpretation der Frage besteht darin, dass ein Verfahren zum Testen der Aufnahme in die Ellipsen angefordert wird, die durch eine bivariate normale Annäherung an die Daten erzeugt werden. Lassen Sie uns zunächst einige Daten generieren, die der Abbildung in der Frage entsprechen:

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

Die Ellipsen werden durch den ersten und zweiten Moment der Daten bestimmt:

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

Die Formel erfordert die Inversion der Varianz-Kovarianz-Matrix:

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

Die Ellipsenfunktion "Höhe" ist das Negative des Logarithmus der bivariaten Normaldichte :

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

(Ich habe eine additive Konstante gleich log ( 2 π √) ignoriertLog((2πdet((Σ)) .)

Um dies zu testen , zeichnen wir einige seiner Konturen. Dazu muss ein Punktgitter in x- und y-Richtung erstellt werden:

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

Berechnen Sie die Höhenfunktion in diesem Raster und zeichnen Sie sie auf:

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)

Konturdiagramm

Offensichtlich funktioniert es. Daher der Test, um festzustellen, ob ein Punkt innerhalb einer elliptischen Kontur auf Ebene c liegt((s,t)c liegt,

ellipse(s,t) <= c

Mathematica erledigt den Job auf die gleiche Weise: Berechnen Sie die Varianz-Kovarianz-Matrix der Daten, invertieren Sie diese, konstruieren Sie die ellipseFunktion, und schon sind Sie fertig.

whuber
quelle
Vielen Dank an alle, besonders an @whuber. Genau das brauche ich.
Matejuh
Übrigens. Gibt es eine einfache Lösung für Konturen zur Schätzung der Kerneldichte? Denn wenn ich strenger sein möchte, sehen meine Daten folgendermaßen aus: github.com/matejuh/doschecker_wiki_images/raw/master/… resp. github.com/matejuh/doschecker_wiki_images/raw/master/…
matejuh
Ich kann in R keine einfache Lösung finden. Erwägen Sie die Verwendung der Funktion "SmoothKernelDistribution" von Mathematica 8.
whuber
2
Entspricht das Niveau dem Konfidenzniveau? Das glaube ich nicht. Wie kann ich das bitte machen?
Matejuh
Dies erfordert eine neue Frage, da Sie angeben müssen, worauf Sie vertrauen möchten, und - anhand Ihrer Diagramme zu urteilen - Bedenken bestehen, ob solche Ellipsen überhaupt eine angemessene Beschreibung der Daten darstellen.
whuber
9

Die Darstellung ist mit der ellipse()Funktion des mixtoolsPakets für R einfach :

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

Geben Sie hier die Bildbeschreibung ein

Stéphane Laurent
quelle
5

Erste Ansatz

Sie können diesen Ansatz in Mathematica ausprobieren.

Lassen Sie uns einige bivariate Daten generieren:

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

Dann müssen wir dieses Paket laden:

Needs["MultivariateStatistics`"]

Und nun:

ellPar=EllipsoidQuantile[data, {0.9}]

gibt eine Ausgabe aus, die eine 90% -Konfidenzellipse definiert. Die Werte, die Sie von dieser Ausgabe erhalten, haben das folgende Format:

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

x1 und x2 geben den Punkt an, an dem die Ellipse zentriert ist, r1 und r2 die Halbachsenradien und d1, d2, d3 und d4 die Ausrichtungsrichtung.

Sie können dies auch zeichnen:

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

Die allgemeine parametrische Form der Ellipse lautet:

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

Und Sie können es so zeichnen:

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

Sie können eine Überprüfung basierend auf reinen geometrischen Informationen durchführen: Wenn der euklidische Abstand zwischen dem Mittelpunkt der Ellipse (ellPar [[1,1]]) und Ihrem Datenpunkt größer ist als der Abstand zwischen dem Mittelpunkt der Ellipse und dem Rand von die Ellipse (offensichtlich in der gleichen Richtung, in der sich Ihr Punkt befindet), dann befindet sich dieser Datenpunkt außerhalb der Ellipse.

Zweiter Ansatz

Dieser Ansatz basiert auf der reibungslosen Kernelverteilung.

Dies sind einige Daten, die auf ähnliche Weise wie Ihre Daten verteilt werden:

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

Wir erhalten eine reibungslose Kernelverteilung für diese Datenwerte:

skd = SmoothKernelDistribution[data];

Wir erhalten für jeden Datenpunkt ein numerisches Ergebnis:

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

Wir legen einen Schwellenwert fest und wählen alle Daten aus, die höher als dieser Schwellenwert sind:

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

Hier erhalten wir die Daten, die außerhalb der Region liegen:

dataOut = Complement[data, dataIn];

Und jetzt können wir alle Daten zeichnen:

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

Die grün gefärbten Punkte befinden sich oberhalb des Schwellenwerts und die rot gefärbten Punkte befinden sich unterhalb des Schwellenwerts.

Geben Sie hier die Bildbeschreibung ein

VLC
quelle
Vielen Dank, Ihr zweiter Ansatz hilft mir sehr bei der Kernel-Verteilung. Ich bin Programmierer, nicht statistisch und ich bin Neuling in Mathmatica und R, daher schätze ich Ihre Hilfe sehr. In Ihrem zweiten Ansatz ist mir klar, wie ich einen Punkt testen kann, an dem er liegt. Aber wie geht das im ersten Ansatz? Ich nehme an, dass ich meinen Punkt mit der Ellipsoiddefinition vergleichen muss. Können Sie bitte angeben, wie? Jetzt muss ich hoffen, dass es die gleichen Definitionen in R gibt, weil ich es in RinRuby verwenden muss ...
matejuh
@matejuh Ich habe gerade ein paar Zeilen mehr über den ersten Ansatz hinzugefügt, der Sie zu einer Lösung führen könnte.
VLC
2

Die ellipseFunktion in derellipse Paket für R generiert diese Ellipsen (tatsächlich ein Polygon, das sich der Ellipse annähert). Sie könnten diese Ellipse verwenden.

ellipseχ2

Greg Snow
quelle
1

Ich fand die Antwort unter: /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)

Geben Sie hier die Bildbeschreibung ein

Guy L.
quelle