Você pode extrair os vetores próprios e os valores via eigen(A)
. No entanto, é mais simples usar a decomposição de Cholesky. Observe que, ao plotar elipses de confiança para dados, os eixos da elipse geralmente são dimensionados para ter comprimento = raiz quadrada dos valores próprios correspondentes, e é isso que a decomposição de Cholesky fornece.
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)
Editar: para plotar os vetores próprios também, você deve usar a abordagem mais complicada. Isso é equivalente à resposta do suncoolsu, apenas usa a notação de matriz para diminuir o código.
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)
Eu acho que esse é o código R que você deseja. Peguei emprestado o código R deste segmento na lista de e-mails. A idéia é basicamente: os semi-diâmetros maior e menor são os dois valores próprios e você gira a elipse pela quantidade de ângulo entre o primeiro vetor eigen e o eixo x
fonte
asp=1
para ter uma proporção de 1 e setas perpendiculares. Alterar seu código paraevs <- sqrt(eigens$values)
dá a mesma elipse da minha resposta.