Czy w przypadku macierzy losowej SVD nie powinien niczego wyjaśniać? Co ja robię źle?

Gdybym zbudował matrycę 2-D złożoną wyłącznie z losowych danych, oczekiwałbym, że komponenty PCA i SVD w zasadzie niczego nie wyjaśnią.

Zamiast tego wydaje się, że pierwsza kolumna SVD wydaje się wyjaśniać 75% danych. Jak to możliwe? Co ja robię źle?

Oto fabuła:

wprowadź opis zdjęcia tutaj

Oto kod R:

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)

Aktualizacja

Dziękuję @Aaron. Jak zauważyłeś, poprawka polegała na dodaniu skalowania do macierzy, aby liczby były wyśrodkowane wokół 0 (tj. Średnia wynosi 0).

m <- scale(m, scale=FALSE)

Oto poprawiony obraz, pokazujący macierz z losowymi danymi, pierwsza kolumna SVD jest bliska 0, zgodnie z oczekiwaniami.

Poprawiony obraz

r pca svd

— Contango
źródło

Twoja macierz jest w przybliżeniu równomiernie rozłożona na kostce jednostkowej w . SVD oblicza momenty bezwładności na temat pochodzenia . W „całkowita wariancja” musi być razy większa od interwału jednostkowego, czyli . Łatwo obliczyć, że moment wzdłuż głównej osi sześcianu (emanujący z początku) wynosi a wszystkie pozostałe momenty - z racji symetrii - równe . Dlatego pierwsza wartość własna wynosi całości. Dla to

[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$ %, widoczne na trzecim wykresie.

— whuber

Odpowiedzi:

Pierwszy komputer wyjaśnia, że zmienne nie są wyśrodkowane wokół zera. Skalowanie najpierw lub centrowanie zmiennych losowych wokół zera da oczekiwany wynik. Na przykład jeden z tych:

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 opuścił Stack Overflow
źródło

Podnosisz dobrą rację, ale myślę, że to tylko część historii. Rzeczywiście, zgaduję, że PO spróbuje każdego z nich i nadal będzie zaskoczony wynikiem. Faktem jest, że ponieważ liczby pojedyncze są uporządkowane w danych wyjściowych, nie pojawią się (a nawet nie będą) równomiernie rozmieszczone, jak można naiwnie oczekiwać na podstawie „losowych” danych. W tym przypadku rozkład Marchenko-Pastur reguluje ich zachowanie.

— kardynał

@Aaron Dziękuję, miałeś absolutną rację. Dodałem wykres skorygowanego wyjścia powyżej, aby pokazać, jak piękny jest wynik.

— Contango,

@cardinal Dziękuję za komentarz, masz absolutną rację (patrz wykres utworzony przez poprawiony kod powyżej). Wierzę, że wartości SVD stałyby się mniej równomiernie rozłożone w miarę zmniejszania się macierzy, ponieważ mniejsza matryca miałaby większe szanse na posiadanie wzorów, które nie byłyby zmiażdżone przez prawo wielkich liczb.

— Contango,

Dodam bardziej wizualną odpowiedź do twojego pytania, poprzez porównanie z zerowym modelem. Procedura losowo tasuje dane w każdej kolumnie, aby zachować ogólną wariancję, podczas gdy kowariancja między zmiennymi (kolumnami) zostaje utracona. Odbywa się to kilka razy, a wynikowy rozkład liczby pojedynczej w losowej macierzy jest porównywany z wartościami pierwotnymi.

Używam prcompzamiast svddo rozkładu macierzy, ale wyniki są podobne:

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

Porównanie modeli zerowych jest wykonywane na macierzy centrowanej poniżej:

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)

Poniżej znajduje się wykres pudełkowy permutowanej macierzy z 95% kwantylem każdej liczby pojedynczej pokazanej jako linia ciągła. Pierwotne wartości PCA mto kropki. wszystkie leżą poniżej linii 95% - dlatego ich amplituda jest nie do odróżnienia od szumu przypadkowego.

Tę samą procedurę można wykonać w niecentralnej wersji mz tym samym wynikiem - Brak znaczących pojedynczych wartości:

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

Dla porównania spójrzmy na zestaw danych z nieprzypadkowym zestawem danych: 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)

Tutaj pierwsza wartość pojedyncza jest znacząca i wyjaśnia ponad 92% całkowitej wariancji:

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

— Marc w pudełku
źródło

+1. Przykład zestawu danych Iris jest interesujący, ponieważ patrząc na pierwsze dwa komputery osobiste (jak np. We własnym poście tutaj stats.stackexchange.com/a/88092 ) jest całkiem jasne, że drugi niesie jakiś sygnał. Test permutacji (inaczej tasowanie) wskazuje, że tylko pierwszy jest „znaczący”. Oczywiste jest, że tasowanie ma tendencję do niedoceniania liczby komputerów: duża wariancja pierwszego prawdziwego komputera zostanie „rozłożona” na tasowane komputery i podniesie je wszystkie, zaczynając od drugiego. Można opracować bardziej czułe testy, które to uwzględniają, ale jest to rzadko wykonywane.

— ameba mówi Przywróć Monikę

@amoeba - Doskonały komentarz. Od jakiegoś czasu zastanawiam się nad efektem „rozprzestrzeniania się”. Przypuszczam, że test weryfikacji krzyżowej może być jednym z bardziej wrażliwych, na który się powołujesz (np. Twoja odpowiedź tutaj )? Byłoby świetnie, gdybyś mógł podać przykład / odniesienie.

— Marc w pudełku

Zwykle wolę korzystać z walidacji krzyżowej (na podstawie błędu rekonstrukcji, zgodnie z moją odpowiedzią tutaj ), ale tak naprawdę nie jestem pewien, czy nie cierpi z powodu tego samego rodzaju niewrażliwości, czy nie. Może warto wypróbować to w zestawie danych Iris. Jeśli chodzi o podejście oparte na tasowaniu, nie znam żadnego odniesienia do rozliczania tego „rozprzestrzeniania się”, po prostu znam ludzi, którzy ostatnio nad tym pracowali. Myślę, że wkrótce chcieli to napisać. Chodzi o to, aby wprowadzić pewne czynniki zmniejszania skali dla wariantów wyższych komputerów z tasowaniem.

— ameba mówi Przywróć Monikę

@amoeba - Dzięki za ten link. Wiele mi to wyjaśnia. Za szczególnie interesujące uważam, że krzyżowa walidacja w PCA wykorzystuje metody, które mogą działać na zestawach danych z brakującymi wartościami. Podjąłem kilka prób zastosowania tego podejścia i (jak twierdzisz) podejście tasowania modelu zerowego rzeczywiście nie docenia liczby znaczących komputerów. Jednak w przypadku zestawu danych tęczówki konsekwentnie zwracam pojedynczy komputer z powodu błędu rekonstrukcji. Interesujące, biorąc pod uwagę to, co wspomniałeś o fabule. Możliwe, że gdybyśmy oceniali błąd w oparciu o prognozy dotyczące gatunków, wyniki mogą być inne.

— Marc w pudełku

Z ciekawości wypróbowałem to na danych Iris. W rzeczywistości otrzymuję dwa znaczące komputery z metodą krzyżowej weryfikacji. Zaktualizowałem mój połączony post, zobacz tam.

— ameba mówi Przywróć Monikę