Estou trabalhando no R através de um excelente tutorial de PCA de Lindsay I Smith e estou ficando paralisado na última etapa. O script R abaixo nos leva ao estágio (na p.19) em que os dados originais estão sendo reconstruídos a partir do (principal neste caso) Componente Principal, que deve gerar um gráfico de linha reta ao longo do eixo PCA1 (dado que os dados possui apenas 2 dimensões, sendo a segunda eliminada intencionalmente).
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)
Isso é o máximo que eu tenho, e tudo bem até agora. Mas não consigo descobrir como os dados são obtidos para o gráfico final - a variação atribuível ao PCA 1 - que Smith plota como:
Isto é o que eu tentei (que ignora a adição dos meios originais):
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)
.. e tem um erro:
.. porque eu perdi uma dimensão de dados de alguma forma na multiplicação da matriz. Ficaria muito grato por ter uma ideia do que está acontecendo de errado aqui.
* Editar *
Gostaria de saber se esta é a fórmula certa:
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')
Mas estou um pouco confuso se sim, porque (a) entendo que as rowVectorFeature
necessidades devem ser reduzidas à dimensionalidade desejada (o vetor próprio para PCA1) e (b) não está alinhado com o abline PCA1:
Qualquer vista muito apreciada.
s1
Respostas:
Você estava muito perto e foi pego por um problema sutil ao trabalhar com matrizes na R. Eu trabalhei com você
final_data
e obtive os resultados corretos independentemente. Então, olhei mais de perto o seu código. Para encurtar uma longa história, onde você escreveuvocê teria ficado bem se tivesse escrito
trans_data
t(feat_vec[1,])
row_orig_data = t(as.matrix(feat_vec[1,],ncol=1,nrow=2) %*% t(trans_data))
non-conformable arguments
final_data
row_orig_data
t(t(p) %*% t(q)) = q %*% t
Escrever
para recuperar seus dados em sua base original, você precisa
Você pode zerar as partes de seus dados projetadas ao longo do segundo componente usando
e você pode transformar como antes
A plotagem no mesmo gráfico, junto com a linha do componente principal em verde, mostra como a aproximação funcionou.
Vamos voltar ao que você tinha. Esta linha estava ok
feat_vec %*% row_data_adj
Então você teve
Tudo bem: você está zerando as partes dos dados projetadas ao longo do segundo componente. Onde isso dá errado é
t(feat_vec[1,]) %*% t(trans_data)
fonte
Eu acho que você tem a idéia certa, mas tropeçou em um recurso desagradável de R. Aqui, novamente, a parte do código relevante como você a declarou:
Essencialmente
final_data
contém as coordenadas dos pontos originais em relação ao sistema de coordenadas definido pelos vetores próprios da matriz de covariância. Para reconstruir os pontos originais, é necessário multiplicar cada vetor próprio pela coordenada transformada associada, por exemploo que produziria as coordenadas originais do primeiro ponto. Na sua pergunta, você define o segundo componente corretamente como zero
trans_data[,2] = 0
,. Se você (como você já editou) calculavocê calcula a fórmula (1) para todos os pontos simultaneamente. Sua primeira abordagem
calcula algo diferente e só funciona porque R elimina automaticamente o atributo de dimensão para
feat_vec[1,]
, portanto, não é mais um vetor de linha, mas tratado como um vetor de coluna. A transposição subsequente o transforma em um vetor de linha novamente e é por isso que pelo menos o cálculo não produz um erro, mas se você passar pela matemática, verá que é algo diferente de (1). Em geral, é uma boa idéia em multiplicações de matrizes suprimir a queda do atributo de dimensão que pode ser alcançado pelodrop
parâmetro, por exemplofeat_vec[1,,drop=FALSE]
.fonte
drop=F
argumento.Depois de explorar este exercício, você pode tentar as maneiras mais fáceis em R. Existem duas funções populares para executar o PCA:
princomp
eprcomp
. Aprincomp
função faz a decomposição do autovalor, como foi feito no seu exercício. Aprcomp
função usa decomposição de valor singular. Ambos os métodos fornecerão os mesmos resultados quase o tempo todo: esta resposta explica as diferenças em R, enquanto que essa resposta explica a matemática . (Obrigado ao TooTone pelos comentários agora integrados a esta postagem.)Aqui usamos os dois para reproduzir o exercício em R. Primeiro, usando
princomp
:Segundo uso
prcomp
:Claramente, os sinais são invertidos, mas a explicação da variação é equivalente.
fonte