Como plotar uma elipse a partir de valores próprios e vetores próprios em R? [fechadas]

15

Alguém poderia inventar o código R para plotar uma elipse a partir dos valores próprios e dos vetores próprios da seguinte matriz

A=(2.20.40.42.8)
MYaseen208
fonte

Respostas:

16

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)

insira a descrição da imagem aqui

caracal
fonte
Você se importaria de plotar os autovalores e autovetores nessa elipse? Graças
MYaseen208
@ MYaseen208 Editei minha resposta para mostrar os autovetores como os eixos da elipse. Metade do comprimento dos eixos é igual à raiz quadrada dos vetores próprios correspondentes.
Caracal
7

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

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

insira a descrição da imagem aqui

suncoolsu
fonte
sinta-se à vontade para me corrigir. Eu não acho que os eigen vecs sejam perpendiculares (eles devem estar na teoria; posso estar tramando algo errado?).
suncoolsu
Obrigado pela sua resposta. Essa elipse parece boa, mas algo está faltando lá. Eu tentei o seu código fornecido com esta matriz UMA=(1-5-51)e deu elipse correta. Gostaria de saber por que ele não está fornecendo a elipse correta para diferentes variações na matriz. Qualquer comentário!
precisa saber é o seguinte
Basta definir asp=1para ter uma proporção de 1 e setas perpendiculares. Alterar seu código para evs <- sqrt(eigens$values)dá a mesma elipse da minha resposta.
Caracal
3
@ MYaseen208 Sua nova matriz não é positiva definida: possui valores próprios negativos e não é uma matriz de covariância possível. Não sei qual elipse desenhar nesse caso.
Caracal
@caracal thanks! ... sim - eu perdi a parte sqrt!
suncoolsu