看到網上有一堆“博客”,明顯是抄襲的,前后矛盾,自己摸索着寫了一個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]]