¿Por qué las funciones R 'princomp' y 'prcomp' dan diferentes valores propios?

22

Puede usar el conjunto de datos de decatlón {FactoMineR} para reproducir esto. La pregunta es por qué los valores propios calculados difieren de los de la matriz de covarianza.

Aquí están los valores propios usando princomp:

> library(FactoMineR);data(decathlon)
> pr <- princomp(decathlon[1:10], cor=F)
> pr$sd^2
      Comp.1       Comp.2       Comp.3       Comp.4       Comp.5       Comp.6 
1.348073e+02 2.293556e+01 9.747263e+00 1.117215e+00 3.477705e-01 1.326819e-01 
      Comp.7       Comp.8       Comp.9      Comp.10 
6.208630e-02 4.938498e-02 2.504308e-02 4.908785e-03 

Y lo mismo usando PCA:

> res<-PCA(decathlon[1:10], scale.unit=FALSE, ncp=5, graph = FALSE)
> res$eig
          eigenvalue percentage of variance cumulative percentage of variance
comp 1  1.348073e+02           79.659589641                          79.65959
comp 2  2.293556e+01           13.552956464                          93.21255
comp 3  9.747263e+00            5.759799777                          98.97235
comp 4  1.117215e+00            0.660178830                          99.63252
comp 5  3.477705e-01            0.205502637                          99.83803
comp 6  1.326819e-01            0.078403653                          99.91643
comp 7  6.208630e-02            0.036687700                          99.95312
comp 8  4.938498e-02            0.029182305                          99.98230
comp 9  2.504308e-02            0.014798320                          99.99710
comp 10 4.908785e-03            0.002900673                         100.00000

¿Puede explicarme por qué los valores propios calculados directamente difieren de esos? (los vectores propios son los mismos):

> eigen(cov(decathlon[1:10]))$values
 [1] 1.381775e+02 2.350895e+01 9.990945e+00 1.145146e+00 3.564647e-01
 [6] 1.359989e-01 6.363846e-02 5.061961e-02 2.566916e-02 5.031505e-03

Además, el prcompmétodo alternativo proporciona los mismos valores propios que el cálculo directo:

> prc <- prcomp(decathlon[1:10])
> prc$sd^2
 [1] 1.381775e+02 2.350895e+01 9.990945e+00 1.145146e+00 3.564647e-01
 [6] 1.359989e-01 6.363846e-02 5.061961e-02 2.566916e-02 5.031505e-03

¿Por qué PCA/ princompy prcompdar diferentes valores propios?

George Dontas
fuente
PCA le dará diferentes resultados dependiendo de si usa la matriz de covarianza o la matriz de correlación.
charles.y.zheng
77
Las diferencias parecen relativamente pequeñas, aunque probablemente demasiado grandes para ser simples problemas numéricos. ¿Podría ser la diferencia entre normalizar por o n - 1 , por ejemplo, al calcular una estimación de la covarianza antes de calcular la SVD o la descomposición de los valores propios? nortenorte-1
Cardenal
77
@cardinal ¡Buena conjetura! Observe que las dos secuencias diferentes de valores propios tienen proporciones sucesivas idénticas. Por lo tanto, un conjunto es un múltiplo constante del otro. El múltiplo es 1.025 = 41/40 ( exactamente ). No me queda claro de dónde viene esto. ¿Quizás el conjunto de datos tiene 41 elementos y el OP revela solo los primeros 10?
whuber
77
@cardinal Indeed: página de ayuda para princomp: "Tenga en cuenta que el cálculo predeterminado utiliza el divisor N para la matriz de covarianza". Página de ayuda para prcomp: "A diferencia de princomp, las variaciones se calculan con el divisor habitual N-1".
caracal
2
@caracal, debe copiar su comentario en una respuesta (y tal vez hacerlo CW) para que pueda aceptarse y la pregunta pueda marcarse como resuelta.
cardenal

Respuestas:

16

princompnorteprcompcovnorte-1norte

Esto se menciona en la sección Detalles de help(princomp):

Tenga en cuenta que el cálculo predeterminado utiliza el divisor 'N' para la matriz de covarianza.

y la sección Detalles de help(prcomp):

A diferencia princomp, las variaciones se calculan con el divisor habitual N - 1.

princompnorte ( n.obs) se usa como denominador al calcular cv.

else if (is.null(covmat)) {
    dn <- dim(z)
    if (dn[1L] < dn[2L]) 
        stop("'princomp' can only be used with more units than variables")
    covmat <- cov.wt(z)
    n.obs <- covmat$n.obs
    cv <- covmat$cov * (1 - 1/n.obs)
    cen <- covmat$center
}

Puede evitar esta multiplicación especificando el covmat argumento en lugar del xargumento.

princomp(covmat = cov(iris[,1:4]))$sd^2

Actualización sobre las puntuaciones de PCA:

cor = TRUEprincompprincompznorte

princomp(scale(data))$scoresprincomp(data, cor = TRUE)$scores(norte-1)/ /norte

Joshua Ulrich
fuente
1
Puede considerar reemplazar "adivinado" por "ya confirmado" (vea la secuencia de comentarios anterior). También puede considerar editar su respuesta para que sea CW. Aclamaciones.
cardenal
@cardinal No vi esos comentarios. Solo vi aquellos que habían sido votados. Gracias. Además, ¿podría explicar la razón detrás de hacer la respuesta CW? ¿Cuáles son las reglas / pautas para eso?
Joshua Ulrich
¿Alguien puede adivinar por qué el código no está simplemente cv <- cov.wt(z, method="ML")haciendo innecesarias las 2 líneas siguientes?
caracal
2
@Joshua: Mi sugerencia con respecto a hacer la respuesta CW se debió al hecho de que la respuesta apareció a través de un flujo de comentarios y fue generada por una discusión "comunitaria". Dado que se resolvió en los comentarios, creo que tiene más sentido reformularlo como una respuesta, marcada como CW para indicar esta colaboración, y esto permite que la respuesta sea aceptada y que la pregunta se marque como resuelta. (De lo contrario, el software lo reactivará automáticamente después de un cierto período de tiempo).
Cardenal
2
@amoeba hubiera sido útil mencionar eso en tu comentario de edición. "Se agregaron 860 caracteres al cuerpo" a una respuesta de ~ 450 caracteres no ayuda a nadie a evaluar si la edición es razonable.
Joshua Ulrich