Krok po kroku wdrożenie PCA w języku R przy użyciu samouczka Lindsay Smith

13

Pracuję w R poprzez doskonały samouczek PCA autorstwa Lindsay I Smith i utknąłem w ostatnim etapie. Poniższy skrypt R przenosi nas do etapu (na str. 19), na którym odtwarzane są oryginalne dane z (w tym przypadku pojedynczego) głównego elementu, który powinien dać wykres linii prostej wzdłuż osi PCA1 (biorąc pod uwagę, że dane ma tylko 2 wymiary, z których drugi jest celowo upuszczany).

d = data.frame(x=c(2.5,0.5,2.2,1.9,3.1,2.3,2.0,1.0,1.5,1.1),
               y=c(2.4,0.7,2.9,2.2,3.0,2.7,1.6,1.1,1.6,0.9))

# mean-adjusted values 
d$x_adj = d$x - mean(d$x)
d$y_adj = d$y - mean(d$y)

# calculate covariance matrix and eigenvectors/values
(cm = cov(d[,1:2]))

#### outputs #############
#          x         y
# x 0.6165556 0.6154444
# y 0.6154444 0.7165556
##########################

(e = eigen(cm))

##### outputs ##############
# $values
# [1] 1.2840277 0.0490834
#
# $vectors
#          [,1]       [,2]
# [1,] 0.6778734 -0.7351787
# [2,] 0.7351787  0.6778734
###########################


# principal component vector slopes
s1 = e$vectors[1,1] / e$vectors[2,1] # PC1
s2 = e$vectors[1,2] / e$vectors[2,2] # PC2

plot(d$x_adj, d$y_adj, asp=T, pch=16, xlab='x', ylab='y')
abline(a=0, b=s1, col='red')
abline(a=0, b=s2)

wprowadź opis zdjęcia tutaj

# PCA data = rowFeatureVector (transposed eigenvectors) * RowDataAdjust (mean adjusted, also transposed)
feat_vec = t(e$vectors)
row_data_adj = t(d[,3:4])
final_data = data.frame(t(feat_vec %*% row_data_adj)) # ?matmult for details
names(final_data) = c('x','y')

#### outputs ###############
# final_data
#              x           y
# 1   0.82797019 -0.17511531
# 2  -1.77758033  0.14285723
# 3   0.99219749  0.38437499
# 4   0.27421042  0.13041721
# 5   1.67580142 -0.20949846
# 6   0.91294910  0.17528244
# 7  -0.09910944 -0.34982470
# 8  -1.14457216  0.04641726
# 9  -0.43804614  0.01776463
# 10 -1.22382056 -0.16267529
############################

# final_data[[1]] = -final_data[[1]] # for some reason the x-axis data is negative the tutorial's result

plot(final_data, asp=T, xlab='PCA 1', ylab='PCA 2', pch=16)

wprowadź opis zdjęcia tutaj

To jest tak daleko, jak mam, i wszystko w porządku do tej pory. Ale nie mogę zrozumieć, w jaki sposób uzyskuje się dane dla ostatecznego wykresu - wariancji przypisywanej PCA 1 - którą Smith wykreśla jako:

wprowadź opis zdjęcia tutaj

Oto, co próbowałem (co ignoruje dodawanie oryginalnych środków):

trans_data = final_data
trans_data[,2] = 0
row_orig_data = t(t(feat_vec[1,]) %*% t(trans_data))
plot(row_orig_data, asp=T, pch=16)

.. i otrzymałem błąd:

wprowadź opis zdjęcia tutaj

.. ponieważ jakoś straciłem wymiar danych podczas mnożenia macierzy. Byłbym bardzo wdzięczny za pomysł, co się tutaj dzieje.


* Edytować *

Zastanawiam się, czy jest to właściwa formuła:

row_orig_data = t(t(feat_vec) %*% t(trans_data))
plot(row_orig_data, asp=T, pch=16, cex=.5)
abline(a=0, b=s1, col='red')

Ale jestem trochę zdezorientowany, jeśli tak, ponieważ (a) rozumiem rowVectorFeaturepotrzebę zredukowania do pożądanej wymiarowości (wektor własny dla PCA1) i (b) nie pokrywa się z linią PCA1:

wprowadź opis zdjęcia tutaj

Wszelkie opinie bardzo cenione.

geotheory
źródło
Krótka notatka (już wspomniana w odpowiedziach poniżej, ale potencjalnie myląca dla kogoś, kto patrzy na twoje pytanie): twoje s1nachylenie zostało obliczone z pomyłką (powinno być , a nie ), dlatego czerwona linia nie jest idealnie dopasowane do danych na pierwszej figurze i rekonstrukcji na ostatniej. x / yy/xx/y
ameba mówi Przywróć Monikę
Jeśli chodzi o rekonstrukcję oryginalnych danych z głównych głównych komponentów, zobacz ten nowy wątek: stats.stackexchange.com/questions/229092 .
ameba mówi Przywróć Monikę

Odpowiedzi:

10

Byłeś bardzo, bardzo blisko i zostałeś złapany przez subtelny problem w pracy z macierzami w R. Pracowałem od ciebie final_datai uzyskałem prawidłowe wyniki niezależnie. Potem przyjrzałem się twojemu kodowi. Krótko mówiąc, gdzie napisałeś

row_orig_data = t(t(feat_vec[1,]) %*% t(trans_data))

byłbyś w porządku, gdybyś napisał

row_orig_data = t(t(feat_vec) %*% t(trans_data))

trans_data2×12×10t(feat_vec[1,])1×2row_orig_data = t(as.matrix(feat_vec[1,],ncol=1,nrow=2) %*% t(trans_data))non-conformable arguments

row_orig_data = t(as.matrix(feat_vec[1,],ncol=1,nrow=2) %*% t(trans_data)[1,])

2×11×10final_data20=2×10row_orig_data12=2×1+1×10

(XY)T=YTXTt(t(p) %*% t(q)) = q %*% t

x/yy/x


pisać

d_in_new_basis = as.matrix(final_data)

następnie, aby odzyskać swoje dane w pierwotnej podstawie, potrzebujesz

d_in_original_basis = d_in_new_basis %*% feat_vec

Za pomocą można wyzerować części danych, które są rzutowane wzdłuż drugiego komponentu

d_in_new_basis_approx = d_in_new_basis
d_in_new_basis_approx[,2] = 0

i możesz następnie przekształcić jak poprzednio

d_in_original_basis_approx = d_in_new_basis_approx %*% feat_vec

Umieszczenie ich na tym samym wykresie wraz z główną linią składową na zielono pokazuje, jak działało przybliżenie.

plot(x=d_in_original_basis[,1]+mean(d$x),
     y=d_in_original_basis[,2]+mean(d$y),
     pch=16, xlab="x", ylab="y", xlim=c(0,3.5),ylim=c(0,3.5),
     main="black=original data\nred=original data restored using only a single eigenvector")
points(x=d_in_original_basis_approx[,1]+mean(d$x),
       y=d_in_original_basis_approx[,2]+mean(d$y),
       pch=16,col="red")
points(x=c(mean(d$x)-e$vectors[1,1]*10,mean(d$x)+e$vectors[1,1]*10), c(y=mean(d$y)-e$vectors[2,1]*10,mean(d$y)+e$vectors[2,1]*10), type="l",col="green")

wprowadź opis zdjęcia tutaj

Wróćmy do tego, co miałeś. Ta linia była w porządku

final_data = data.frame(t(feat_vec %*% row_data_adj))

feat_vec %*% row_data_adjY=STXSXYYXYX

Więc miałeś

trans_data = final_data
trans_data[,2] = 0

Jest w porządku: zerujesz tylko te części danych, które są rzutowane wzdłuż drugiego komponentu. Gdzie idzie źle

row_orig_data = t(t(feat_vec[1,]) %*% t(trans_data))

Y^Ye1t(feat_vec[1,]) %*% t(trans_data)e1Y^

2×12×10Y^Yy1e1y1ie1y1e1i

TooTone
źródło
Dzięki TooTone jest to bardzo obszerne i rozwiązuje dwuznaczności w moim rozumieniu obliczeń macierzy i roli FeatureVector na ostatnim etapie.
geotheory
Wspaniały :). Odpowiedziałem na to pytanie, ponieważ obecnie studiuję teorię SVD / PCA i chciałem zrozumieć, jak to działa na przykładzie: twoje pytanie było w odpowiednim momencie. Po przejściu wszystkich obliczeń macierzy byłem nieco zaskoczony, że okazał się to problem R - więc cieszę się, że doceniłeś również aspekt macierzy.
TooTone
4

Myślę, że masz dobry pomysł, ale natknąłeś się na nieprzyjemną cechę R. Ponownie odpowiedni fragment kodu, jak już powiedziałeś:

trans_data = final_data
trans_data[,2] = 0
row_orig_data = t(t(feat_vec[1,]) %*% t(trans_data))
plot(row_orig_data, asp=T, pch=16)

Zasadniczo final_datazawiera współrzędne pierwotnych punktów w odniesieniu do układu współrzędnych określonego przez wektory własne macierzy kowariancji. Aby zrekonstruować oryginalne punkty, należy zatem pomnożyć każdy wektor własny z powiązaną transformowaną współrzędną, np

(1) final_data[1,1]*t(feat_vec[1,] + final_data[1,2]*t(feat_vec[2,])

co dałoby oryginalne współrzędne pierwszego punktu. W swoim pytaniu drugi składnik prawidłowo ustawić na zero trans_data[,2] = 0. Jeśli następnie (jak już edytowałeś) oblicz

(2) row_orig_data = t(t(feat_vec) %*% t(trans_data))

obliczasz wzór (1) dla wszystkich punktów jednocześnie. Twoje pierwsze podejście

row_orig_data = t(t(feat_vec[1,]) %*% t(trans_data))

oblicza coś innego i działa tylko dlatego, że R automatycznie upuszcza atrybut wymiaru feat_vec[1,], więc nie jest to już wektor wiersza, ale traktowany jak wektor kolumny. Kolejna transpozycja sprawia, że ​​znów jest to wektor wiersza i dlatego przynajmniej obliczenia nie powodują błędu, ale jeśli przejdziesz przez matematykę, zobaczysz, że jest to coś innego niż (1). Zasadniczo dobrym pomysłem w przypadku mnożenia macierzy jest tłumienie upuszczania atrybutu wymiaru, który można osiągnąć za pomocą dropparametru, np feat_vec[1,,drop=FALSE].

Δy/Δx

s1 = e$vectors[2,1] / e$vectors[1,1] # PC1
s2 = e$vectors[2,2] / e$vectors[1,2] # PC2
Georg Schnabel
źródło
Dziękuję bardzo Georg. Masz rację co do stoku PCA1. Bardzo przydatna wskazówka również na temat drop=Fargumentu.
geotheory
4

Po zapoznaniu się z tym ćwiczeniem możesz wypróbować łatwiejsze sposoby w R. Istnieją dwie popularne funkcje do wykonywania PCA: princompi prcomp. Ta princompfunkcja rozkłada wartość własną jak w ćwiczeniu. prcompWykorzystuje rozkładu wartości pojedynczej. Obie metody dają te same wyniki prawie przez cały czas: ta odpowiedź wyjaśnia różnice w R, podczas gdy ta odpowiedź wyjaśnia matematykę . (Dzięki TooTone za komentarze teraz zintegrowane z tym postem).

Tutaj używamy obu do odtworzenia ćwiczenia w R. Po pierwsze princomp:

d = data.frame(x=c(2.5,0.5,2.2,1.9,3.1,2.3,2.0,1.0,1.5,1.1), 
               y=c(2.4,0.7,2.9,2.2,3.0,2.7,1.6,1.1,1.6,0.9))

# compute PCs
p = princomp(d,center=TRUE,retx=TRUE)

# use loadings and scores to reproduce with only first PC
loadings = t(p$loadings[,1]) 
scores = p$scores[,1] 

reproduce = scores %*% loadings  + colMeans(d)

# plots
plot(reproduce,pch=3,ylim=c(-1,4),xlim=c(-1,4))
abline(h=0,v=0,lty=3)
mtext("Original data restored using only a single eigenvector",side=3,cex=0.7)

biplot(p)

wprowadź opis zdjęcia tutaj wprowadź opis zdjęcia tutaj

Drugie użycie prcomp:

d = data.frame(x=c(2.5,0.5,2.2,1.9,3.1,2.3,2.0,1.0,1.5,1.1), 
               y=c(2.4,0.7,2.9,2.2,3.0,2.7,1.6,1.1,1.6,0.9))

# compute PCs
p = prcomp(d,center=TRUE,retx=TRUE)

# use loadings and scores to reproduce with only first PC
loadings = t(p$rotation[,1])
scores = p$x[,1]

reproduce = scores %*% loadings  + colMeans(d)

# plots
plot(reproduce,pch=3,ylim=c(-1,4),xlim=c(-1,4))
abline(h=0,v=0,lty=3)
mtext("Original data restored using only a single eigenvector",side=3,cex=0.7)

biplot(p)

wprowadź opis zdjęcia tutaj wprowadź opis zdjęcia tutaj

Oczywiście znaki są odwrócone, ale wyjaśnienie zmienności jest równoważne.

mrbcuda
źródło
Dzięki mrbcuda. Twój dwupłat wygląda identycznie jak Lindsay Smith, więc zakładam, że zastosował tę samą metodę 12 lat temu! Zdaję sobie również sprawę z niektórych innych metod wyższego poziomu , ale jak słusznie zauważyłeś, jest to ćwiczenie, aby wyjaśnić podstawowe matematyki PCA.
geotheory