OX01 什么是主成分分析
20个变量试图理解其中之间的关系：
=190种可能；
PCA就是用来简化多变量复杂关系的常用方法；比如，20个相关(冗余)的变量，经过PCA后，就转化成了5个线性无关的成分变量，并尽可能多的保住原始信息；
PCA就是将原有的特征进行重组，得到新的特征。PCA的主要思想是将n维特征映射到k维上，这k维是全新的正交特征也被称为主成分，是在原有n维特征的基础上重新构造出来的k维特征。
PCA作用: 降维，去除噪声，加快计算。
如第一主成分为：
它是k个观测变量的加权组合，对初始变量集的方差解释性大。
OX02 如何做主成分分析
主成分分析的数理基础：(构造两两指标间的相关系数矩阵，借助拉格朗日方程求出特征值与特征向量)---一行一观测，一列一指标参考附录1
主成分分析的步骤
“矩阵论”：
1采用标准化的方法对实际值进行无量纲化处理(有必要做)
2计算无量纲数据的协差阵/相关系数矩阵(标准化之后)
3计算相关系数矩阵的特征值并按大小排序
4根据累计贡献率，确定特征值个数；
5根据特征值对应的特征向量确定权重(特征向量即权重)
6PCA评价(可选项，pca还可以做评价)标准化方法：z-score 标准化
协方差，有单位，不方便进行比较；相关系数无单位，我们经常说的衡量两个变量之间相关性用的是皮尔逊相关系数。
相关系数与相关系数矩阵：有一个p*n维矩阵，任意两个变量之间可以计算一个相关系数，那么最终形成
*
阶相关系数方阵，我们就可以争对此方阵求特征值与特征向量。
得到特征向量(主成分)后，可以用来做线性回归(主成分之间线性无关)
“代数论”：
求主成分同样转化成了最优化问题，我们的目标是：选择出来的主成分尽可能的涵盖原始数据信息，那么就应该选择方差作为衡量指标，方差越大，含有信息越多。一般我们会使用方差(Variance)来定义样本之间的间距
对于如何找到一个轴，使得样本空间的所有点映射到这个轴的方差最大。就转化成了数学问题.图片来自 数据科学家联盟
其过程分为两步：样本归0
找到样本点映射后方差最大的单位向量w
OX03 如何求w
3.1 代数论
有了目标函数(方差)，求得w的梯度(一行一样本，一列一特征：对所有特征求偏导可得梯度向量)，应用梯度上升法就可以很快找到目标函数的最大值。参考附录2
3.1.1 求第一主成分
下面我们以二维数据为例，将其映射到一维，即求出一系列样本点的第一主成分。
#构建二维数据集
import numpy as np
import matplotlib.pyplot as plt
X = 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()将样本进行均值归0(demean)，即所有样本将去样本的均值。样本的分布没有改变，只是将坐标轴进行了移动。
def demean(X):
# axis=0按列计算均值，即每个属性的均值，1则是计算行的均值
return (X - np.mean(X, axis=0))
X_demean = demean(X)
# 注意看数据分布没变，但是坐标已经以原点为中心了
plt.scatter(X_demean[:, 0], X_demean[:, 1])
plt.show()
接下来就是对目标(方差)函数和梯度(导数)函数的定义。
首先定义目标函数：
def f(w,X):
return np.sum((X.dot(w)**2))/len(X)
然后根据梯度公式求梯度w：
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] += epsilon
w_2 = w.copy()
w_2[i] -= epsilon
res[i] = (f(w_1, X) - f(w_2, X)) / (2 * epsilon)
return res
在这里，有一处需要注意的地方：对于要求的轴，向量，实际上一个单位向量w，即要让向量w的模为1
因此我们使用np中的线性代数库linalg里面的norm()方法，求第二范数，也就是求算术平方根
def direction(w):
return w / np.linalg.norm(w)
然后就实现梯度上升代码的流程了：gradient_ascent()方法为梯度上升法的过程，在n_iters次循环中，每次获得新的w，如果新的w和旧的w对应的目标函数f值的差值满足精度要求，就表示找到了合适的w。此处要注意，对于每一次计算向量w之前都要进行单位化计算，使其模为1。
# 梯度上升法代码
def gradient_ascent(df, X, initial_w, eta, n_iters=1e4, epsilon=1e-8):
w = direction(initial_w)
cur_iter = 0
while cur_iter < n_iters:
gradient = df_math(w,X)
last_w = w
w = last_w + eta * gradient
w = direction(w) # 将w转换成单位向量
if (abs(f(w,X) - f(last_w,X)) < epsilon):
break
cur_iter += 1
return w
使用梯度上升法需要先设定一个初始值w和学习率
这里面还需要注意一点：线性回归中，通常将特征系数θ的值设为全部为0的向量。但在主成分分析中w的初始值不能为0！！！如果w初值为0，带入梯度求导的公式中，也是每次求梯度得到的都是没有任何方向的0。所以，要设置一组随机数：
initial_w = np.random.random(X.shape[1])
eta = 0.001
然后执行gradient_ascent()方法，计算梯度上升的结果，即要求的w(单位向量)
w = gradient_ascent(df_debug, X_demean, initial_w, eta)
# 输出
array([0.76567331, 0.64322965])
进行可视化，轴对应的方向，即将样本映射到该轴上的方差最大，这个轴就是一个主成分(第一主成分)：
plt.scatter(X_demean[:,0],X_demean[:,1])
plt.plot([0,w[0]*30],[0,w[1]*30], color='red')
plt.show()
3.1.2 求第n主成分
下一个主成分就是将数据在第一主成分上的分量去掉，再对新的数据求第一主成分。怎么理解这句话？
这样就得到第二主成分，依此循环就可以得到第n主成分。
求第二主成分的步骤：
首先，我们要将数据集在第一个主成分上的分量去掉。即用样本
减去样本在第一主成分分量上的映射，得到
。阵(的矩阵)与向量(的矩阵)点乘，得到的是一个的向量，这个元素表示中的每一个样本映射到向量方向上的模长
然后对向量reshape，使其成真正的的矩阵
再乘以，得到的矩阵，也就是把矩阵中每个样本的每个维度在方向上的分量。即
用样本减去样本在第一主成分上的映射
X_new = X - X.dot(w).reshape(-1,1) * w
# 看一下
plt.scatter(X_new[:,0], X_new[:,1])
plt.show()
所有数据把第一主成分的分量去掉，得到的就是与第一主成分的分量垂直的向量。
由于整个数据只有两个维度，第一个维度去掉之后，剩下的第二维度就是所有的内容了。因此样本在第二维度上的分布完全在一条直线上，不再有其他方向上分布的方差。
如何第二主成分的轴？就是将去掉第一主成分的分量的
带入到gradient_ascent()函数中，即可以得到新的轴。
w_new = gradient_ascent(df_math, X_new, initial_w, eta)
# 输出：
array([-0.64320916, 0.76569052])
求第n主成分
first_n_component()方法用于求X的前n个主成分。
首先进行数据归零操作，使用列表存储前n个主成分的方向向量。在for循环中每次求主成分时都设定一个随机的w(initial_w)，使用梯度上升法得到一个主成分后就去掉在这个主成分上的分量X_pca = X_pca - X_pca.dot(w).reshape(-1, 1) * w：
def first_n_component(n, X, eta=0.001, 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 = gradient_ascent(df_math, X_pca, initial_w, eta)
res.append(w)
X_pca = X_pca - X_pca.dot(w).reshape(-1, 1) * w
return res
3.2 矩阵论
相比于“代数论”求w，采用矩阵论求主成分会显得更加简洁。
from numpy.linalg import inv,det,eig,qr,svd
import numpy as np
import pandas as pd
X=pd.read_csv("C://Users//baihua//Desktop//principal.csv",encoding='utf-8') # 3列
# 定义Z-score函数
def autoNorm(dataSet):
norm_data=(dataSet-dataSet.mean(axis=0))/dataSet.std(axis=0)
return norm_data
# 计算特征值与单位特征向量
data=pd.DataFrame(autoNorm(X))
x=data.corr() # 计算相关系数矩阵
eig(x) # 计算特征值及单位特征向量
# 输出
(array([1.40749131, 0.99795365, 0.59455504]),
array([[-0.26177284, 0.93111755, 0.25395883],
[ 0.65783606, 0.36467863, -0.65898499],
[ 0.70620584, 0.00544109, 0.70798566]]))求出单位特征向量后，还需要计算主成份：
Zi= Xw，i=1,2,3...m
OX04 PCA的实现
在之前已经学习了如何求一个数据集的前n个主成分，但是数据集本身已经是n维的，并没有进行降维度。那么PCA是如何降维的呢？如何从高维数据向低维数据映射？
主成分分析的作用就是选出能使样本方差最大的维度，选择完维度之后，进入对数据降维的操作。将高维数据映射为低维数据。
假设经过主成分分析之后，左侧X还是数据样本，一个mn的矩阵，m个样本n个特征。根据主成分分析法求出了前k个主成分，得到矩阵，即有k个主成分向量，每个主成分的坐标系有n个维度(与转换前的维度相同)，形成一个k*n的矩阵。
如何将样本X从n维转换为k维呢？矩阵乘法：(m*n)
(n*k)=m*k
4.1 代码实现
import numpy as np
class PCA:
def __init__(self, n_components):
# 主成分的个数n
self.n_components = n_components
# 具体主成分
ponents_ = None
def fit(self, X, eta=0.001, n_iters=1e4):
'''均值归零'''
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, n_iters, epsilon=1e-8):
w = direction(initial_w)
cur_iter = 0
while cur_iter < n_iters:
gradient = df(w, X)
last_w = w
w = w + eta * gradient
w = direction(w)
if(abs(f(w, X) - f(last_w, X)) < epsilon):
break
cur_iter += 1
return w
# 过程如下：
# 归0操作
X_pca = demean(X)
# 初始化空矩阵，行为n个主成分，列为样本列数
ponents_ = np.empty(shape=(self.n_components, X.shape[1]))
# 循环执行每一个主成分
for i in range(self.n_components):
# 每一次初始化一个方向向量w
initial_w = np.random.random(X_pca.shape[1])
# 使用梯度上升法，得到此时的X_PCA所对应的第一主成分w
w = first_component(X_pca, initial_w, eta, n_iters)
# 存储起来
ponents_[i:] = w
# X_pca减去样本在w上的所有分量，形成一个新的X_pca，以便进行下一次循环
X_pca = X_pca - X_pca.dot(w).reshape(-1, 1) * w
return self
# 将X数据集映射到各个主成分分量中
def transform(self, X):
assert X.shape[1] == ponents_.shape[1]
return X.dot(ponents_.T)
def inverse_transform(self, X):
return X.dot(ponents_)
4.2 sklearn中的PCA
pca 对象的属性
components_ ：返回具有最大方差的成分。
explained_variance_ratio_：返回 所保留的n个成分各自的方差百分比。
n_components_：返回所保留的成分个数n。
pca对象的方法fit(X,y=None)
fit()可以说是scikit-learn中通用的方法，每个需要训练的算法都会有fit()方法，它其实就是算法中的“训练”这一步骤。
因为PCA是无监督学习算法，此处y自然等于None。fit_transform(X)
用X来训练PCA模型，同时返回降维后的数据。inverse_transform()
将降维后的数据转换成原始数据具体参考：链接
# sklearn 中的pca
from sklearn.decomposition import PCA
pca = PCA() # 初始实例化对象
pca.fit(x) # 训练模型
print("返回模型的各个特征向量",ponents_) #返回模型的各个特征向量(原始数据)
print("返回各个成分各自的方差百分比(贡献率)",pca.explained_variance_ratio_ )# 返回各个成分各自的方差百分比(贡献率)
# 根据以上贡献率，选取前n个主成份，累积贡献率达到95%左右。然后重新建立模型
pca = PCA(3)#这里假设前3个主成份已经足够代表数据集的大部分信息。
pca.fit(data)
low_d = pca.transform(data) #newX=pca.fit_transform(X)，newX就是降维后的数据
print("降维",low_d)
inv_d=pca.inverse_transform(data) #将降维后的数据转换成原始数据
print("inverse",inv_d)
4.3 降维算法pca的应用
我们对手写数据集digits使用主成分分析法进行降维，再用kNN算法进行分类，观察前后的结果有何不同。
import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
digits = datasets.load_digits()
X = digits.data
y = digits.target
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
接下来直接使用knn进行分类预测，
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train, y_train)
knn_clf.score(X_test, y_test)
0.9866666666666667
在使用pca后，继续使用knn
#下面用PCA算法对数据进行降维：
from sklearn.decomposition import PCA
pca = PCA(n_components=2) # 预先设置主成分有2个
pca.fit(X_train)
X_train_reduction = pca.transform(X_train) # 训练数据集降维结果 。X都需要降维，y不需要
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)
0.6066666666666667显然主成分个数是个超参数，那么如何设置主成分个数n_components了？
print(pca.explained_variance_ratio_) # pca.explained_variance_ratio_(解释方差比例)，我们可以使用这个指标找到某个数据集保持多少的精度：
[0.14566817 0.13735469]
对于pca = PCA(n_components=2) 来说，0.14566817表示第一个轴能够解释14.56%数据的方差；0.13735469表示第二个轴能够解释13.73%数据的方差。对于这两个维度来说，[ 0.14566817, 0.13735469]涵盖了原数据的总方差的28%左右的信息，剩下72%的方差信息就丢失了，显然丢失的信息过多。
# 横轴是是样本X的i个特征数，纵轴是前i个轴解释方差比例的和
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()
随着n_components 轴的增加，数据的方差会变大，直到n_components大到一定值后，解释比例保持不变，这个n_components就是我们想找的最优主成分！实际sklearn使用pca = PCA(0.95) 就可以保证传入95%以上的信息。
pca = PCA(0.95)
pca.fit(X_train)
# 输出：
PCA(copy=True, iterated_power='auto', n_components=0.95, random_state=None,
svd_solver='auto', tol=0.0, whiten=False)
# 查看主成分个数
pca.n_components_
28
然后用这种pca去重新使用kNN做训练，得到的结果较好。
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)
# 输出
0.97999999999999998降维的好处缩减了计算时间
数据降维还有一个作用是可视化，降到2维数据之后：
#数据降维还有一个作用是可视化，降到2维数据之后：
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)
plt.show()
OX05 PCA的去噪
5.1 数值去噪
前面提到，PCA的作用大致有三个：降低维度(将X从高纬度映射到低维度)、加快运算(维度约简后运算量降低)；这两个好理解，那么去噪又是怎么回事了？
PCA通过选取主成分将原有数据映射到低维数据再映射回高维数据的方式进行一定程度的降噪。
import numpy as np
import matplotlib.pyplot as plt
X = np.empty((100, 2))
X[:,0] = np.random.uniform(0., 100., size=100)
X[:,1] = 0.75 * X[:,0] + 3. + np.random.normal(0, 5, size=100) # 线性方程加一个正态分布的随机扰动项
plt.scatter(X[:,0], X[:,1])
plt.show()
# np.random.normal(0, 5, size=100) 是我们添加的噪音:随机数的均值、方差以及输出的size(可数可(2,2))
我们降噪，这就需要使用PCA中的一个方法：X_ori=pca.inverse_transform(X_pca)，将降维后的数据转换成与维度相同数据。要注意，还原后的数据，不等同于原数据！
这是因为在使用PCA降维时，已经丢失了部分的信息(忽略了解释方差比例)。因此在还原时，只能保证维度相同。会尽最大可能返回原始空间，但不会跟原来的数据一样。
from sklearn.decomposition import PCA
pca = PCA(n_components=1)
pca.fit(X)
X_reduction = pca.transform(X)
X_restore = pca.inverse_transform(X_reduction)
plt.scatter(X_restore[:,0], X_restore[:,1])
plt.show()
transform降维成一维数据，再inverse_transform返回成二维数据。此时数据成为了一条直线。这个过程可以理解为将原有的噪音去除了。
当然，我们丢失的信息也不全都是噪音。我们可以将PCA的过程定义为：降低了维度，丢失了信息，也去除了一部分噪音。
去除噪音，不仅仅去除的是数字噪音，图片中的噪音也可以通过PCA方法进行去噪。
5.2 手写数字去噪参考附录4
from sklearn import datasets
digits = datasets.load_digits()
X = digits.data
y = digits.target
# 添加一个正态分布的噪音矩阵
noisy_digits = X + np.random.normal(0, 4, size=X.shape)
# 绘制噪音数据：从y==0数字中取10个，进行10次循环；
# 依次从y==num再取出10个，将其与原来的样本拼在一起
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.shape
# 输出：
(100, 64) # 即含有100个数字(0～9各10个)，每个数字有64维
# 手写数字 ，用的是包里面的数据；工程中如何构建满足pca模型的原始图库，这个很耗费时间？！
def 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 = PCA(0.5).fit(noisy_digits)
pca.n_components_
# 输出：12，即原始数据保存了50%的信息，需要12维
# 进行降维、还原过程
components = pca.transform(example_digits)
filtered_digits = pca.inverse_transform(components)
plot_digits(filtered_digits)前后对比，明亮很多
附录 (获权)
目前知乎数据分析专栏，已经将近1600小伙伴关注了，点击下方链接关注专栏：数据分析
-------------------------正文结束--------------------------
#划重点#同步一个重要的事情，小编我近期开通了微信公号，
微信搜索 求知鸟，即可关注。
公号里有一篇自我介绍，欢迎来撩。
开通公号的原因：文章越来越多，很多号主开始转载我的文章，为了交流方便就开通了公号。
通过公号，可以很方便的获得代码，csv文件。
通过公号，可以更快的找到我。(现在加v，可以交流秋招经历哦)。
