1 基本概念
PCA即主成分分析技术。主成分分析也称主分量分析，旨在利用降维的思想，把多指标转化为少数几个综合指标。在统计学中，主成分分析PCA是一种简化数据集的技术。它是一个线性变换。这个变换把数据变换到一个新的坐标系统中，使得任何数据投影的第一大方差在第一个坐标(称为第一主成分)上，第二大方差在第二个坐标(第二主成分)上，依次类推。主成分分析经常用于减少数据集的维数，同时保持数据集的对方差贡献最大的特征。这是通过保留低阶主成分，忽略高阶主成分做到的。这样低阶成分往往能够保留住数据的最重要方面。但是，这也不是一定的，要视具体应用而定。2 原理与数学推导
1.主成分分析使用的是梯度上升法。
特点
原理：
数学推导：
3 自己实现PCA算法
3.1 利用梯度上升法模拟实现PCA
1.数据的模拟准备
import numpy as npimport matplotlib.pyplot as pltX = np.empty((100, 2))X[:,0] = np.random.uniform(0., 100., size=100)X[:,1] = 0.75 * X[:,0] + 3. + np.random.normal(0, 10., size=100)plt.scatter(X[:,0], X[:,1])plt.show()
2.demean操作（对每一特征，按X的列进行均值归0操作）
def demean(X):return X - np.mean(X, axis=0)X_demean = demean(X)
3.使用梯度上升法
def f(w, X):return np.sum((X.dot(w)**2)) / len(X)def df_math(w, X):return X.T.dot(X.dot(w)) * 2. / len(X)def df_debug(w, X, epsilon=0.0001):res = np.empty(len(w))for i in range(len(w)):w_1 = w.copy()w_1[i] += epsilonw_2 = w.copy()w_2[i] -= epsilonres[i] = (f(w_1, X) - f(w_2, X)) / (2 * epsilon)return resdef direction(w):return w / np.linalg.norm(w) # 对向量进行单位化def gradient_ascent(df, X, initial_w, eta, n_iters = 1e4, epsilon=1e-8):w = direction(initial_w) cur_iter = 0while cur_iter < n_iters:gradient = df(w, X)last_w = ww = w + eta * gradientw = direction(w) # 注意1：每次求一个单位方向if(abs(f(w, X) - f(last_w, X)) < epsilon):breakcur_iter += 1return winitial_w =np.random.random(X.shape[1]) # 注意2 不能从0向量开始w = gradient_ascent(df_math,X_demean,initial_w,eta)plt.scatter(X_demean[:,0],X_demean[:,1])plt.plot([0,w[0]*30],[0,w[1]*30],color = 'r')plt.show()
注意3：不能使用StandardScaler标准化数据
4.上面求出的是第一主成分分量，若求其他主成分分量需减掉第一主成分分量再放入其中进行求解。
#其他主成分分量X2 = np.empty(X.shape)for i in range(len(X)):X2[i] = X[i] - X[i].dot(w)*w #去掉第一主成分分量#或可以表示为X2 = X - X.dot(w).reshape(-1,1)*w
5.求前n个主成分分量
def first_n_components(n,X,eta=0.01,n_iters=1e4,epsilon=1e-8):X_pca = X.copy()X_pca = demean(X_pca)res = [] # 存放主成分for i in range(n):initial_w = np.random.random(X_pca.shape[1])w = first_component(df,X_pca,initial_w,eta)res.append(w)X_pca = X_pca - X_pca.dot(w).reshape(-1,1)*wreturn resfirst_n_components(2,X) #调用
5.自己实现的PCA并将其封装
import numpy as npclass PCA:def __init__(self, n_components):"""初始化PCA"""assert n_components >= 1, "n_components must be valid"self.n_components = ponents_ = Nonedef fit(self, X, eta=0.01, n_iters=1e4):"""获得数据集X的前n个主成分"""assert self.n_components <= X.shape[1], \"n_components must not be greater than the feature number of X"def demean(X):return X - np.mean(X, axis=0)def f(w, X):return np.sum((X.dot(w) ** 2)) / len(X)def df(w, X):return X.T.dot(X.dot(w)) * 2. / len(X)def direction(w):return w / np.linalg.norm(w)def first_component(X, initial_w, eta=0.01, n_iters=1e4, epsilon=1e-8):w = direction(initial_w)cur_iter = 0while cur_iter < n_iters:gradient = df(w, X)last_w = ww = w + eta * gradientw = direction(w)if (abs(f(w, X) - f(last_w, X)) < epsilon):breakcur_iter += 1return wX_pca = demean(X)ponents_ = np.empty(shape=(self.n_components, X.shape[1]))for i in range(self.n_components):initial_w = np.random.random(X_pca.shape[1])w = first_component(X_pca, initial_w, eta, n_iters)ponents_[i,:] = wX_pca = X_pca - X_pca.dot(w).reshape(-1, 1) * wreturn selfdef transform(self, X):"""将给定的X，映射到各个主成分分量中"""assert X.shape[1] == ponents_.shape[1]return X.dot(ponents_.T)def inverse_transform(self, X):"""将给定的X，反向映射回原来的特征空间"""assert X.shape[1] == ponents_.shape[0]return X.dot(ponents_)def __repr__(self):return "PCA(n_components=%d)" % self.n_components
数据的降维
import numpy as npimport matplotlib.pyplot as pltX = np.empty((100, 2))X[:,0] = np.random.uniform(0., 100., size=100)X[:,1] = 0.75 * X[:,0] + 3. + np.random.normal(0, 10., size=100)from playML.PCA import PCApca = PCA(n_components=2)pca.fit(X)ponents_pca = PCA(n_components=1)pca.fit(X)X_reduction = pca.transform(X)X_reduction.shapeX_restore = pca.inverse_transform(X_reduction)X_restore.shapeplt.scatter(X[:,0], X[:,1], color='b', alpha=0.5)plt.scatter(X_restore[:,0], X_restore[:,1], color='r', alpha=0.5)plt.show()
n_components代表要求的几个主成分。X_reduction表示降维后的矩阵。X_restore表示降维后又恢复的矩阵。
4 sklearn中的PCA算法
1.利用pca实现降维
from sklearn.decomposition import PCApca = PCA(n_components=1)pca.fit(X)ponents_X_reduction = pca.transform(X)X_restore = pca.inverse_transform(X_reduction)
2.以具体的手写数字识别来用PCA进行降维之后，然后用KNN分类
import numpy as npimport matplotlib.pyplot as pltfrom sklearn import datasetsdigits = datasets.load_digits()X = digits.datay = digits.targetfrom sklearn.model_selection import train_test_splitX_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)from sklearn.decomposition import PCApca = PCA(n_components=2)pca.fit(X_train)X_train_reduction = pca.transform(X_train)X_test_reduction = pca.transform(X_test)knn_clf = KNeighborsClassifier()knn_clf.fit(X_train_reduction, y_train)knn_clf.score(X_test_reduction, y_test)
3.取的主成分分量的个数来解释代表主要的特征（pca.explained_variance_ratio_主成分所解释的方差）即重要程度所需的主成分数
from sklearn.decomposition import PCApca = PCA(n_components=X_train.shape[1])pca.fit(X_train)plt.plot([i for i in range(X_train.shape[1])], [np.sum(pca.explained_variance_ratio_[:i+1]) for i in range(X_train.shape[1])])plt.show()
或使用以下：
pca = PCA(0.95)pca.fit(X_train)pca.n_components_X_train_reduction = pca.transform(X_train)X_test_reduction = pca.transform(X_test)knn_clf = KNeighborsClassifier()knn_clf.fit(X_train_reduction, y_train)knn_clf.score(X_test_reduction, y_test)
pca.n_components_为对应比例是0.95所有的主成分分量的个数
4.对数据进行降维之后的数据的可视化
%%time pca = PCA(n_components=2)pca.fit(X)X_reduction = pca.transform(X)for i in range(10):plt.scatter(X_reduction[y==i,0], X_reduction[y==i,1], alpha=0.8)
5 PCA的其他应用
1. 利用PCA进行降噪
from sklearn import datasetsdigits = datasets.load_digits()X = digits.datay = digits.targetnoisy_digits = X + np.random.normal(0, 4, size=X.shape)example_digits = noisy_digits[y==0,:][:10]for num in range(1,10):example_digits = np.vstack([example_digits, noisy_digits[y==num,:][:10]])example_digits.shapedef plot_digits(data):fig, axes = plt.subplots(10, 10, figsize=(10, 10),subplot_kw={'xticks':[], 'yticks':[]},gridspec_kw=dict(hspace=0.1, wspace=0.1)) for i, ax in enumerate(axes.flat):ax.imshow(data[i].reshape(8, 8),cmap='binary', interpolation='nearest',clim=(0, 16))plt.show()plot_digits(example_digits)pca = PCA(0.5).fit(noisy_digits)pca.n_components_components = pca.transform(example_digits)filtered_digits = pca.inverse_transform(components)plot_digits(filtered_digits)components = pca.transform(example_digits)filtered_digits = pca.inverse_transform(components)plot_digits(filtered_digits)
利用PCA的正反变换，去除存在的噪音
2.特征脸
import numpy as npimport matplotlib.pyplot as pltfrom sklearn.datasets import fetch_lfw_peoplefaces = fetch_lfw_people()faces.keys()faces.data.shapefaces.target_namesfaces.images.shaperandom_indexes = np.random.permutation(len(faces.data))X = faces.data[random_indexes]example_faces = X[:36,:]example_faces.shapedef plot_faces(faces):fig, axes = plt.subplots(6, 6, figsize=(10, 10),subplot_kw={'xticks':[], 'yticks':[]},gridspec_kw=dict(hspace=0.1, wspace=0.1)) for i, ax in enumerate(axes.flat):ax.imshow(faces[i].reshape(62, 47), cmap='bone')plt.show()plot_faces(example_faces)# 特征脸from sklearn.decomposition import PCA pca = PCA(svd_solver='randomized')pca.fit(X)ponents_.shapeplot_faces(ponents_[:36,:])faces2 = fetch_lfw_people(min_faces_per_person=60)faces2.data.shapefaces2.target_nameslen(faces2.target_names)len(faces2.target_names)
