Sollte eine SVD für eine Zufallsmatrix überhaupt nichts erklären? Was mache ich falsch?

Wenn ich eine 2-D-Matrix konstruiere, die vollständig aus Zufallsdaten besteht, würde ich erwarten, dass die PCA- und SVD-Komponenten im Wesentlichen nichts erklären.

Stattdessen scheint die erste SVD-Spalte 75% der Daten zu erklären. Wie kann das möglich sein? Was mache ich falsch?

Hier ist die Handlung:

Bildbeschreibung hier eingeben

Hier ist der R-Code:

set.seed(1)
rm(list=ls())
m <- matrix(runif(10000,min=0,max=25), nrow=100,ncol=100)
svd1 <- svd(m, LINPACK=T)
par(mfrow=c(1,4))
image(t(m)[,nrow(m):1])
plot(svd1$d,cex.lab=2, xlab="SVD Column",ylab="Singluar Value",pch=19)

percentVarianceExplained = svd1$d^2/sum(svd1$d^2) * 100
plot(percentVarianceExplained,ylim=c(0,100),cex.lab=2, xlab="SVD Column",ylab="Percent of variance explained",pch=19)

cumulativeVarianceExplained = cumsum(svd1$d^2/sum(svd1$d^2)) * 100
plot(cumulativeVarianceExplained,ylim=c(0,100),cex.lab=2, xlab="SVD column",ylab="Cumulative percent of variance explained",pch=19)

Aktualisieren

Vielen Dank an Aaron. Wie Sie bemerkt haben, bestand die Korrektur darin, die Matrix so zu skalieren, dass die Zahlen um 0 zentriert sind (dh der Mittelwert ist 0).

m <- scale(m, scale=FALSE)

Hier ist das korrigierte Bild, das für eine Matrix mit Zufallsdaten zeigt, dass die erste SVD-Spalte erwartungsgemäß nahe bei 0 liegt.

Bild korrigiert

r pca svd Contango
quelle

[0, 1]^{100}

$[0,1]^{100}$

R^{100}

$\mathbb{R}^{100}$

R^{n}

$\mathbb{R}^n$

n

$n$

1 / 3

$1/3$

n / 3 - (n - 1) / 12

$n/3-(n-1)/12$

1 / 12

$1/12$

(n / 3 - (n - 1) / 12) / (n / 3) = 3 / 4 + 1 / (4 n)

$(n/3-(n-1)/12)/(n/3)=3/4+1/(4n)$

n = 100

$n=100$

75.25

$75.25$

Antworten:

Der erste PC erklärt, dass die Variablen nicht um Null zentriert sind. Wenn Sie zuerst skalieren oder Ihre Zufallsvariablen um Null zentrieren, erhalten Sie das erwartete Ergebnis. Zum Beispiel:

m <- matrix(runif(10000,min=0,max=25), nrow=100,ncol=100)
m <- scale(m, scale=FALSE)

m <- matrix(runif(10000,min=-25,max=25), nrow=100,ncol=100)

Aaron hat Stack Overflow verlassen
quelle

Sie sprechen einen guten Punkt an, aber ich denke, das erzählt nur einen Teil der Geschichte. In der Tat würde ich vermuten, dass das OP jede dieser Möglichkeiten ausprobiert und immer noch vom Ergebnis überrascht ist. Tatsache ist, dass die singulären Werte, weil sie in der Ausgabe geordnet sind, nicht gleichmäßig verteilt erscheinen (und tatsächlich nicht), wie es naiv von "zufälligen" Daten erwartet werden könnte. Die Marchenko-Pastur-Distribution regelt in diesem Fall ihr Verhalten.

Kardinal

@ Aaron Danke, du hattest absolut recht. Ich habe oben ein Diagramm der korrigierten Ausgabe hinzugefügt, um zu zeigen, wie schön das Ergebnis ist.

Contango

@cardinal Vielen Dank für Ihren Kommentar, Sie haben absolut Recht (siehe das Diagramm, das durch den korrigierten Code oben erstellt wurde). Ich glaube, die SVD-Werte würden mit abnehmender Matrix weniger gleichmäßig verteilt, da eine kleinere Matrix mit größerer Wahrscheinlichkeit Muster haben würde, die nicht durch das Gesetz der großen Zahlen gequetscht würden.

Contango

Ich werde Ihrer Frage eine visuellere Antwort hinzufügen, indem ich einen Nullmodellvergleich verwende. Die Prozedur mischt die Daten in jeder Spalte nach dem Zufallsprinzip, um die Gesamtvarianz beizubehalten, während die Kovarianz zwischen Variablen (Spalten) verloren geht. Dies wird mehrmals durchgeführt und die resultierende Verteilung der Singularwerte in der randomisierten Matrix wird mit den Originalwerten verglichen.

Ich benutze prcompstatt svdfür die Matrixzerlegung, aber die Ergebnisse sind ähnlich:

set.seed(1)
m <- matrix(runif(10000,min=0,max=25), nrow=100,ncol=100)

S <- svd(scale(m, center = TRUE, scale=FALSE))
P <- prcomp(m, center = TRUE, scale=FALSE)
plot(S$d, P$sdev) # linearly related

Der Nullmodellvergleich wird in der folgenden zentrierten Matrix durchgeführt:

library(sinkr) # https://github.com/marchtaylor/sinkr

# centred data
Pnull <- prcompNull(m, center = TRUE, scale=FALSE, nperm = 100)
Pnull$n.sig
boxplot(Pnull$Lambda[,1:20], ylim=range(Pnull$Lambda[,1:20], Pnull$Lambda.orig[1:20]), outline=FALSE, col=8, border="grey50", log="y", main=paste("m (center=FALSE); n sig. =", Pnull$n.sig))
lines(apply(Pnull$Lambda, 2, FUN=quantile, probs=0.95))
points(Pnull$Lambda.orig[1:20], pch=16)

Das Folgende ist ein Boxplot der permutierten Matrix, wobei das 95% -Quantil jedes Singularwerts als durchgezogene Linie dargestellt ist. Die ursprünglichen PCA-Werte msind die Punkte. alle liegen unterhalb der 95% -Linie - daher ist ihre Amplitude nicht von zufälligem Rauschen zu unterscheiden.

Dieselbe Prozedur kann mit der nicht zentrierten Version von mmit demselben Ergebnis durchgeführt werden - Keine signifikanten singulären Werte:

# centred data
Pnull <- prcompNull(m, center = FALSE, scale=FALSE, nperm = 100)
Pnull$n.sig
boxplot(Pnull$Lambda[,1:20], ylim=range(Pnull$Lambda[,1:20], Pnull$Lambda.orig[1:20]), outline=FALSE, col=8, border="grey50", log="y", main=paste("m (center=TRUE); n sig. =", Pnull$n.sig))
lines(apply(Pnull$Lambda, 2, FUN=quantile, probs=0.95))
points(Pnull$Lambda.orig[1:20], pch=16)

Schauen wir uns zum Vergleich einen Datensatz mit einem nicht zufälligen Datensatz an: iris

# iris dataset example
m <- iris[,1:4]
Pnull <- prcompNull(m, center = TRUE, scale=FALSE, nperm = 100)
Pnull$n.sig
boxplot(Pnull$Lambda, ylim=range(Pnull$Lambda, Pnull$Lambda.orig), outline=FALSE, col=8, border="grey50", log="y", main=paste("m (center=FALSE); n sig. =", Pnull$n.sig))
lines(apply(Pnull$Lambda, 2, FUN=quantile, probs=0.95))
points(Pnull$Lambda.orig[1:20], pch=16)

Hier ist der 1. Singularwert signifikant und erklärt über 92% der Gesamtvarianz:

P <- prcomp(m, center = TRUE)
P$sdev^2 / sum(P$sdev^2)
# [1] 0.924618723 0.053066483 0.017102610 0.005212184

Marc in der Kiste
quelle

+1. Das Beispiel des Iris-Datensatzes ist interessant, da bei Betrachtung der ersten beiden PCs (wie z. B. in Ihrem eigenen Beitrag hier stats.stackexchange.com/a/88092 ) klar ist, dass der zweite ein Signal überträgt . Der Permutationstest (auch Mischen genannt) zeigt jedoch, dass nur der erste Test "signifikant" ist. Es ist klar, dass das Mischen dazu neigt, die Anzahl der PCs zu unterschätzen: Die große Varianz des ersten realen PCs wird sich über die gemischten PCs "ausbreiten" und alle ab dem zweiten erhöhen. Man kann empfindlichere Tests entwickeln, die das erklären, aber dies wird selten gemacht.

Amöbe sagt Reinstate Monica

@amoeba - Ausgezeichneter Kommentar. Ich wundere mich schon seit einiger Zeit über den "Ausbreitungseffekt". Ich nehme an, ein Kreuzvalidierungstest ist einer der heikeleren Tests, auf die Sie verweisen (z. B. Ihre Antwort hier ). Wäre toll, wenn Sie ein Beispiel / Referenz zur Verfügung stellen könnten.

Marc in der Box

Normalerweise bevorzuge ich die Kreuzvalidierung (basierend auf dem Rekonstruktionsfehler, wie in meiner Antwort hier angegeben ), aber ich bin mir nicht sicher, ob es nicht irgendwie unter der gleichen Art von Unempfindlichkeit leidet oder nicht. Könnte sinnvoll sein, es mit dem Iris-Datensatz auszuprobieren. Was die auf Mischen basierenden Ansätze angeht, so kenne ich keine Anhaltspunkte für diese "Ausbreitung". Ich kenne nur einige Leute, die in letzter Zeit daran gearbeitet haben. Ich glaube, sie wollten es bald aufschreiben. Die Idee ist, einige Downscaling-Faktoren für die Varianzen von höher gemischten PCs einzuführen.

Amöbe sagt Reinstate Monica

@amoeba - Danke für diesen Link. Das erklärt mir viel. Ich fand es besonders interessant zu sehen, dass die Kreuzvalidierung in PCA Methoden verwendet, die mit Datensätzen mit fehlenden Werten arbeiten können. Ich habe einige Versuche mit diesem Ansatz unternommen, und (wie Sie sagen) der Nullmodell-Shuffling-Ansatz unterschätzt in der Tat die Anzahl der wichtigen PCs. Für das Iris-Dataset gebe ich jedoch konsistent einen einzelnen PC für Rekonstruktionsfehler zurück. Interessant angesichts dessen, was Sie über die Handlung gesagt haben. Könnte sein, dass, wenn wir Fehler basierend auf Artenvorhersagen messen, die Ergebnisse unterschiedlich sein könnten.

Marc in der Box

Aus Neugier probierte ich es an den Iris-Daten aus. Tatsächlich bekomme ich zwei bedeutende PCs mit der Kreuzvalidierungsmethode. Ich habe meinen verlinkten Beitrag aktualisiert, siehe dort.

Amöbe sagt Reinstate Monica