复制 Shalizi 的纽约时报 PCA 示例

机器算法验证 主成分分析 Python 降维 svd
2022-03-28 12:42:53

我正在尝试复制 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 的图片中。

有任何想法吗?

1个回答

事实证明,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)