Langkah demi langkah implementasi PCA di R menggunakan tutorial Lindsay Smith

Saya bekerja di R melalui tutorial PCA yang sangat baik oleh Lindsay I Smith dan saya terjebak di tahap terakhir. Script R di bawah ini membawa kita ke tahap (pada hal.19) di mana data asli sedang direkonstruksi dari Komponen Utama (tunggal dalam hal ini), yang akan menghasilkan plot garis lurus sepanjang sumbu PCA1 (mengingat bahwa data hanya memiliki 2 dimensi, yang kedua sedang sengaja dihapus).

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)

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

Sejauh ini yang saya miliki, dan semuanya baik-baik saja sejauh ini. Tapi saya tidak tahu bagaimana data diperoleh untuk plot terakhir - varian yang dikaitkan dengan PCA 1 - yang oleh Smith diplot sebagai:

Inilah yang saya coba (yang mengabaikan menambahkan cara asli):

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)

.. dan mendapat erron:

.. karena saya telah kehilangan dimensi data entah bagaimana dalam perkalian matriks. Saya akan sangat berterima kasih atas ide apa yang salah di sini.

* Edit *

Saya ingin tahu apakah ini formula yang tepat:

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

Tapi saya agak bingung jika demikian karena (a) Saya mengerti rowVectorFeaturekebutuhan harus direduksi ke dimensi yang diinginkan (vektor eigen untuk PCA1), dan (b) tidak sejalan dengan PCA1 abline:

Setiap pandangan sangat dihargai.

Anda hampir sangat di sana dan tertangkap oleh masalah halus dalam bekerja dengan matriks di R. Saya bekerja melalui Anda final_datadan mendapatkan hasil yang benar secara mandiri. Kemudian saya melihat lebih dekat kode Anda. Singkatnya, di mana Anda menulis

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

Anda akan baik-baik saja jika Anda telah menulis

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

sebagai gantinya (karena Anda akan memusatkan perhatian pada bagian trans_datayang diproyeksikan pada vektor eigen kedua). Karena Anda mencoba mengalikan matriks dengan matriks tetapi R tidak memberi Anda kesalahan. Masalahnya adalah bahwa diperlakukan sebagai . Mencoba akan memberi Anda kesalahan. Berikut ini, mungkin lebih sesuai dengan apa yang Anda maksudkan, juga akan berhasil$2\times1$$2\times10$t(feat_vec[1,])$1\times2$row_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,])

karena mengalikan matriks dengan matriks (perhatikan bahwa Anda bisa menggunakan matriks asli di sini). Tidak perlu melakukannya dengan cara ini, tetapi lebih baik secara matematis karena ini menunjukkan bahwa Anda mendapatkan nilai dari nilai di sisi kanan.$2\times1$$1\times10$final_data$20=2\times10$row_orig_data$12=2\times1 + 1\times10$

Saya telah meninggalkan jawaban asli saya di bawah, karena seseorang mungkin menganggapnya berguna, dan itu menunjukkan mendapatkan plot yang diperlukan. Ini juga menunjukkan bahwa kode bisa sedikit lebih sederhana dengan menyingkirkan beberapa transpos yang tidak perlu: begitu .$(XY)^T=Y^TX^T$t(t(p) %*% t(q)) = q %*% t

Re edit Anda, saya telah menambahkan garis komponen utama berwarna hijau ke plot saya di bawah ini. Dalam pertanyaan Anda, Anda memiliki kemiringan sebagai bukan .$x/y$$y/x$

Menulis

d_in_new_basis = as.matrix(final_data)

kemudian untuk mendapatkan data Anda kembali ke basis aslinya yang Anda butuhkan

d_in_original_basis = d_in_new_basis %*% feat_vec

Anda dapat menghapus bagian-bagian dari data Anda yang diproyeksikan di sepanjang komponen kedua menggunakan

d_in_new_basis_approx = d_in_new_basis
d_in_new_basis_approx[,2] = 0

dan Anda kemudian dapat mengubah seperti sebelumnya

d_in_original_basis_approx = d_in_new_basis_approx %*% feat_vec

Memplot ini pada plot yang sama, bersama dengan garis komponen utama berwarna hijau, menunjukkan kepada Anda bagaimana perkiraannya bekerja.

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

Mari kita mundur ke apa yang Anda miliki. Baris ini ok

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

Bit penting di sini adalah feat_vec %*% row_data_adjyang setara dengan mana adalah matriks vektor eigen dan adalah matriks data Anda dengan data Anda di baris, dan adalah data dalam basis baru. Apa yang dikatakan ini adalah bahwa baris pertama adalah jumlah (baris ditimbang oleh vektor eigen pertama). Dan baris kedua adalah jumlah (deretan ditimbang oleh vektor eigen kedua).$Y=S^TX$$S$$X$$Y$$Y$$X$$Y$$X$

Lalu kamu punya

trans_data = final_data
trans_data[,2] = 0

Ini tidak masalah: Anda hanya memusatkan perhatian pada bagian data Anda yang diproyeksikan di sepanjang komponen kedua. Di mana itu salah adalah

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

Menulis untuk matriks data di basis baru, dengan nol di baris kedua, dan menulis untuk vektor eigen pertama, bagian bisnis dari kode ini turun ke .$\hat Y$$Y$$\mathbf{e}_1$t(feat_vec[1,]) %*% t(trans_data)$\mathbf{e}_1 \hat Y$

Seperti dijelaskan di atas (ini adalah di mana saya menyadari masalah R halus dan menulis bagian pertama dari jawaban saya), secara matematis Anda mencoba untuk memperbanyak suatu vektor oleh matriks. Ini tidak berfungsi secara matematis. Yang harus Anda lakukan adalah mengambil baris pertama = baris pertama : panggil ini . Kemudian, gandakan dan bersamaan. The th kolom hasil adalah vektor eigen ditimbang dengan 1 koordinat hanya dari titik th di dasar baru, yang adalah apa yang Anda inginkan.$2\times1$$2\times10$$\hat Y$$Y$$\mathbf{y}_1$$\mathbf{e}_1$$\mathbf{y}_1$$i$$\mathbf{e}_1\mathbf{y}_1$$\mathbf{e}_1$$i$

Saya pikir Anda memiliki ide yang tepat tetapi tersandung fitur jahat R. Di sini sekali lagi potongan kode yang relevan seperti yang Anda nyatakan:

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)

Pada dasarnya final_databerisi koordinat titik-titik asli sehubungan dengan sistem koordinat yang ditentukan oleh vektor eigen dari matriks kovarian. Untuk merekonstruksi titik-titik asli kita harus mengalikan masing-masing vektor eigen dengan koordinat transformasi yang terkait, misalnya

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

yang akan menghasilkan koordinat asli dari titik pertama. Di pertanyaan Anda, Anda mengatur komponen kedua dengan benar ke nol trans_data[,2] = 0,. Jika Anda kemudian (seperti yang sudah diedit) menghitung

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

Anda menghitung rumus (1) untuk semua poin secara bersamaan. Pendekatan pertama Anda

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

menghitung sesuatu yang berbeda dan hanya berfungsi karena R secara otomatis menjatuhkan atribut dimensi untuk feat_vec[1,], jadi itu bukan vektor baris lagi tetapi diperlakukan sebagai vektor kolom. Transpos berikutnya membuatnya menjadi vektor baris lagi dan itulah alasan mengapa setidaknya perhitungan tidak menghasilkan kesalahan, tetapi jika Anda membaca matematika, Anda akan melihat bahwa itu adalah sesuatu yang berbeda dari (1). Secara umum itu adalah ide yang baik dalam perkalian matriks untuk menekan menjatuhkan atribut dimensi yang dapat dicapai oleh dropparameter, misalnya feat_vec[1,,drop=FALSE].

Solusi Anda yang diedit tampaknya benar, tetapi Anda menghitung kemiringan jika PCA1 salah. Kemiringan diberikan oleh , karenanya$\Delta y / \Delta x$

s1 = e$vectors[2,1] / e$vectors[1,1] # PC1
s2 = e$vectors[2,2] / e$vectors[1,2] # PC2

Setelah menjelajahi latihan ini, Anda dapat mencoba cara-cara yang lebih mudah di R. Ada dua fungsi populer untuk melakukan PCA: princompdan prcomp. The princompfungsi melakukan dekomposisi eigenvalue seperti yang dilakukan dalam latihan Anda. The prcompFungsi menggunakan dekomposisi nilai singular. Kedua metode akan memberikan hasil yang sama hampir sepanjang waktu: jawaban ini menjelaskan perbedaan R, sedangkan jawaban ini menjelaskan matematika . (Terima kasih kepada TooTone untuk komentar yang sekarang terintegrasi ke dalam posting ini.)

Di sini kami menggunakan keduanya untuk mereproduksi latihan di R. Pertama menggunakan 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)

Kedua menggunakan 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)

Jelas tanda-tanda dibalik tetapi penjelasan variasi sama.