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]]