複製 Shalizi 的紐約時報 PCA 示例

我正在嘗試複製 Shalizi 在他的Advanced Data Analysis with an Elementary Point of View書中找到的 NY Times PCA 示例。我在這裡找到了示例代碼和數據

http://www.stat.cmu.edu/~cshalizi/490/pca/

我加載了 R 工作區 pca-examples.Rdata，並能夠從中提取一個 CSV，並將其放入

https://github.com/burakbayramli/classnotes/blob/master/app-math-tr/pca/nytimes.csv

為了使用 SVD 進行 PCA，我做了

from pandas import *
import numpy.linalg as lin
nyt = read_csv ("nytimes.csv")
labels = nyt['class.labels']
nyt = nyt.drop(["class.labels"],axis=1)
nyt_demeaned = nyt - nyt.mean(0)
u,s,v = lin.svd(nyt_demeaned.T,full_matrices=False)

然後我嘗試投影到由特徵向量定義的空間上

nyt2 = np.dot(nyt_demeaned, u[:,0:2])

然後繪製

plot(nyt2[:,0],nyt2[:,1],'.')

和

arts = nyt2[labels == 'art']
music = nyt2[labels == 'music']
plot(arts[:,0],arts[:,1],'r.')
plot(music[:,0],music[:,1],'b.')

我明白了

這張圖片與沙利子書中的圖片不同（使用以下R代碼）

load("pca-examples.Rdata")
nyt.pca = prcomp(nyt.frame[,-1])
nyt.latent.sem = nyt.pca$rotation
plot(nyt.pca$x[,1:2],type="n")
points(nyt.pca$x[nyt.frame[,"class.labels"]=="art",1:2],
      pch="A",col="red")
points(nyt.pca$x[nyt.frame[,"class.labels"]=="music",1:2],
      pch="M",col="blue")

也許有一個 90 度的錯誤，以一種或另一種方式翻轉，等等。然後你可能會讓兩個圖像有點接近，但即使這樣分佈也不完全正確，藝術的數據點也沒有那麼清楚地分開來自音樂數據點，因為它們在 Shalizi 的圖片中。

有任何想法嗎？

事實證明，Shalizi 使用此代碼進行了一些額外的歸一化，稱為逆文檔頻率加權

scale.cols <- function(x,s) {
 return(t(apply(x,1,function(x){x*s})))
}

div.by.euc.length <- function(x) {
 scale.rows(x,1/sqrt(rowSums(x^2)+1e-16))
}

idf.weight <- function(x) {
 # IDF weighting
 doc.freq <- colSums(x>0)
 doc.freq[doc.freq == 0] <- 1
 w <- log(nrow(x)/doc.freq)
 return(scale.cols(x,w))
}

load("pca-examples.Rdata")
nyt.frame.raw$class.labels <- NULL
nyt.frame <- idf.weight(nyt.frame.raw)
nyt.frame <- div.by.euc.length(nyt.frame)

當我在 nyt.frame 而不是 nyt.frame.raw 上運行 SVD 時（加載後，不需要上面的代碼，它只是用於演示），這兩個類是明顯可分離的。點的分佈仍然不像沙裡子的照片，但我想我可以用這個。

注意：僅僅讀取 nytimes4.csv 標準化文件是不夠的，數據仍然需要以 0 為中心（貶義）。

注意：在 nytimes4.csv 中跳過刪除虛假列並貶低後，輸出與 Shalizi 的完全相同。

Python / Pandas 中的相同內容

from pandas import *
nyt = read_csv ("nytimes.csv")
nyt = nyt.drop(['class.labels'],axis=1)
freq = nyt.astype(bool).sum(axis=0)
freq = freq.replace(0,1)
w = np.log(float(nyt.shape[0])/freq)
nyt = nyt.apply(lambda x: x*w,axis=1)
nyt = nyt.apply(lambda x: x / np.sqrt(np.sum(np.square(x))+1e-16), axis=1)

引用自：https://stats.stackexchange.com/questions/48214