标签: pca

假设我有一个数据矩阵d

pc = prcomp(d)

# pc1 and pc2 are the principal components  
pc1 = pc$rotation[,1] 
pc2 = pc$rotation[,2]

那么这应该适合线性回归模型吗？

r = lm(y ~ pc1+pc2)

但后来我得到了这个错误:

Errormodel.frame.default(formula = y ~ pc1+pc2, drop.unused.levels = TRUE) : 
   unequal dimensions('pc1')

我想那里有一个自动执行此操作的软件包,但这也应该有用吗？

我正在使用Python,我已经使用本教程实现了PCA .

一切都很好,我得到了协方差,我做了一个成功的转换,使它对原始尺寸没有问题.

但是我该如何进行美白？我尝试用特征值划分特征向量:

S, V = numpy.linalg.eig(cov)
V = V / S[:, numpy.newaxis]

并使用V来转换数据,但这导致了奇怪的数据值.请问有人可以对此有所了解吗？

假设我有一些文本句子,我想用kmeans进行聚类.

sentences = [
    "fix grammatical or spelling errors",
    "clarify meaning without changing it",
    "correct minor mistakes",
    "add related resources or links",
    "always respect the original author"
]

from sklearn.feature_extraction.text import CountVectorizer
from sklearn.cluster import KMeans

vectorizer = CountVectorizer(min_df=1)
X = vectorizer.fit_transform(sentences)
num_clusters = 2
km = KMeans(n_clusters=num_clusters, init='random', n_init=1,verbose=1)
km.fit(X)

现在我可以预测新文本会落入哪个类,

new_text = "hello world"
vec = vectorizer.transform([new_text])
print km.predict(vec)[0]

但是,我说应用PCA将10,000个功能减少到50个.

from sklearn.decomposition import RandomizedPCA

pca = RandomizedPCA(n_components=50,whiten=True)
X2 = pca.fit_transform(X)
km.fit(X2)

我不能再做同样的事情来预测新文本的集群,因为矢量化器的结果不再相关

new_text = "hello world"
vec = …

我正在创建一个可视化来说明主成分分析如何工作,通过绘制一些实际数据的特征值(为了说明的目的,我将子集化为2维).

我想要从这个梦幻般的PCA教程中结合这两个图,仅用于我的真实数据.

我可以绘制矢量并且一切正常:

Person1 <- c(-3,1,1,-3,0,-1,-1,0,-1,-1,3,4,5,-2,1,2,-2,-1,1,-2,1,-3,4,-6,1,-3,-4,3,3,-5,0,3,0,-3,1,-2,-1,0,-3,3,-4,-4,-7,-5,-2,-2,-1,1,1,2,0,0,2,-2,4,2,1,2,2,7,0,3,2,5,2,6,0,4,0,-2,-1,2,0,-1,-2,-4,-1)
Person2 <- c(-4,-3,4,-5,-1,-1,-2,2,1,0,3,2,3,-4,2,-1,2,-1,4,-2,6,-2,-1,-2,-1,-1,-3,5,2,-1,3,3,1,-3,1,3,-3,2,-2,4,-4,-6,-4,-7,0,-3,1,-2,0,2,-5,2,-2,-1,4,1,1,0,1,5,1,0,1,1,0,2,0,7,-2,3,-1,-2,-3,0,0,0,0)
df <- data.frame(cbind(Person1, Person2))
g <- ggplot(data = df, mapping = aes(x = Person1, y = Person2))
g <- g + geom_point(alpha = 1/3)  # alpha b/c of overplotting
g <- g + geom_smooth(method = "lm")  # just for comparsion
g <- g + coord_fixed()  # otherwise, the angles of vectors are off
corre <- cor(x = df$Person1, y = df$Person2, method = …

我正在尝试使用KernelPCA将数据集的维度降低到2D(用于可视化目的和进一步的数据分析).

我尝试使用不同Gamma值的RBF内核计算KernelPCA,但结果不稳定:

(每帧的Gamma值略有不同,其中Gamma从0到1连续变化)

看起来它不是确定性的.

有没有办法稳定它/使其确定性？

用于生成转换数据的代码:

def pca(X, gamma1):
    kpca = KernelPCA(kernel="rbf", fit_inverse_transform=True, gamma=gamma1)
    X_kpca = kpca.fit_transform(X)
    #X_back = kpca.inverse_transform(X_kpca)
    return X_kpca

我有一个100,000行×27,000列的csv,我试图在PCA上生成100,000行X 300列矩阵.csv大9GB.这是我正在做的事情:

from sklearn.decomposition import PCA as RandomizedPCA
import csv
import sys
import numpy as np
import pandas as pd

dataset = sys.argv[1]
X = pd.DataFrame.from_csv(dataset)
Y = X.pop("Y_Level")
X = (X - X.mean()) / (X.max() - X.min())
Y = list(Y)
dimensions = 300
sklearn_pca = RandomizedPCA(n_components=dimensions)
X_final = sklearn_pca.fit_transform(X)

当我运行上面的代码时,我的程序在执行.from_csv步骤时被杀死.我已经能够通过将csv分成10,000组来解决这个问题; 逐个读取它们,然后调用pd.concat.这允许我在被杀之前进入标准化步骤(X-X.mean()).... 我的数据对我的macbook空间来说太大了吗？或者有更好的方法来做到这一点.我真的很想将我拥有的所有数据用于我的机器学习应用程序.

如果我想按照下面的答案建议使用增量PCA,我会这样做吗？:

from sklearn.decomposition import IncrementalPCA
import csv
import sys
import numpy as np
import pandas as pd

dataset = sys.argv[1]
chunksize_ = 10000
#total_size is 100000
dimensions …

我想N=10,000从矩阵中选择最主要的组件.在完成pca之后,MATLAB应该返回一个pxp矩阵,但它不会!

>> size(train_data)
ans =
         400      153600

>> [coefs,scores,variances] = pca(train_data);
>> size(coefs)
ans =
      153600         399

>> size(scores)
ans =

   400   399
>> size(variances)
ans =
    399     1

应该是coefs:153600 x 153600？和scores:400 X 153600？

当我使用下面的代码时,它给我一个Out of Memory错误::

>> [V D] = eig(cov(train_data));
Out of memory. Type HELP MEMORY for your options.

Error in cov (line 96)
    xy = (xc' * xc) / (m-1);

我不明白为什么MATLAB会返回一个较小的维矩阵.它应该返回错误与pca:153600*153600*8字节= 188 GB

eigs错误:

>> eigs(cov(train_data));
Out of memory. …

我哪里错了？我试图通过prcomp和我自己进行PCA,我得到不同的结果,你能帮助我吗？

由我自己做的:

>database <- read.csv("E:/R/database.csv", sep=";", dec=",") #it's a 105 rows x 8 columns, each column is a variable
>matrix.cor<-cor(database)
>standardize<-function(x) {(x-mean(x))/sd(x)}
>values.standard<-apply(database, MARGIN=2, FUN=standardize)
>my.eigen<-eigen(matrix.cor)
>loadings<-my.eigen$vectors
>scores<-values.standard %*% loadings
>head (scores, n=10) # I m just posting here the first row scores for the first 6 pc

[,1]       [,2]       [,3]        [,4]       [,5]        [,6]        

2.3342586  2.3426398 -0.9169527  0.80711713  1.1409138 -0.25832090    

>sd <-sqrt (my.eigen$values)
>sd

[1] 1.5586078 1.1577093 1.1168477 0.9562853 0.8793033 0.8094500 0.6574788
0.4560247

与PRCOMP一起做:

>database.pca<-prcomp(database, retx=TRUE, center= TRUE, scale=TRUE)
>sd1<-database.pca$sdev …

我希望能够使用其加载来构建主成分分析的分数,但我无法弄清楚princomp函数在计算数据集的分数时实际上在做什么.玩具示例:

cc <- matrix(1:24,ncol=4)
PCAcc <- princomp(cc,scores=T,cor=T)
PCAcc$loadings

Loadings:
     Comp.1 Comp.2 Comp.3 Comp.4
[1,]  0.500  0.866              
[2,]  0.500 -0.289  0.816       
[3,]  0.500 -0.289 -0.408 -0.707
[4,]  0.500 -0.289 -0.408  0.707

PCAcc$scores

       Comp.1        Comp.2        Comp.3 Comp.4
[1,] -2.92770 -6.661338e-16 -3.330669e-16      0
[2,] -1.75662 -4.440892e-16 -2.220446e-16      0
[3,] -0.58554 -1.110223e-16 -6.938894e-17      0
[4,]  0.58554  1.110223e-16  6.938894e-17      0
[5,]  1.75662  4.440892e-16  2.220446e-16      0
[6,]  2.92770  6.661338e-16  3.330669e-16      0

我的理解是分数是负载和重新缩放的原始数据的线性组合.尝试"手":

rescaled <- t(t(cc)-apply(cc,2,mean))
rescaled%*%PCAcc$loadings

     Comp.1        Comp.2        Comp.3 Comp.4
[1,]     -5 -1.332268e-15 -4.440892e-16      0
[2,] …

我有一个300万x 900万的稀疏矩阵,有数十亿个非零项.R和Python不允许稀疏矩阵具有超过MAXINT非零的条目,因此我发现自己使用Julia.

虽然使用标准偏差来缩放这些数据是微不足道的,但是贬低当然是一种天真的方式,因为这会产生一个密集的200+太字节矩阵.

有关svd的相关代码是julia,请访问https://github.com/JuliaLang/julia/blob/343b7f56fcc84b20cd1a9566fd548130bb883505/base/linalg/arnoldi.jl#L398

从我的阅读来看,这段代码的一个关键元素是AtA_or_AAt结构和围绕这些结构的几个函数,特别是A_mul_B!.为方便起见,下面复制

struct AtA_or_AAt{T,S} <: AbstractArray{T, 2}
    A::S
    buffer::Vector{T}
end

function AtA_or_AAt(A::AbstractMatrix{T}) where T
    Tnew = typeof(zero(T)/sqrt(one(T)))
    Anew = convert(AbstractMatrix{Tnew}, A)
    AtA_or_AAt{Tnew,typeof(Anew)}(Anew, Vector{Tnew}(max(size(A)...)))
end

function A_mul_B!(y::StridedVector{T}, A::AtA_or_AAt{T}, x::StridedVector{T}) where T
    if size(A.A, 1) >= size(A.A, 2)
        A_mul_B!(A.buffer, A.A, x)
        return Ac_mul_B!(y, A.A, A.buffer)
    else
        Ac_mul_B!(A.buffer, A.A, x)
        return A_mul_B!(y, A.A, A.buffer)
    end
end
size(A::AtA_or_AAt) = ntuple(i -> min(size(A.A)...), Val(2))
ishermitian(s::AtA_or_AAt) = true

这被传递到eigs函数中,在那里发生了一些魔术,然后输出被处理到SVD的相关组件.

我认为使这个"快速居中"类型设置工作的最佳方法是使用AtA_or_AAT_centered版本的子类AtA_or_AAT,或多或少模仿行为,但也存储列方式,并重新定义A_mul_B!功能恰当.

但是,我并没有非常使用朱莉娅,并且已经遇到了一些难以修改的问题.在我再次尝试深入研究之前,我想知道如果这被认为是一个合适的攻击计划,或者如果在这么大的矩阵上进行SVD​​的简单方法,我是否能得到反馈(我没有看到它,但我可能错过了一些东西).

编辑:我没有修改基础Julia,而是尝试编写一个"中心稀疏矩阵"包来保持输入稀疏矩阵的稀疏结构,但是在各种计算中适当地进入列.它的实施范围有限,而且有效.不幸的是,它仍然太慢,尽管为了优化事物做了一些非常广泛的努力.