机器学习—PCA
6.PCA
6.1什么是PCA
6.1.1 PCA的实现思想
主成分分析法(Principal Componentm Analysis)是一个非监督的机器学习算法，主要用于数据降维、去噪声等，简称PCA。
问题1：对于一个二维的特征空间，如何将其降维到一维的特征空间？
解决1：在二维特征空间中找到一条直线，将特征空间中所有的样本点映射到这条直线上。
问题2：降维过程中，数据会发生丢失，如何保持数据的完整性？即如何找到这条直线能最大程度的保证数据的丢失率最低？
解决2：对于二维特征空间，样本可以映射到任意一条直线；在映射直线上，样本间的间距最大时，表示数据保存最完整；这条映射直线称为最佳映射直线。在数据统计中，使用方差来表示数据的离散程度，因此在映射直线上，其样本特征的方差越大，表示样本越离散，说明数据越完整。方差：
问题3：如何找到这条映射直线，使样本特征的方差达到最大值？
解决3：(1)对所有样本特征的均值归零，即demean，特征空间的样本分布如下所示：
(2)求映射直线的方向w=(w1,w2)，使样本映射到w后，样本的方差达到最大值，即：
由于对所有样本进行demean处理，因此样本的均值为0；Xproject是一个向量，因此其距离需要求向量的模；Xproject是在映射直线上的向量，而不是在二维特征空间上的向量，如下所示。
(3)Xproject=(Xpro1,Xpro2)，X=(X1,X2)，w=(w1,w2)并且w是单位向量，推导如下：
即求取w，使有最大值。
问题4：在n维特征空间，该如何进行降维？
解决4：将问题3中的公式推导到n维，如下所示：
求取上式中，w为何值时，有最大值，可以使用梯度上升法来解决。
问题5：如何使用梯度上升法？
解决5：使用迭代公式，其中参数都和梯度下降法一致。
问题6：如何使用梯度上升法来求取
的最大值？
解决6：推导公式如下
其中，
6.1.2自己实现PCA算法
#程序6-1
import numpy as np
import matplotlib.pyplot as plt
X = np.empty(shape=(200,2))
X[:,0] = np.random.randint(0,1000,size=200)
X[:,1] = X[:,0]*0.5 + 600 + np.random.normal(6,30,size=200)
#demean
def demean(X):
return X - np.mean(X,axis=0)
X_demean = demean(X)
#梯度上升法
def f_x(X,w):
return sum(X.dot(w)**2)/len(X)
def df_x(X,w):
return 2.*(X.T.dot(X.dot(w))) / len(X)
#向量单位化
def unit_vector(w):
return w/np.linalg.norm(w)
#求第一主成分:梯度上升法
def first_component(X,init_w,eta,max_iters=1e4):
#w是单位向量
w = unit_vector(init_w)
cur_iters = 0
while cur_iters < max_iters:
grad = df_x(X,w)
last_w = w
w = w + eta*grad
w = unit_vector(w)
if(abs(f_x(X,w)-f_x(X,last_w)) < 1e-8):
break
cur_iters += 1
return w
#不需要进行数据归一化
eta = 0.01
#w不能是0向量
init_w = np.random.random(size=X_demean.shape[1])
w = first_component(X_demean,init_w,eta)
print(w)
#画图
plt.scatter(X_demean[:,0],X_demean[:,1])
plt.plot([-w[0]*550, w[0]*550], [-w[1]*550, w[1]*550], color='r')
plt.show()
运行结果：
[0.89588346 0.44428912]
前面公式使用的X都是去均值化的样本数据集，而w是过原点的单位向量。
6.2数据降维和降噪
6.2.1主成分
在n维的特征空间中，求出的第一个w称为第一主成分；而n维特征空间有n个主成分，如何求出剩下的主成分？
以二维特征空间为例，将样本数据集在第一主成分上的投影称为Xproject=||Xproject||*w=(Xw)*w。其中X表示去均值化的样本数据集，大小为(m,2)；w是过原点的单位向量，大小为(2,1)。将数据集在第一主成分上的分量减去，用变换后的数据集来计算剩下的主成分，即X = X - Xproject = X - (Xw)*w；再用这个X使用梯度上升法计算第二主成分。在n维特征空间，同样如此。
#程序6-2
import numpy as np
import matplotlib.pyplot as plt
X = np.empty(shape=(200,2))
X[:,0] = np.random.randint(0,1000,size=200)
X[:,1] = X[:,0]*0.5 + 600 + np.random.normal(6,30,size=200)
#demean
def demean(X):
return X - np.mean(X,axis=0)
X_demean = demean(X)
#梯度上升法
def f_x(X,w):
return sum(X.dot(w)**2)/len(X)
def df_x(X,w):
return 2.*(X.T.dot(X.dot(w))) / len(X)
#向量单位化
def unit_vector(w):
return w/np.linalg.norm(w)
#求第一主成分:梯度上升法
def first_component(X,init_w,eta,max_iters=1e4):
w = unit_vector(init_w)
cur_iters = 0
while cur_iters < max_iters:
grad = df_x(X,w)
last_w = w
w = w + eta*grad
w = unit_vector(w)
if(abs(f_x(X,w)-f_x(X,last_w)) < 1e-8):
break
cur_iters += 1
return w
#求所有主成分
def all_component(X,init_w,eta,max_iters=1e4):
max_component = X.shape[1]
component = np.empty(shape=(max_component,max_component))
for i in range(max_component):
w = first_component(X, init_w, eta, max_iters)
component[i,:] = w
X = X - X.dot(w).reshape(-1,1)*(w.reshape(1,-1))
return component
eta = 0.01
init_w = np.random.random(size=X_demean.shape[1])
component = all_component(X_demean,init_w,eta)
print(component)
#画图
plt.scatter(X_demean[:,0],X_demean[:,1])
plt.plot([-component[0][0]*550, component[0][0]*550], [-component[0][1]*550, component[0][1]*550], color='r')
plt.plot([-component[1][0]*550, component[1][0]*550], [-component[1][1]*550, component[1][1]*550], color='black')
plt.show()
运行结果：
[[ 0.89500351 0.44605909]
[-0.446059 0.89500356]]
6.2.2高维数据向低维数据映射
在Component列表中依次存放了第一主成分到第n主成分，取前k个主成分存放在Wk中，将去均值化的数据集X从n维变换到k维称为降维。
降维公式如下：
Xk是大小为(m,k)的矩阵，若要将Xk反向映射回原来的特征空间，则X’=XkWk。相比于X，X’会有一定的数据丢失。
#程序6-3
import numpy as np
import matplotlib.pyplot as plt
X = np.empty(shape=(200,2))
X[:,0] = np.random.randint(0,1000,size=200)
X[:,1] = X[:,0]*0.5 + 600 + np.random.normal(6,30,size=200)
#demean
def demean(X):
return X - np.mean(X,axis=0)
X_demean = demean(X)
#梯度上升法
def f_x(X,w):
return sum(X.dot(w)**2)/len(X)
def df_x(X,w):
return 2.*(X.T.dot(X.dot(w))) / len(X)
#向量单位化
def unit_vector(w):
return w/np.linalg.norm(w)
#求第一主成分:梯度上升法
def first_component(X,init_w,eta,max_iters=1e4):
w = unit_vector(init_w)
cur_iters = 0
while cur_iters < max_iters:
grad = df_x(X,w)
last_w = w
w = w + eta*grad
#w是单位向量
w = unit_vector(w)
if(abs(f_x(X,w)-f_x(X,last_w)) < 1e-8):
break
cur_iters += 1
return w
#求所有主成分
def all_component(X,init_w,eta,max_iters=1e4):
max_component = X.shape[1]
component = np.empty(shape=(max_component,max_component))
for i in range(max_component):
w = first_component(X, init_w, eta, max_iters)
component[i,:] = w
X = X - X.dot(w).reshape(-1,1)*(w.reshape(1,-1))
return component
#将数据集X映射到第一主成分的分量中
def mapping_first(X_demean,W_k):
return X_demean.dot(W_k.T)
#将数据集反向映射到原来的特征空间
def mapping_reverse(X_reduction,W_k):
return X_reduction.dot(W_k)
#不需要进行数据归一化
eta = 0.01
#w不能是0向量
init_w = np.random.random(size=X_demean.shape[1])
component = all_component(X_demean,init_w,eta)
#取前k个主成分
#取一个主成分时，需要变换维度，1行n列
W_k = component[0].reshape(1,-1)
#降维的数据集必须去均值化
X_reduction = mapping_first(X_demean,W_k)
#反向映射回原来的特征空间，会损失一部分数据
X_reverse = mapping_reverse(X_reduction,W_k)
#为了和原数据集匹配，需要加上平均值
X_reverse = X_reverse + np.mean(X,axis=0)
#画图
plt.scatter(X[:,0],X[:,1])
plt.scatter(X_reverse[:,0],X_reverse[:,1],color='r')
plt.show()
运行结果：
6.2.3使用sklearn来实现数据降维
#程序6-4
import numpy as np
import matplotlib.pyplot as plt
X = np.empty(shape=(200,2))
X[:,0] = np.random.randint(0,1000,size=200)
X[:,1] = X[:,0]*0.5 + 600 + np.random.normal(6,30,size=200)
#使用sklearn中的PCA
from sklearn.decomposition import PCA
#初始化主成分的个数
pca = PCA(n_components=1)
#训练数据，求出主成分
pca.fit(X)
#打印主成分
print(pca.components_)
#主成分可以解释原数据集的方差率，可以理解为数据集的信息保存率
print(pca.explained_variance_ratio_)
#映射到主成分的分量中,类似mapping_first
#transform的参数X是未去均值化的数据集
X_reduction = pca.transform(X)
#将数据集反向映射到原来的特征空间,类似mapping_reverse
#X_restore也不需要再加上均值
X_restore = pca.inverse_transform(X_reduction)
#画图
plt.scatter(X[:,0],X[:,1])
plt.scatter(X_restore[:,0],X_restore[:,1],color='r')
plt.show()
运行结果：
[[0.89370014 0.44866474]]
[0.99312511]
6.2.4降噪
噪声数据：指数据中存在着错误或异常(偏离期望值)的数据，这些数据对数据的分析造成了干扰。在前面的代码，X[:,1] = X[:,0]*0.5 + 600 + np.random.normal(6,30,size=200)中的np.random.normal(6,30,size=200)就是为了模拟真实数据所加的噪声数据。
降噪：(1)将大小为(m,n)的数据集X映射到主成分分量上(k个主成分)，即Xk；(2)再将Xk反向映射回原来的特征空间。
在程序6-3和6-4已经使用了降噪，以程序6-4为例。对数据集X使用transform映射到主成分分量上；再使用inverse_transform反向映射回原来的特征空间；在图中以红色的点表示。
我们可以通过降噪的方式使数据集更真实、准确。
6.3 PCA的实际使用
6.3.1数据降维后的准确率
#程序6-5
import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
#数字0-9数据集
digits = datasets.load_digits()
print(digits.keys())
print(digits.target_names)
X = digits.data
y = digits.target
print(X.shape)
from sklearn.model_selection import train_test_split
X_train,X_test,y_train,y_test = train_test_split(X,y,test_size=0.2,random_state=1234)
from sklearn.neighbors import KNeighborsClassifier
knn_clf = KNeighborsClassifier(n_neighbors=5)
knn_clf.fit(X_train,y_train)
score = knn_clf.score(X_test,y_test)
print('维数为64的knn算法准确率：',score)
from sklearn.decomposition import PCA
#或 pca = PCA(n_components=0.488)
#n_components为int型时，表示主成分的个数，即所降维数
#n_components为float型时，表示数据信息的保存完整率，
# 0.488正好对应前4个主成分的信息保存率的总和
pca = PCA(n_components=4)
pca.fit(X_train)
X_train_reduction = pca.transform(X_train)
X_test_reduction = pca.transform(X_test)
knn_clf.fit(X_train_reduction,y_train)
score = knn_clf.score(X_test_reduction,y_test)
print('维数为4的knn算法准确率：',score)
print(pca.components_)
print(pca.explained_variance_ratio_)
#当信息保存率为95%时，需要从64维降到28维
pca = PCA(n_components=0.95)
pca.fit(X_train)
X_train_reduction = pca.transform(X_train)
X_test_reduction = pca.transform(X_test)
knn_clf.fit(X_train_reduction,y_train)
score = knn_clf.score(X_test_reduction,y_test)
print('维数为28的knn算法准确率：',score)
print(pca.components_.shape)
print(pca.explained_variance_ratio_)
运行结果：
dict_keys(['data', 'target', 'target_names', 'images', 'DESCR'])
[0 1 2 3 4 5 6 7 8 9]
(1797, 64)
维数为64的knn算法准确率： 0.9833333333333333
维数为4的knn算法准确率： 0.8333333333333334
[[ 4.82562984e-18 1.61193679e-02 2.11800070e-01 1.24654753e-01
4.42768170e-02 1.13706571e-01 1.78431666e-02 -1.42293212e-03
3.91675767e-04 1.17403395e-01 2.35831597e-01 -1.52578968e-01
4.39366029e-02 2.40650510e-01 3.12623216e-02 -2.97158883e-03
8.78247941e-06 8.22096033e-02 -8.85105263e-02 -2.22953850e-01
1.72825600e-01 1.92013824e-01 -1.48423527e-02 -3.14271450e-03
-5.62903998e-05 -6.18731058e-02 -2.59481583e-01 2.17377156e-02
2.16604019e-01 6.30904875e-02 -4.05271458e-02 -1.66167595e-04
-0.00000000e+00 -1.56121295e-01 -3.74186928e-01 -1.72588730e-01
-7.37342398e-02 -2.04592902e-02 -1.44705764e-02 -0.00000000e+00
-1.09454857e-03 -1.02295267e-01 -3.16962922e-01 -2.48265000e-01
-1.84156762e-01 -9.00422121e-03 2.19167345e-02 -1.51636869e-03
-1.80296602e-04 5.75372861e-03 2.94210656e-02 -9.94828373e-02
-9.58973482e-02 1.16975668e-01 4.39728593e-02 -1.92076396e-03
-9.16736174e-06 1.31476957e-02 2.22761634e-01 1.32637904e-01
-4.27879518e-03 6.50081705e-02 2.37774759e-02 9.24669464e-03]
[-7.81093072e-18 -8.71259385e-03 -4.16298279e-02 4.35250748e-05
-5.36347026e-02 -1.27008617e-01 -6.78403835e-02 -8.15510304e-03
-2.23123680e-05 -1.17793766e-02 7.11545444e-02 -9.11317566e-03
-8.70508399e-02 -4.37921320e-02 -6.23123201e-02 -4.14936560e-03
1.03383828e-05 4.52293765e-02 1.94856964e-01 -5.71009989e-02
-2.19698643e-01 -6.73037417e-03 2.30227984e-02 -8.30457796e-04
-1.44588808e-05 7.79321304e-02 1.76523210e-01 -1.39397371e-01
-2.58002218e-01 3.57924923e-02 5.99363601e-02 -2.51060152e-05
0.00000000e+00 8.40001540e-02 7.17492212e-02 -2.74433430e-01
-2.85902712e-01 1.48998660e-01 1.24767609e-01 0.00000000e+00
9.58081387e-05 4.59994120e-02 1.20393908e-01 -2.74275035e-01
-3.11975686e-01 2.28051019e-01 2.19157888e-01 1.25525157e-03
-1.29706123e-04 8.89511641e-03 1.59304233e-01 -1.19757796e-01
-1.08853926e-01 2.87314134e-01 1.56592602e-01 5.23374495e-04
-4.28235061e-05 -9.32706434e-03 -6.41892917e-02 1.77886883e-02
1.99232482e-01 1.88997241e-01 2.63513167e-02 -5.41372043e-03]
[-2.31579643e-17 1.82938675e-02 1.35014862e-01 1.37644805e-01
-1.27776708e-01 -2.57734665e-01 -1.19096722e-01 -1.86647573e-02
2.80064979e-04 8.91519492e-02 9.56730762e-02 4.97254526e-02
8.98865913e-03 -2.62632357e-01 -1.57645319e-01 -1.70572573e-02
1.25657672e-04 4.86020935e-02 -6.58300126e-03 5.94583773e-02
1.27880646e-01 -2.86297187e-01 -1.51330396e-01 -6.98975861e-03
5.90866004e-05 -5.10211906e-02 -9.73928208e-02 8.04672173e-02
2.92994946e-02 -3.44579004e-01 -1.56251149e-01 -2.77618486e-04
-0.00000000e+00 -6.86764811e-02 -6.61290227e-02 1.32928799e-01
3.13123309e-02 -2.64435852e-01 -1.02507713e-01 -0.00000000e+00
-2.71583589e-04 -3.47323450e-02 8.76067505e-02 1.59504792e-01
-8.30554771e-02 -1.43479158e-01 7.22070470e-02 2.43807491e-03
-9.40666482e-05 1.49205959e-02 2.27665142e-01 1.63011913e-01
-5.35935166e-02 9.89787946e-02 2.38630276e-01 2.34266669e-02
-1.63480923e-05 1.72702284e-02 1.33053368e-01 9.53767099e-02
1.39658844e-01 2.43223336e-01 1.74727703e-01 3.69855786e-02]
[ 3.97844758e-18 2.00712866e-02 1.77105191e-01 2.02437232e-01
2.50643073e-02 7.84022598e-02 6.74905655e-02 7.71084131e-03
-2.76372809e-05 7.15668496e-02 2.99572247e-01 3.76844504e-02
-1.40625001e-01 2.62618716e-02 3.77112617e-02 7.96035662e-04
1.22849746e-04 1.14113414e-01 1.96870096e-01 -2.75881870e-01
-2.58601469e-01 -5.38680476e-02 2.10919663e-02 -2.10095211e-03
8.49638587e-05 9.19902999e-02 9.41889308e-02 -2.09966988e-01
-2.23958735e-01 -9.22213694e-02 5.39831025e-02 -3.99159466e-05
-0.00000000e+00 5.12200688e-02 9.49609749e-02 1.66160928e-03
-1.48316958e-01 -1.80012278e-01 8.68433570e-03 -0.00000000e+00
-2.03205915e-04 1.31740745e-02 7.57298443e-02 7.49433131e-02
-2.83085547e-02 -2.07539701e-01 -8.29324805e-02 -2.08117140e-03
1.13032539e-05 1.45488953e-02 1.32182627e-01 1.52051845e-01
-2.50413420e-02 -2.21646931e-01 -1.35739707e-01 -1.61097466e-02
-9.28378963e-06 2.17564236e-02 2.05339983e-01 2.11633880e-01
-2.22823204e-01 -3.02801410e-01 -9.77060246e-02 -2.48306888e-02]]
[0.14859084 0.13632551 0.11814992 0.08496234]
维数为28的knn算法准确率： 0.9833333333333333
(28, 64)
[0.14859084 0.13632551 0.11814992 0.08496234 0.0560296 0.04939048
0.04345098 0.03590105 0.03382174 0.03120917 0.02461127 0.02213363
0.01795666 0.01755989 0.01477137 0.01397352 0.01334286 0.01264867
0.01043573 0.00918058 0.00897625 0.00815701 0.00780193 0.00724554
0.00696956 0.00606637 0.00574502 0.00513423]
6.3.2人脸识别
#程序6-6
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.datasets import fetch_lfw_people
faces = fetch_lfw_people()
#对数据集乱序操作
random_indexes = np.random.permutation(len(faces.data))
X = faces.data[random_indexes]
#取出数据集前9个数据
X_faces = X[:9,:]
def plot_pictures(pictures):
#官网：https://matplotlib.org/api/pyplot_api.html
#nrows和ncols表示子图分布的行列数，即3*3个子图
#figsize表示画布的长宽，1000*1000
#subplot_kw表示字典类型(可选参数),把字典的关键字传递给add_subplot()来创建每个子图,其中xticks表示行坐标，yticks表示列坐标
#fig表示画布的长宽，即1000*1000
#axes表示子图的对象，使用axes.flat来访问子图
fig, axes = plt.subplots(3, 3, figsize=(10, 10),subplot_kw={'xticks':[], 'yticks':[]})
for i, ax in list(enumerate(axes.flat)):
#使用reshape来确定子图的像素分布2914 = 62*47
#参数cmap用于设置热图的Colormap
ax.imshow(pictures[i].reshape(62, 47), cmap='bone')
plt.show()
plot_pictures(X_faces)
#特征脸
pca = PCA()
pca.fit(X)
plot_pictures(pca.components_[:9,:])
运行结果：