Wie zeichnet man eine Ellipse aus Eigenwerten und Eigenvektoren in R? [geschlossen]

15

Könnte sich jemand einen R- Code einfallen lassen, um eine Ellipse aus den Eigenwerten und den Eigenvektoren der folgenden Matrix zu zeichnen:

A=(2.20.40.42.8)
MYaseen208
quelle

Antworten:

16

Sie können die Eigenvektoren und -werte über extrahieren eigen(A). Es ist jedoch einfacher, die Cholesky-Zerlegung zu verwenden. Beachten Sie, dass beim Zeichnen von Vertrauensellipsen für Daten die Ellipsenachsen normalerweise so skaliert werden, dass die Länge = Quadratwurzel der entsprechenden Eigenwerte ist, und dies ist, was die Cholesky-Zerlegung ergibt.

ctr    <- c(0, 0)                               # data centroid -> colMeans(dataMatrix)
A      <- matrix(c(2.2, 0.4, 0.4, 2.8), nrow=2) # covariance matrix -> cov(dataMatrix)
RR     <- chol(A)                               # Cholesky decomposition
angles <- seq(0, 2*pi, length.out=200)          # angles for ellipse
ell    <- 1 * cbind(cos(angles), sin(angles)) %*% RR  # ellipse scaled with factor 1
ellCtr <- sweep(ell, 2, ctr, "+")               # center ellipse to the data centroid
plot(ellCtr, type="l", lwd=2, asp=1)            # plot ellipse
points(ctr[1], ctr[2], pch=4, lwd=2)            # plot data centroid

library(car)  # verify with car's ellipse() function
ellipse(c(0, 0), shape=A, radius=0.98, col="red", lty=2)

Bearbeiten: Um auch die Eigenvektoren zu zeichnen, müssen Sie den komplizierteren Ansatz verwenden. Dies entspricht der Antwort von suncoolsu, es wird lediglich die Matrixnotation verwendet, um den Code zu verkürzen.

eigVal  <- eigen(A)$values
eigVec  <- eigen(A)$vectors
eigScl  <- eigVec %*% diag(sqrt(eigVal))  # scale eigenvectors to length = square-root
xMat    <- rbind(ctr[1] + eigScl[1, ], ctr[1] - eigScl[1, ])
yMat    <- rbind(ctr[2] + eigScl[2, ], ctr[2] - eigScl[2, ])
ellBase <- cbind(sqrt(eigVal[1])*cos(angles), sqrt(eigVal[2])*sin(angles)) # normal ellipse
ellRot  <- eigVec %*% t(ellBase)                                          # rotated ellipse
plot((ellRot+ctr)[1, ], (ellRot+ctr)[2, ], asp=1, type="l", lwd=2)
matlines(xMat, yMat, lty=1, lwd=2, col="green")
points(ctr[1], ctr[2], pch=4, col="red", lwd=3)

Bildbeschreibung hier eingeben

caracal
quelle
Würde es Ihnen etwas ausmachen, die Eigenwerte und Eigenvektoren auf dieser Ellipse zu zeichnen? Vielen Dank
MYaseen208
@ MYaseen208 Ich habe meine Antwort so bearbeitet, dass die Eigenvektoren als Ellipsenachsen angezeigt werden. Die halbe Länge der Achsen entspricht der Quadratwurzel der entsprechenden Eigenvektoren.
Karakal
7

Ich denke, das ist der R-Code, den Sie wollen. Ich habe den R-Code von diesem Thread in der R-Mailing-Liste ausgeliehen. Die Grundidee ist: Der Haupt- und der Nebenhalbdurchmesser sind die beiden Eigenwerte, und Sie drehen die Ellipse um den Winkelbetrag zwischen dem ersten Eigenvektor und der x-Achse

mat <- matrix(c(2.2, 0.4, 0.4, 2.8), 2, 2)
eigens <- eigen(mat)
evs <- sqrt(eigens$values)
evecs <- eigens$vectors

a <- evs[1]
b <- evs[2]
x0 <- 0
y0 <- 0
alpha <- atan(evecs[ , 1][2] / evecs[ , 1][1])
theta <- seq(0, 2 * pi, length=(1000))

x <- x0 + a * cos(theta) * cos(alpha) - b * sin(theta) * sin(alpha)
y <- y0 + a * cos(theta) * sin(alpha) + b * sin(theta) * cos(alpha)


png("graph.png")
plot(x, y, type = "l", main = expression("x = a cos " * theta * " + " * x[0] * " and y = b sin " * theta * " + " * y[0]), asp = 1)
arrows(0, 0, a * evecs[ , 1][2], a * evecs[ , 1][2])
arrows(0, 0, b * evecs[ , 2][3], b * evecs[ , 2][2])
dev.off()

Bildbeschreibung hier eingeben

suncoolsu
quelle
Bitte zögern Sie nicht, mich zu korrigieren. Ich denke nicht, dass die Eigen-VECs senkrecht sind (sie müssen theoretisch sein; kann ich etwas falsches zeichnen?).
Suncoolsu
Danke für deine Antwort. Diese Ellipse sieht gut aus, aber dort fehlt etwas. Ich habe den angegebenen Code mit dieser Matrix ausprobiert und es wurde die richtige Ellipse angegeben. Ich frage mich, warum es nicht die richtige Ellipse für verschiedene Varianzen in der Matrix liefert. Jeglicher Kommentar! A=(1551)
MYaseen208
Stellen Sie einfach asp=1ein Seitenverhältnis von 1 und senkrechte Pfeile ein. Ändern Sie Ihren Code, um evs <- sqrt(eigens$values)die gleiche Ellipse wie meine Antwort zu erhalten.
Karakal
3
@ MYaseen208 Deine neue Matrix ist nicht eindeutig positiv: Sie hat negative Eigenwerte und ist keine mögliche Kovarianzmatrix. Ich weiß nicht, welche Ellipse ich in diesem Fall zeichnen soll.
Karakal
@caracal danke! ... ja - ich habe den sqrt Teil verpasst!
Suncoolsu