PCA程序

资源来自网上搜集。

python代码，调用的是sklearn库里面的PCA函数。

#! /usr/bin/env python
#coding=utf-8
import numpy as np
from sklearn.decomposition import PCA
X = np.array([[-1,2,66,-1], [-2,6,58,-1], [-3,8,45,-2], [1,9,36,1], [2,10,62,1], [3,5,83,2]])  #导入数据，维度为4
pca = PCA(n_components=2)   #降到2维
pca.fit(X)                  #训练
newX=pca.fit_transform(X)   #降维后的数据
# PCA(copy=True, n_components=2, whiten=False)
print("导入数据，维度为4")
print(X)
print("输出贡献率")
print(pca.explained_variance_ratio_)  #输出贡献率
print("输出降维后的数据")
print(newX)                  #输出降维后的数据

运行结果：

导入数据，维度为4
[[-1  2 66 -1]
 [-2  6 58 -1]
 [-3  8 45 -2]
 [ 1  9 36  1]
 [ 2 10 62  1]
 [ 3  5 83  2]]
输出贡献率
[ 0.95713353  0.03398198]
输出降维后的数据
[[  7.96504337   4.12166867]
 [ -0.43650137   2.07052079]
 [-13.63653266   1.86686164]
 [-22.28361821  -2.32219188]
 [  3.47849303  -3.95193502]
 [ 24.91311585  -1.78492421]]

 python不调用封装好的PCA函数

#! /usr/bin/env python
#coding=utf-8
from numpy import *
import numpy
def pca(X,CRate):
    #矩阵X每行是一个样本
     #对样本矩阵进行中心化样本矩阵
     meanValue=mean(X,axis=0)#计算每列均值
     X=X-meanValue#每个维度元素减去对应维度均值
     #协方差矩阵
     C=cov(X,rowvar=0)
     #特征值，特征向量
     eigvalue,eigvector=linalg.eig(mat(C))#特征值，特征向量
     #根据贡献率，来决定取多少个特征向量构成变换矩阵
     sumEigValue=sum(eigvalue)#所有特征值之和
     sortedeigvalue= numpy.sort(eigvalue)[::-1]    #对特征值从大到小排序
     for i in range(sortedeigvalue.size):
        j=i+1
        rate=sum(eigvalue[0:j])/sumEigValue
        if rate>CRate:
            break
     #取前j个列向量构成变换矩阵
     indexVec=numpy.argsort(-eigvalue)    #对covEigenVal从大到小排序，返回索引
     nLargestIndex=indexVec[:j] #取出最大的特征值的索引
     T=eigvector[:,nLargestIndex] #取出最大的特征值对应的特征向量
     newX=numpy.dot(X,T)#将X矩阵降维得到newX
     return newX,T,meanValue#返回降维后矩阵newX，变换矩阵T，每列的均值构成的数组



if __name__ == "__main__":
    
    X = numpy.array([[-1,2,66,-1], [-2,6,58,-1], [-3,8,45,-2], [1,9,36,1], [2,10,62,1], [3,5,83,2]])  #导入数据，维度为4
    rate1 = 0.96
    [newX,T,meanValue] = pca(X,rate1)
    
    print("newX")
    print(newX)
    print("T")
    print(T)
    print("meanValue")
    print(meanValue)

运行结果：

newX
[[  7.96504337  -4.12166867]
 [ -0.43650137  -2.07052079]
 [-13.63653266  -1.86686164]
 [-22.28361821   2.32219188]
 [  3.47849303   3.95193502]
 [ 24.91311585   1.78492421]]
T
[[ 0.06761155  0.60143206]
 [-0.09934165  0.68922828]
 [ 0.99207082  0.01304165]
 [ 0.03681576  0.40382394]]
meanValue
[  0.           6.66666667  58.33333333   0.        ]

 matlab代码实现PCA

function [newX,T,meanValue] = pca_row(X,CRate)
% 来自博客：https://www.cnblogs.com/simon-c/p/4902651.html
%每行是一个样本
%newX  降维后的新矩阵
%T 变换矩阵
%meanValue  X每列均值构成的矩阵，用于将降维后的矩阵newX恢复成X
%CRate 贡献率
%计算中心化样本矩阵
meanValue=ones(size(X,1),1)*mean(X);
X=X-meanValue;%每个维度减去该维度的均值
C=X'*X/(size(X,1)-1);%计算协方差矩阵

%计算特征向量，特征值
[V,D]=eig(C);
%将特征向量按降序排序
[dummy,order]=sort(diag(-D));
V=V(:,order);%将特征向量按照特征值大小进行降序排列
d=diag(D);%将特征值取出，构成一个列向量
newd=d(order);%将特征值构成的列向量按降序排列

%取前n个特征向量，构成变换矩阵
sumd=sum(newd);%特征值之和
for j=1:length(newd)
    i=sum(newd(1:j,1))/sumd;%计算贡献率，贡献率=前n个特征值之和/总特征值之和
    if i>CRate%当贡献率大于95%时循环结束,并记下取多少个特征值
        cols=j;
        break;
    end
end
T=V(:,1:cols);%取前cols个特征向量，构成变换矩阵T
newX=X*T;%用变换矩阵T对X进行降维
end
%测试程序
%test=[10 15 29;15 46 13;23 21 30;11 9 35;42 45 11;9 48 5;11 21 14;8 5 15;11 12 21;21 20 25]
%[newX,T,meanValue]=pca_row(test,0.9)
% 将降维后得到的新矩阵newX恢复：
% 
% 公式为X=newX*T'+meanValue

在命令行窗口运行并输出：

>> test=[10 15 29;15 46 13;23 21 30;11 9 35;42 45 11;9 48 5;11 21 14;8 5 15;11 12 21;21 20 25]

test =

    10    15    29
    15    46    13
    23    21    30
    11     9    35
    42    45    11
     9    48     5
    11    21    14
     8     5    15
    11    12    21
    21    20    25

>> [newX,T,meanValue]=pca_row(test,0.9)

newX =

   13.4627   -0.1472
  -21.2616   -6.1205
    4.7222   11.1751
   20.7366    4.1128
  -29.3539   16.6403
  -24.3452  -15.3551
    2.0237   -6.9416
   17.2018   -7.6807
   12.5972   -2.8162
    4.2167    7.1330


T =

   -0.3025    0.8750
   -0.8672   -0.0881
    0.3956    0.4760


meanValue =

   16.1000   24.2000   19.8000
   16.1000   24.2000   19.8000
   16.1000   24.2000   19.8000
   16.1000   24.2000   19.8000
   16.1000   24.2000   19.8000
   16.1000   24.2000   19.8000
   16.1000   24.2000   19.8000
   16.1000   24.2000   19.8000
   16.1000   24.2000   19.8000
   16.1000   24.2000   19.8000

参考博客：https://blog.csdn.net/reticent_man/article/details/82633214

参考博客：https://www.cnblogs.com/simon-c/p/4934305.html

参考博客：https://www.cnblogs.com/simon-c/p/4902651.html