Dlaczego funkcje R „princomp” i „prcomp” dają różne wartości własne?

22

Aby to odtworzyć, możesz użyć zestawu danych decathlon {FactoMineR}. Pytanie brzmi, dlaczego obliczone wartości własne różnią się od wartości macierzy kowariancji.

Oto wartości własne przy użyciu 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 

I to samo za pomocą 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

Czy możesz mi wyjaśnić, dlaczego bezpośrednio obliczone wartości własne różnią się od nich? (wektory własne są takie same):

> 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

Ponadto prcompmetoda alternatywna daje te same wartości własne, co bezpośrednie obliczenia:

> 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

Dlaczego PCA/ princompi prcomppodać różne wartości własne?

George Dontas
źródło
PCA da różne wyniki w zależności od tego, czy używasz macierzy kowariancji, czy macierzy korelacji.
charles.y.zheng
7
nn-1
7
@cardinal Zgadnij! Zauważ, że dwie różne sekwencje wartości własnych mają identyczne kolejne stosunki. Zatem jeden zestaw jest stałą wielokrotnością drugiego. Wielokrotność wynosi 1,025 = 41/40 ( dokładnie ). Nie jest dla mnie jasne, skąd to pochodzi. Może zestaw danych zawiera 41 elementów, a PO ujawnia tylko pierwsze 10?
whuber
7
@ cardinal Rzeczywiście: Strona pomocy dla princomp: „Zauważ, że domyślne obliczenia używają dzielnika N dla macierzy kowariancji”. Strona pomocy dla prcomp: „W przeciwieństwie do princomp, wariancje są obliczane za pomocą zwykłego dzielnika N-1”.
karakal
2
@caracal, powinieneś skopiować swój komentarz do odpowiedzi (i być może uczynić go CW), aby mógł zostać zaakceptowany, a pytanie oznaczone jako rozwiązane.
kardynał

Odpowiedzi:

16

princompNprcompcovN1N

Jest to wspomniane zarówno w sekcji Szczegółyhelp(princomp) :

Zauważ, że domyślne obliczenia używają dzielnika „N” dla macierzy kowariancji.

oraz sekcja Szczegółyhelp(prcomp) :

W przeciwieństwie do princompwariancji oblicza się zwykły dzielnik N - 1.

princompNn.obscv

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
}

Możesz uniknąć tego mnożenia, podając covmatargument zamiast xargumentu.

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

Aktualizacja dotycząca wyników PCA:

Możesz ustawić cor = TRUEw swoim wywołaniu princomp, aby wykonać PCA na macierzy korelacji (zamiast macierzy kowariancji). Spowoduje to, princompabyzN

W rezultacie princomp(scale(data))$scoresprincomp(data, cor = TRUE)$scores(N1)/N. .

Joshua Ulrich
źródło
1
Możesz rozważyć zamianę „zgadł” na „już potwierdzony” (patrz strumień komentarzy powyżej.) Możesz również rozważyć edycję swojej odpowiedzi, aby była CW. Twoje zdrowie.
kardynał
@cardinal Nie widziałem tych komentarzy. Widziałem tylko te, które zostały poddane pod głosowanie. Dzięki. Czy mógłbyś również wyjaśnić uzasadnienie udzielenia odpowiedzi CW? Jakie są zasady / wytyczne?
Joshua Ulrich
Czy ktoś może zgadnąć, dlaczego kod nie cv <- cov.wt(z, method="ML")powoduje po prostu, że 2 następujące linie nie są potrzebne?
caracal
2
@Joshua: Moja sugestia dotycząca udzielenia odpowiedzi CW wynikała z faktu, że odpowiedź pojawiła się w strumieniu komentarzy i została wygenerowana w wyniku dyskusji „społeczności”. Ponieważ zostało to rozwiązane w komentarzach, myślę, że najbardziej sensowne jest przeformułowanie go jako odpowiedzi, oznaczonej jako CW, aby wskazać tę współpracę, a to pozwala na zaakceptowanie odpowiedzi i pytanie jako oznaczone jako rozwiązane. (W przeciwnym razie oprogramowanie zostanie automatycznie przywrócone po pewnym czasie).
kard.
2
@amoeba warto wspomnieć o tym w komentarzu do edycji. „Dodano 860 znaków do treści” w odpowiedzi na ~ 450 znaków nie pomaga nikomu ocenić, czy edycja jest rozsądna.
Joshua Ulrich