Пошаговое внедрение PCA в R с использованием учебника Линдси Смит

Я работаю в R с помощью превосходного учебника по PCA Линдси и Смита, и застреваю на последнем этапе. Сценарий R, приведенный ниже, выводит нас на этап (на стр.19), на котором исходные данные восстанавливаются из (в данном случае, единственного) основного компонента, который должен давать прямую линию вдоль оси PCA1 (учитывая, что данные имеет только 2 измерения, второе из которых намеренно отбрасывается).

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)

Это насколько я знаю, и пока все в порядке. Но я не могу понять, как данные получены для окончательного графика - дисперсия, относящаяся к PCA 1 - которую Смит готовит как:

Это то, что я пробовал (что игнорирует добавление оригинальных средств):

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)

.. и получил ошибочный

.. потому что я как-то потерял измерение данных в умножении матриц. Я был бы очень благодарен за идею, что здесь происходит не так.

* Редактировать *

Интересно, это правильная формула:

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

Но я немного сбит с толку, если это так, потому что (а) я понимаю, что rowVectorFeatureнеобходимо уменьшить до желаемой размерности (собственный вектор для PCA1), и (b) он не совпадает с абзацем PCA1:

Любые мнения высоко ценятся.

Вы были очень близко к этому и попали в ловушку тонкой проблемы при работе с матрицами в R. Я работал с вашим final_dataи получил правильные результаты независимо. Затем я присмотрелся к вашему коду. Короче говоря, где вы написали

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

ты был бы в порядке, если бы написал

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

trans_data$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,])

$2\times1$$1\times10$final_data$20=2\times10$row_orig_data$12=2\times1 + 1\times10$

$(XY)^T=Y^TX^T$t(t(p) %*% t(q)) = q %*% t

$x/y$$y/x$

Написать

d_in_new_basis = as.matrix(final_data)

затем, чтобы вернуть ваши данные в исходную базу вам нужно

d_in_original_basis = d_in_new_basis %*% feat_vec

Вы можете обнулить части ваших данных, которые проецируются вдоль второго компонента, используя

d_in_new_basis_approx = d_in_new_basis
d_in_new_basis_approx[,2] = 0

и вы можете преобразовать, как раньше

d_in_original_basis_approx = d_in_new_basis_approx %*% feat_vec

Нанесение их на один и тот же график вместе с зеленой линией главного компонента показывает, как работает аппроксимация.

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

Давайте вернемся к тому, что у вас было. Эта линия была в порядке

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

feat_vec %*% row_data_adj$Y=S^TX$$S$$X$$Y$$Y$$X$$Y$$X$

Тогда у вас было

trans_data = final_data
trans_data[,2] = 0

Это нормально: вы просто обнуляете части своих данных, которые проецируются вдоль второго компонента. Где это идет не так

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

$\hat Y$$Y$$\mathbf{e}_1$t(feat_vec[1,]) %*% t(trans_data)$\mathbf{e}_1 \hat Y$

$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$

Я думаю, что у вас есть правильная идея, но вы наткнулись на неприятную особенность R. Вот снова соответствующий фрагмент кода, как вы это сформулировали:

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)

По существу final_dataсодержит координаты исходных точек относительно системы координат, определяемой собственными векторами ковариационной матрицы. Поэтому для восстановления исходных точек необходимо умножить каждый собственный вектор на соответствующую преобразованную координату, например

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

который дал бы исходные координаты первой точки. В своем вопросе вы установите второй компонент правильно к нулю trans_data[,2] = 0. Если вы тогда (как вы уже редактировали) рассчитать

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

Вы рассчитываете формулу (1) для всех точек одновременно. Ваш первый подход

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

вычисляет что-то другое и работает только потому, что R автоматически удаляет атрибут измерения feat_vec[1,], поэтому он больше не является вектором строки, а рассматривается как вектор столбца. Последующее преобразование снова превращает его в вектор строки, и по этой причине, по крайней мере, вычисление не приводит к ошибке, но если вы пройдете математику, вы увидите, что это нечто иное, чем (1). В целом, при умножении матриц хорошей идеей является подавление отбрасывания атрибута измерения, которое может быть достигнуто dropпараметром, например feat_vec[1,,drop=FALSE].

$\Delta y / \Delta x$

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

Изучив это упражнение, вы можете попробовать более простые способы в R. Существуют две популярные функции для выполнения PCA: princompи prcomp. princompФункция делает собственное значение разложения , как это делалось в тренировках. prcompФункция использует разложение по сингулярным значениям. Оба метода будут давать одинаковые результаты почти все время: этот ответ объясняет различия в R, тогда как этот ответ объясняет математику . (Спасибо TooTone за комментарии, теперь интегрированные в этот пост.)

Здесь мы используем оба, чтобы воспроизвести упражнение в R. Сначала используем 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)

Второе использование 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)

Очевидно, что знаки перевернуты, но объяснение вариации эквивалентно.