Warum geben die R-Funktionen 'princomp' und 'prcomp' unterschiedliche Eigenwerte an?

22

Sie können den Zehnkampf-Datensatz {FactoMineR} verwenden, um dies zu reproduzieren. Die Frage ist, warum sich die berechneten Eigenwerte von denen der Kovarianzmatrix unterscheiden.

Hier sind die Eigenwerte mit 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 

Und das gleiche mit 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

Können Sie mir erklären, warum sich die direkt berechneten Eigenwerte von denen unterscheiden? (Die Eigenvektoren sind die gleichen):

> 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

Die alternative prcompMethode liefert auch die gleichen Eigenwerte wie die direkte Berechnung:

> 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

Warum tun PCA/ princompund prcompgibt verschiedene Eigenwerte?

George Dontas
quelle
PCA liefert unterschiedliche Ergebnisse, je nachdem, ob Sie die Kovarianzmatrix oder die Korrelationsmatrix verwenden.
charles.y.zheng
7
Die Unterschiede scheinen relativ klein zu sein, obwohl sie wahrscheinlich zu groß sind, um einfache numerische Probleme zu lösen. Könnte es der Unterschied zwischen der Normalisierung um oder , wenn beispielsweise eine Schätzung der Kovarianz vor der Berechnung der SVD oder der Eigenwertzerlegung berechnet wird? n - 1nn1
Kardinal
7
@ Kardinal Schöne Vermutung! Beachten Sie, dass die beiden unterschiedlichen Folgen von Eigenwerten identische aufeinanderfolgende Verhältnisse haben. Ein Satz ist also ein konstantes Vielfaches des anderen. Das Vielfache ist 1,025 = 41/40 ( genau ). Mir ist unklar, woher das kommt. Vielleicht hat der Datensatz 41 Elemente und das OP enthüllt nur die ersten 10?
whuber
7
@ Cardinal Indeed: Hilfeseite für princomp: "Beachten Sie, dass die Standardberechnung den Divisor N für die Kovarianzmatrix verwendet." prcompHilfeseite für : "Im Gegensatz zu princomp werden Varianzen mit dem üblichen Divisor N-1 berechnet."
Karakal
2
@caracal, Sie sollten Ihren Kommentar in eine Antwort kopieren (und möglicherweise in CW ändern), damit er akzeptiert und die Frage als gelöst markiert werden kann.
Kardinal

Antworten:

16

Wie in den Kommentaren erwähnt, princompwird für den Divisor verwendet, aber für die direkte Berechnung mit beiden wird anstelle von .N - 1 NNprcompcovN1N

Dies wird sowohl im Detailabschnitt von erwähnt help(princomp):

Beachten Sie, dass die Standardberechnung den Divisor 'N' für die Kovarianzmatrix verwendet.

und der Detailbereich von help(prcomp):

Im Gegensatz dazu princompwerden Varianzen mit dem üblichen Divisor N - 1 berechnet.

Sie können dies auch in der Quelle sehen. Das folgende princompQuelltext- Snippet zeigt beispielsweise, dass ( ) als Nenner für die Berechnung verwendet wird .Nn.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
}

Sie können diese Multiplikation vermeiden, indem Sie das covmatArgument anstelle des xArguments angeben.

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

Update bezüglich der PCA-Scores:

Sie können cor = TRUEin Ihrem Aufruf festlegen princomp, dass PCA für die Korrelationsmatrix (anstelle der Kovarianzmatrix) ausgeführt werden soll. Dies führt dazu princomp, dass die Daten bewertet werden, der Nenner wird jedoch weiterhin .NzN

Als Ergebnis princomp(scale(data))$scoresund princomp(data, cor = TRUE)$scoreswird sich durch den Faktor .(N1)/N

Joshua Ulrich
quelle
1
Sie könnten erwägen, "erraten" durch "bereits bestätigt" zu ersetzen (siehe Kommentar-Stream oben). Prost.
Kardinal
@ Kardinal Ich habe diese Kommentare nicht gesehen. Ich habe nur diejenigen gesehen, für die abgestimmt wurde. Vielen Dank. Können Sie auch erklären, warum die Antwort CW lautet? Was sind die Regeln / Richtlinien dafür?
Joshua Ulrich
Kann jemand erraten, warum der Code cv <- cov.wt(z, method="ML")die 2 folgenden Zeilen nicht unnötig macht?
Karakal
2
@Joshua: Mein Vorschlag, die Antwort CW zu machen, beruhte auf der Tatsache, dass die Antwort über einen Strom von Kommentaren erschien und von einer "Community" -Diskussion generiert wurde. Da es in den Kommentaren gelöst wurde, ist es meiner Meinung nach am sinnvollsten, es als Antwort umzuformulieren, die als CW markiert ist, um auf diese Zusammenarbeit hinzuweisen. Dadurch kann die Antwort akzeptiert und die Frage als gelöst markiert werden. (Andernfalls wird es nach einer bestimmten Zeit automatisch von der Software
Kardinal
1
@amoeba es wäre hilfreich gewesen, das in deinem Bearbeitungskommentar zu erwähnen. "860 Zeichen zum Text hinzugefügt" zu einer Antwort mit ~ 450 Zeichen hilft niemandem bei der Beurteilung, ob die Bearbeitung sinnvoll ist.
Joshua Ulrich