PCA（主成分分析）

最新推荐文章于 2024-06-05 16:00:00 发布

negative_sun

最新推荐文章于 2024-06-05 16:00:00 发布

阅读量2.6k

点赞数 1

文章标签： python 机器学习

本文链接：https://blog.csdn.net/negative_sun/article/details/123582680

版权

PCA（主成分分析）

1.概述

通过计算，计算出主要的“特征值”，达到减少特征值个数的目的。（这里的特征值并不是原来的特征值，而是经过投影后的新特征值）
对数据的要求：PCA成立的假设前提是数据应该满足正交分布并且是高信噪比的（低通滤波特性），而这两个特点都正是高斯分布的特点，所以，PCA输入需要数据满足高斯分布。或者说，高斯分布是PCA适用的充分条件，但不是必要条件
优点：理论依据充分
缺点：由于降维过程是纯数学分析，可能丢失掉数据的实际代表含义
作用：分类、对指标进行建模时候的客观赋权（如最小风险贝叶斯中的风险函数）、构造新特征（因为得到了原来数据的投影）、异常数据监测（因为计算出了源数据的方差，方差大的是好数据，方差小的是噪声数据）

2.原理

采用了yale人脸数据集作为手动计算案例
步骤
1.计算均值向量u
2.计算中心化矩阵C
3.计算协方差矩阵Sigma
4.计算变换矩阵W
5.计算投影Y
6.计算测试集投影Z
7.遍历搜索匹配，预测Z的标签

ps.步骤6和步骤7是预测步骤，降维结果是Y

2.1 手动计算见：pca主成分分析.pdf

3.手动实现PCA算法

采用了Iris数据集和DryBean数据集
也进行了均衡标签/不均衡标签的划分

import pandas as pd
import numpy as np
import struct
import os
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.decomposition import PCA

#引入数据集
iris = datasets.load_iris()
x= iris.data
y= iris.target

#拆分训练集和测试集
split = 0.7
#trianset : testset = 7:3

np.random.seed(0) #固定随机结果
train_indices = np.random.choice(len(x),round(len(x) * split),replace=False) 
test_indices = np.array(list(set(range(len(x))) - set(train_indices)))

np.random.seed(0) #固定随机结果
train_indices = np.random.choice(len(x),round(len(x) * split),replace=False) 
test_indices = np.array(list(set(range(len(x))) - set(train_indices)))

train_indices = sorted(train_indices)
test_indices =sorted(test_indices)
train_x = x[train_indices]
test_x = x[test_indices]
train_y = y[train_indices]
test_y = y[test_indices]
# print(train_indices)
# print(test_indices)

#PCA核心计算过程
def MyPca(X,k):
    U = np.mean(X,axis=0)
    C = X - U
    n = len(X)
#     print(n)
    Cov = np.dot(C.T,C) / n
    eigval, eigvec = np.linalg.eig(Cov)
    indexes = np.argsort(eigval)[-k:]
    W = eigvec[:, indexes]
    Y = np.dot(C, W).T
#     print(pd.DataFrame(W).shape)
#     print(pd.DataFrame(Y).shape)
    return Y,W

#使用了矩阵乘法计算欧几里得距离，缩短计算时间
#计算原理见附件：矩阵间的欧氏距离.pdf
def test(test_x,Y,test_y,train_y,W):
    n = len(test_x)
    Z = test_x
    U = np.mean(Z, axis=0)
    C = Z - U
    ZZ = np.dot(C, W).T # k * n
    #求出ZZ的每一行与Y的每一行的欧式距离
    #ZZ的第一行表示第一个待测样本
    #对ZZ的第一行进行广播，与Y的每一行求欧式距离
#     print(ZZ.T.shape,Y.shape)
    dists = np.sqrt(-2 * np.dot(ZZ.T, Y) + np.repeat([np.sum(np.square(ZZ.T), axis=1)],repeats=Y.shape[1],axis=0).T + np.repeat(np.transpose([np.sum(np.square(Y), axis=0)]),repeats=ZZ.shape[1],axis=1).T)
    # print(dists.shape)
    # print(np.argmin(dists, axis=1))
    pre = np.argmin(dists, axis=1)
    right = 0
    test_y = pd.DataFrame(test_y)
    train_y = pd.DataFrame(train_y)
    for i in range(len(pre)):
#         print(i,pre[i],len(test_y),len(train_y))
#         print(test_y[i],train_y.iloc[1212])
        #ub_train_y.iloc[1212][0]
        if test_y.iloc[i][0] == train_y.iloc[pre[i]][0]:
            right=right+1
    return right/n

#抽取不均衡标签的数据集
from numpy import *
i_train_x = train_x
i_train_y = train_y
i_test_x = test_x
i_test_y = test_y
ui_train_x = concatenate((concatenate((x[:45],x[50:95]),axis=0),x[100:115]),axis=0) #axis=0表示竖向拼接
ui_train_y = concatenate((concatenate((y[:45],y[50:95]),axis=0),y[100:115]),axis=0)
ui_test_x = concatenate((concatenate((x[:5],x[50:55]),axis=0),x[100:135]),axis=0)
ui_test_y = concatenate((concatenate((y[:5],y[50:55]),axis=0),y[100:135]),axis=0)

#比较均衡标签和不均衡标签的效果
klist = []
i_score = []
ui_score = []
for k in range(1,5):  
    klist.append(k)
    Y,W = MyPca(i_train_x,k)
    ret = test(i_test_x,Y,i_test_y,i_train_y,W)
    i_score.append(ret)
    Y,W = MyPca(ui_train_x,k)
    ret = test(ui_test_x,Y,ui_test_y,ui_train_y,W)
    ui_score.append(ret)

plt.plot(klist, i_score, marker = 'o', label = 'banlanced Iris')
plt.plot(klist, ui_score, marker = '*', label = 'unbanlanced Iris')
plt.legend() #让图例生效
plt.xlabel('k-value')
plt.ylabel('accuracy-value')
plt.title(u'Iris map')
plt.show()

在这里插入图片描述

3.1 接下来验证自我实现的结果与调包的结果一样

取k=3，画出3D图来展示
发现计算结果仅有方向上的差别

from mpl_toolkits.mplot3d import Axes3D
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from matplotlib import cm
klist = []
i_score = []
ui_score = []
for k in range(3,4):  
    klist.append(k)
    Y,W = MyPca(i_train_x,k)
    fig = plt.figure()
    ax = fig.add_subplot(projection='3d')
    ax.scatter3D(Y[0,: ], Y[1,: ] ,Y[2,:] * -1)  # 三个数组对应三个维度（三个数组中的数一一对应）
    plt.xlabel('X')
    plt.ylabel('Y', rotation=38)  # y 轴名称旋转 38 度
    ax.set_zlabel('Z')
    
    pca = PCA(n_components=k)
    X_new = pca.fit_transform(i_train_x)
   
    fig = plt.figure()
    ax = fig.add_subplot(projection='3d')
    ax.scatter3D(X_new[:, 2], X_new[:, 1] * -1,X_new[:,0] * -1)  # 三个数组对应三个维度（三个数组中的数一一对应）
    # print(Y[0,:].flatten().shape,Y[1,:].flatten().shape,Y[2,:].flatten().shape)
    # surf = ax.plot_trisurf(Y[0,:].tolist(),Y[1,:].tolist(),Y[2,:].tolist(),cmap=cm.coolwarm,linewidth=0, antialiased=False)
    plt.xlabel('X')
    plt.ylabel('Y', rotation=38)  # y 轴名称旋转 38 度
    ax.set_zlabel('Z')  # 因为 plt 不能设置 z 轴坐标轴名称，所以这里只能用 ax 轴来设置（当然，x 轴和 y 轴的坐标轴名称也可以用 ax 设置）
    plt.savefig('3D-2.jpg', bbox_inches='tight', dpi=2400)  # 保存图片，如果不设置 bbox_inches='tight'，保存的图片有可能显示不全
    plt.show()

在这里插入图片描述

3.2 用DryBean数据集试一试，看看效果

import operator
from sklearn.preprocessing import StandardScaler  # 均值归一化
from sklearn.metrics import confusion_matrix  # 生成混淆矩阵
from sklearn.metrics import classification_report  # 分类报告

def openfile(filename):
    """
    打开数据集，进行数据处理
    :param filename:文件名
    :return:特征集数据、标签集数据
    """    
    # 打开excel
    sheet = pd.read_excel(filename,sheet_name='Dry_Beans_Dataset')
    
    data = sheet.iloc[:,:16].values
    target = sheet['Class'].values
#     print(data.shape)
#     print(target.shape)                                            
    return data, target, sheet.columns

def split_data_set(data_set, target_set, rate=0.7):
    """
    说明：分割数据集，默认数据集的30%是测试集

    :param data_set: 数据集
    :param target_set: 标签集
    :param rate: 测试集所占的比率
    :return: 返回训练集数据、训练集标签、测试集数据、测试集标签
    """
    # 计算训练集的数据个数
    train_size = len(data_set)
    # 随机获得数据的下标
    train_index = sorted(np.random.choice(train_size,round(train_size * rate), replace=False))
    test_index = sorted(np.array(list(set(range(train_size)) - set(train_index)))) #不用排序也行，强迫症，为了上面保持一致就排序了
    # 分割数据集（X表示数据，y表示标签）
    x_train = data_set.iloc[train_index,:] #因为这里的data_set和target_set变成DataFrame，而不是ndarray了，所以要用iloc访问
    x_test = data_set.iloc[test_index,:]
    y_train = target_set.iloc[train_index,:]
    y_test = target_set.iloc[test_index,:]
    return x_train, y_train, x_test, y_test

filename = r'D:\jjq\code\jupyterWorkSpace\datasets\DryBeanDataset\Dry_Bean_Dataset.xlsx'
o_bean_dataset = openfile(filename)
#每个类别的种子抽取400条数据，这个是每个类别的起始索引
step = 400
start_index = [0,1322,1844,3474,7020,8948,10975] #一共7类
# bean_dataset_x = pd.DataFrame(columns=o_bean_dataset[2])
# bean_dataset_y  =pd.DataFrame(columns=o_bean_dataset[2])
bean_dataset_x = pd.DataFrame(columns=range(16))
bean_dataset_y  =pd.DataFrame(columns=range(1))
bean_dataset_x.drop(bean_dataset_x.index,inplace=True)
bean_dataset_y.drop(bean_dataset_y.index,inplace=True)
for i in range(7):
    bean_dataset_x = pd.concat((bean_dataset_x, pd.DataFrame(o_bean_dataset[0][start_index[i]:(step+start_index[i])])),axis=0)
    bean_dataset_y = pd.concat((bean_dataset_y, pd.DataFrame(o_bean_dataset[1][start_index[i]:(step+start_index[i])])),axis=0)
# bean_dataset_y.to_excel("./123.xlsx")

#按照均衡和不均衡的方式，划分训练集和测试集

#均衡
b_train_x, b_train_y, b_test_x, b_test_y = split_data_set(bean_dataset_x,bean_dataset_y)
# print(b_train_x.shape,b_train_y.shape)
# print(b_test_x.shape,b_test_y.shape)
#不均衡
steps_train = [318,318,318,318,318,212,159]
steps_test = [68,91,136,136,136,136,136]
now = 0
#初始化不均衡数组
ub_train_x = pd.DataFrame(columns=range(16))
ub_test_x = pd.DataFrame(columns=range(16))
ub_train_y = pd.DataFrame(columns=range(1))
ub_test_y = pd.DataFrame(columns=range(1))
#保证添加数据之前数组为空
ub_train_x.drop(ub_train_x.index,inplace=True)
ub_test_x.drop(ub_test_x.index,inplace=True)
ub_train_y.drop(ub_train_y.index,inplace=True)
ub_test_y.drop(ub_test_y.index,inplace=True)

#开始添加数据
for i in range(7):
    ub_train_x = pd.concat((ub_train_x, bean_dataset_x[now:(now+steps_train[i])]),axis=0) 
    ub_train_y = pd.concat((ub_train_y, bean_dataset_y[now:(now+steps_train[i])]),axis=0)
    now = now+steps_train[i]
    ub_test_x = pd.concat((ub_test_x, bean_dataset_x[now:(now+steps_test[i])]),axis=0)
    ub_test_y = pd.concat((ub_test_y, bean_dataset_y[now:(now+steps_test[i])]),axis=0)
    now = now+steps_test[i]

b_score = []
ub_score = []
klist = []
for k in range(1,10):  
    klist.append(k)
    Y,W = MyPca(b_train_x,k)
    ret = test(b_test_x,Y,b_test_y,b_train_y,W)
    b_score.append(ret)
    Y,W = MyPca(ub_train_x,k)
    ret = test(ub_test_x,Y,ub_test_y,ub_train_y,W)
    ub_score.append(ret)

plt.plot(klist, b_score, marker = 'o', label = 'banlanced DryBean')
plt.plot(klist, ub_score,marker = '*', label = 'unbalanced DryBean')
plt.legend() #让图例生效
plt.xlabel('k-value')
plt.ylabel('accuracy-value')
plt.title(u'DryBean map')
plt.show()

在这里插入图片描述

#这里发现，Iris数据集和DryBean数据集都是均衡标签的数据效果更好
#接下来寻找原因

4.寻找效果不同的原因

4.1首先查看两个数据集的分布

发现分布一样

from sklearn.preprocessing import StandardScaler  # 均值归一化
scaler = StandardScaler()
b_train_x = pd.DataFrame(scaler.fit_transform(b_train_x))

plt.figure(figsize=(20,30))
p = 1
k=14
for i in range(k):
    a = plt.subplot(k,2,p)
    a.plot(range(b_train_x.shape[0]),b_train_x.iloc[:,i].tolist(),label="balance "+str(i),color='red')
    p = p+1
    plt.legend() #让图例生效
    b = plt.subplot(k,2,p)
    b.plot(range(ub_train_x.shape[0]),ub_train_x.iloc[:,i].tolist(),label="unbalance "+str(i))
    p = p+1
    plt.legend() #让图例生效

plt.xlabel('index')
plt.ylabel('value')
plt.title(u'DryBean map')
plt.show()

在这里插入图片描述

# plt.subplots_adjust(hspace=100, wspace=10)
plt.figure(figsize=(15,40))
p = 1
for i in range(14):
    a = plt.subplot(14,2,p)
    a.set_title(b_train_x.columns.values.tolist()[i])
    a.hist(x = b_train_x.iloc[:,i].tolist(), # 指定绘图数据
          bins = 50, # 指定直方图中条块的个数
           color='red'
    )
    p=p+1
    a = plt.subplot(14,2,p)
    a.set_title(ub_train_x.columns.values.tolist()[i])
    a.hist(x = ub_train_x.iloc[:,i].tolist(), # 指定绘图数据
          bins = 50, # 指定直方图中条块的个数
    )
    p=p+1
plt.show()

在这里插入图片描述

plt.figure(figsize=(15,5))
b = plt.subplot(1,2,1)
b.hist(x = b_test_y.iloc[:,:].values[:,0].tolist(), # 指定绘图数据
          bins = 50, # 指定直方图中条块的个数
           color="red"
      )
c = plt.subplot(1,2,2)
c.hist(x = ub_test_y.iloc[:,:].values[:,0].tolist(), # 指定绘图数据
          bins = 50, # 指定直方图中条块的个数
      )
plt.show()

在这里插入图片描述

4.2根据理论，pca需要符合高斯分布的数据

使用函数，计算两个数据集的“正态分数”，发现效果更好的数据确实更接近正态分布一些（强行解释）

#DryBean数据集
from scipy.stats import normaltest
k2, pp = normaltest(b_train_x)
k3,ppp = normaltest(ub_train_x)
plt.plot(range(len(pp)),pp,label="balance")
plt.plot(range(len(ppp)),ppp,label="ubbalance")
plt.legend()
plt.show()
# print(k2-k3,pp-ppp)
#统计量p越接近0就越表明数据和标准正态分布拟合的越好

在这里插入图片描述

#Iris数据集
from scipy.stats import normaltest
k2, pp = normaltest(i_train_x)
k3,ppp = normaltest(ui_train_x)
plt.plot(range(len(pp)),pp,label="balance")
plt.plot(range(len(ppp)),ppp,label="ubbalance")
plt.legend()
plt.show()

在这里插入图片描述

5.sklearn中PCA

原型：sklearn.decomposition.PCA(n_components=None, copy=True, whiten=False)

5.1属性

components_ ：返回具有最大方差的成分。
explained_variance_ratio_：返回所保留的n个成分各自的方差百分比。
n_components_：返回所保留的成分个数n。
mean_：均值
noise_variance_：噪声比

5.2 方法

fit(X,y=None)：fit(X)，表示用数据X来训练PCA模型。函数返回值：调用fit方法的对象本身。比如pca.fit(X)，表示用X对pca这个对象进行训练
fit_transform(X)：用X来训练PCA模型，同时返回降维后的数据
inverse_transform()：将降维后的数据转换成原始数据，X=pca.inverse_transform(newX)
transform(X)：将数据X转换成降维后的数据。当模型训练好后，对于新输入的数据，都可以用transform方法来降维

negative_sun

关注

1
点赞
踩
4

收藏

觉得还不错? 一键收藏
0
评论
PCA（主成分分析）

PCA（主成分分析）1.概述通过计算，计算出主要的“特征值”，达到减少特征值个数的目的。（这里的特征值并不是原来的特征值，而是经过投影后的新特征值）对数据的要求：PCA成立的假设前提是数据应该满足正交分布并且是高信噪比的（低通滤波特性），而这两个特点都正是高斯分布的特点，所以，PCA输入需要数据满足高斯分布。或者说，高斯分布是PCA适用的充分条件，但不是必要条件优点：理论依据充分缺点：由于降维过程是纯数学分析，可能丢失掉数据的实际代表含义作用：分类、对指标进行建模时候的客观赋权（如最小风险贝叶
复制链接

扫一扫