• PCA的计算方法


    看到网上有一堆“博客”,明显是抄袭的,前后矛盾,自己摸索着写了一个PCA的计算过程。
    假设有5个学生的6门功课:语文、数学、地理、化学、英语、历史,成绩如下:

    X = np.array(
        [[84,65,61,72,79,81],
         [64,77,77,76,55,70],
         [65,67,63,49,57,67],
         [74,80,69,75,63,74],
         [84,74,70,80,74,82]])
    

    注意,行是样本(表示一个学生),列是特征(表示一门课)。

    首先要搞明白什么是协方差。定义:(下面的n是样本数)

    均值(假设权重概率都为1):

    [mu = frac{1}{n}sum_{i=1}^n x_i ]

    标准差(除以n-1表示无偏估计):

    [std = sqrt{frac{1}{n-1}sum_{i=1}^n(x_i-mu)^2} ]

    方差:

    [var = std^2 = frac{1}{n-1}sum_{i=1}^n(x_i-mu)^2 ]

    协方差:两个特征之间的方差

    [cov(X,Y) = frac{1}{n-1}sum_{i=1}^n (x_i - mu_x)(y_i-mu_y) ]

    也就是计算所有样本的语文成绩与数学成绩之间的方差,或者化学成绩与英语成绩之间的方差。

    用python实现方差协方差计算

    def my_mean(data):
        return np.sum(data) / len(data)
    
    def cov(a,b):
        assert(len(a) == len(b))
        mean_a = my_mean(a)
        mean_b = my_mean(b)
        p = (a - mean_a)
        q = (b - mean_b)
        r = np.dot(p,q.T)
        return r/(len(a)-1)
    

    协方差矩阵:多个特征之间的方差的矩阵

    [c = egin{pmatrix} cov(x,x) & cov(x,y) & cov(x,z) \ cov(y,x) & cov(y,y) & cov(y,z) \ cov(z,x) & cov(z,y) & cov(z,z) \ end{pmatrix} ]

    可以用上面的函数来计算,当然numpy也有现成的函数:np.cov()。使用这个函数时,注意矩阵的形状。以上面的X的数据为例,一共5个样本6个特征,按照协方差矩阵的定义,应该是一个6x6的方阵。
    如果用np.cov(X)计算出来的是一个5x5的方阵,不满足条件。仔细阅读cov()函数的文档,X的例子是列为特征,所以应该用:

    cc = np.cov(X, rowvar=False)
    #或者 
    dd = np.cov(X.T)
    
    cc= 
    [[ 95.2  -13.9  -23.75  62.15 100.35  63.05]
     [-13.9   41.3   32.75  44.95 -26.95  -5.1 ]
     [-23.75  32.75  40.    42.5  -33.    -8.5 ]
     [ 62.15  44.95  42.5  151.3   53.7   53.85]
     [100.35 -26.95 -33.    53.7  110.8   65.9 ]
     [ 63.05  -5.1   -8.5   53.85  65.9   43.7 ]]
    
    dd= 
    [[ 95.2  -13.9  -23.75  62.15 100.35  63.05]
     [-13.9   41.3   32.75  44.95 -26.95  -5.1 ]
     [-23.75  32.75  40.    42.5  -33.    -8.5 ]
     [ 62.15  44.95  42.5  151.3   53.7   53.85]
     [100.35 -26.95 -33.    53.7  110.8   65.9 ]
     [ 63.05  -5.1   -8.5   53.85  65.9   43.7 ]]
    

    当然,也可以用协方差矩阵的定义来求:

    c1 = np.array(
        [[cov(x[:,0],x[:,0]),cov(x[:,0],x[:,1]),cov(x[:,0],x[:,2]),cov(x[:,0],x[:,3]),cov(x[:,0],x[:,4]),cov(x[:,0],x[:,5])],
         [cov(x[:,1],x[:,0]),cov(x[:,1],x[:,1]),cov(x[:,1],x[:,2]),cov(x[:,1],x[:,3]),cov(x[:,1],x[:,4]),cov(x[:,1],x[:,5])],
         [cov(x[:,2],x[:,0]),cov(x[:,2],x[:,1]),cov(x[:,2],x[:,2]),cov(x[:,2],x[:,3]),cov(x[:,2],x[:,4]),cov(x[:,2],x[:,5])],
         [cov(x[:,3],x[:,0]),cov(x[:,3],x[:,1]),cov(x[:,3],x[:,2]),cov(x[:,3],x[:,3]),cov(x[:,3],x[:,4]),cov(x[:,3],x[:,5])],
         [cov(x[:,4],x[:,0]),cov(x[:,4],x[:,1]),cov(x[:,4],x[:,2]),cov(x[:,4],x[:,3]),cov(x[:,4],x[:,4]),cov(x[:,4],x[:,5])],
         [cov(x[:,5],x[:,0]),cov(x[:,5],x[:,1]),cov(x[:,5],x[:,2]),cov(x[:,5],x[:,3]),cov(x[:,5],x[:,4]),cov(x[:,5],x[:,5])]])
    print("c1=",c1)
    
    c1= 
    [[ 95.2  -13.9  -23.75  62.15 100.35  63.05]
     [-13.9   41.3   32.75  44.95 -26.95  -5.1 ]
     [-23.75  32.75  40.    42.5  -33.    -8.5 ]
     [ 62.15  44.95  42.5  151.3   53.7   53.85]
     [100.35 -26.95 -33.    53.7  110.8   65.9 ]
     [ 63.05  -5.1   -8.5   53.85  65.9   43.7 ]]
    

    可以看到cc,dd,c1三者是一样的。

    得到协方差矩阵后,使用函数求其特征值和特征向量:

    eig_value, eig_vector = np.linalg.eig(cc)
    
    eig_value: 
    [3.06293191e+02 1.63510310e+02 9.89302953e+00 2.60347035e+00 3.77194607e-15 1.37974973e-14]
    
    eig_vector: 
    [[-0.53264253  0.20279107 -0.34433806  0.39437042 -0.61869481 -0.55543331]
     [ 0.00876193 -0.46059524 -0.81597078  0.02185232  0.25842516  0.34848844]
     [ 0.04593605 -0.47328385  0.37877077  0.70892582 -0.03144886  0.21014772]
     [-0.51955599 -0.64238594  0.24891406 -0.45230979 -0.15412561 -0.22434743]
     [-0.55131936  0.32775478  0.09651389 -0.13044526  0.29446728  0.67491022]
     [-0.37445103  0.05145202  0.0297077   0.34614812  0.66255449  0.14160509]]
    

    因为要降维,原来是6维,我们降为2维,所以取eig_value中两个最大的值:

    eig_sort_index = np.argsort(eig_value)
    # 得到结果 [4,5,3,2,1,0] 是个从小到大的排序数组下标,表示第0个值最大,第1个值其次,第4个值最小
    eig_big_index = eig_sort_index[:-3:-1] # 两个最大值,[0,1]
    eig_big_vector = eig_vector[:,eig_big_index]
    

    得到eig_big_vector的结果:

    array([[-0.53264253,  0.20279107],
           [ 0.00876193, -0.46059524],
           [ 0.04593605, -0.47328385],
           [-0.51955599, -0.64238594],
           [-0.55131936,  0.32775478],
           [-0.37445103,  0.05145202]])
    

    这实际上是两个特征向量的组合,每一列是一组特征向量。

    然后得到X的mean均值矩阵:

    mean = np.mean(x, axis=0)
    print(mean)
    x_mean = x - mean
    print("x_mean=",x_mean)
    
    mean=[74.2 72.6 68.  70.4 65.6 74.8]
    x_mean= 
    [[  9.8  -7.6  -7.    1.6  13.4   6.2]
     [-10.2   4.4   9.    5.6 -10.6  -4.8]
     [ -9.2  -5.6  -5.  -21.4  -8.6  -7.8]
     [ -0.2   7.4   1.    4.6  -2.6  -0.8]
     [  9.8   1.4   2.    9.6   8.4   7.2]]
    

    最后用特征向量与x_mean相乘,得到PCA的降维结果:

    r = np.dot(x_mean, eig_big_vector)
    print("r=", r)
    
    r= 
    [[-16.14860528  12.48396235]
     [ 10.61676743 -15.67317428]
     [ 23.40212697  13.607117  ]
     [ -0.43966353  -7.77054621]
     [-17.43062559  -2.64735885]]
    
  • 相关阅读:
    知识搜索
    使用 getopt() 进行命令行处理
    【新提醒】夏新大v安卓4.1尝鲜最新更新版本发布(包含进步版)1124更新 大V综合交流区 360论坛
    搜狗知立方高调亮相 开启知识计算新时代
    socat: Linux / UNIX TCP Port Forwarder
    Crontab 每两周执行一次
    python 命令行解析 optionparser
    crontab jojo's blog--快乐忧伤都与你同在 BlogJava
    搜索引擎开始「实体搜索」新领域竞争,Google、百度分别发力实体搜索产品
    netcat(nc)的替代产品 Socat
  • 原文地址:https://www.cnblogs.com/woodyh5/p/12299914.html
Copyright © 2020-2023  润新知