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 = \begin{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]]


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM