主成分分析（理论及Matlab实现）

最新推荐文章于 2024-05-10 13:10:06 发布

Xcee

最新推荐文章于 2024-05-10 13:10:06 发布

阅读量8.1k

点赞数 10

分类专栏：算法学习和使用文章标签： matlab python 算法数据分析

本文链接：https://blog.csdn.net/xuce_sir/article/details/125960933

版权

算法学习和使用专栏收录该内容

1 篇文章 1 订阅

订阅专栏

应用场景

主成分分析（Principal Component Analysis，PCA），主要用于对特征变量进行降维，在面对特征变量过多的时候使用。由于变量之间具有一定的相关关系，此时可以解释为这两个变量反映的信息有一定的重叠。针对这些重复的变量（关系紧密的变量）可以进行删减融合，建立尽可能少的新变量，使得这些新变量是两两不相关的，而且这些新变量在反映的信息方面尽可能保持原有的信息。

将原来变量重新组合成一组新的互相无关的几个综合变量，同时根据实际需要从中可以取出几个较少的综合变量尽可能多地反映原来变量的信息的统计方法。

个人使用较多的场景是对高维数据预处理和回归特征的提取。

原理

主要思想

PCA的主要思想是将n维特征映射到k维上，这k维是全新的正交特征也被称为主成分，是在原有n维特征的基础上重新构造出来的k维特征。PCA的工作就是从原始的空间中顺序地找一组相互正交的坐标轴，新的坐标轴的选择与数据本身是密切相关的。其中，第一个新坐标轴选择是原始数据中方差最大的方向，第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的，第三个轴是与第1,2个轴正交的平面中方差最大的。依次类推，可以得到n个这样的坐标轴。通过这种方式获得的新的坐标轴，我们发现，大部分方差都包含在前面k个坐标轴中，后面的坐标轴所含的方差几乎为0。于是，我们可以忽略余下的坐标轴，只保留前面k个含有绝大部分方差的坐标轴。事实上，这相当于只保留包含绝大部分方差的维度特征，而忽略包含方差几乎为0的特征维度，实现对数据特征的降维处理。

通过计算数据矩阵的协方差矩阵，然后得到协方差矩阵的特征值特征向量，选择特征值最大(即方差最大)的k个特征所对应的特征向量组成的矩阵。这样就可以将数据矩阵转换到新的空间当中，实现数据特征的降维。

计算步骤

假设有 n个样本，每个样本有p 个指标，则可构成大小为 n × p 的样本矩阵 X :

$\mathrm{X}=\left[\begin{array}{cccc} \mathrm{x}_{11} & \mathrm{x}_{12} & \cdots & \mathrm{x}_{1 \mathrm{p}} \\ \mathrm{x}_{21} & \mathrm{x}_{22} & \cdots & \mathrm{x}_{2 \mathrm{p}} \\ \vdots & \vdots & \ddots & \vdots \\ \mathrm{x}_{\mathrm{n} 1} & \mathrm{x}_{\mathrm{n} 2} & \cdots & \mathrm{x}_{\mathrm{np}} \end{array}\right]$

1. 首先对数据进行标准化处理：

$X_{ij} = \frac{x_{ij}-\bar{x} _{j}}{S_{j}}$

其中： $\bar{x} _{j}$ 为每列的均值； $S_{j}$ 为每列的标准差，公示如下：

$\mathrm{S}_{\mathrm{j}}=\sqrt{\frac{\sum_{\mathrm{i}=1}^{\mathrm{n}}\left(\mathrm{x}_{\mathrm{ij}}-\overline{\mathrm{x}}_{\mathrm{j}}\right)^{2}}{\mathrm{n}-1}}$

将X更新为标准化后的矩阵：

$\mathrm{X}=\left[\begin{array}{cccc} \mathrm{X}_{11} & \mathrm{X}_{12} & \cdots & \mathrm{X}_{1 \mathrm{p}} \\ \mathrm{X}_{21} & \mathrm{X}_{22} & \cdots & \mathrm{X}_{2 \mathrm{p}} \\ \vdots & \vdots & \ddots & \vdots \\ \mathrm{X}_{\mathrm{n} 1} & \mathrm{X}_{\mathrm{n} 2} & \cdots & \mathrm{X}_{\mathrm{np}} \end{array}\right]$

2. 计算标准化样本的协方差矩阵

$\mathrm{R}=\left[\begin{array}{cccc} \mathrm{r}_{11} & \mathrm{r}_{12} & \cdots & \mathrm{r}_{1 \mathrm{p}} \\ \mathrm{r}_{21} & \mathrm{r}_{22} & \cdots & \mathrm{r}_{2 \mathrm{p}} \\ \vdots & \vdots & \ddots & \vdots \\ \mathrm{r}_{\mathrm{p} 1} & \mathrm{r}_{\mathrm{p} 2} & \cdots & \mathrm{r}_{\mathrm{pp}} \end{array}\right]$

协方差则一般用来刻画两个随机变量的相似程度，协方差计算公式为：

$\mathrm{r}_{\mathrm{ij}}=\frac{1}{\mathrm{n}-1} \sum_{\mathrm{k}=1}^{\mathrm{n}}\left(\mathrm{X}_{\mathrm{ki}}-\overline{\mathrm{X}}_{\mathrm{i}}\right)\left(\mathrm{X}_{\mathrm{ki}}-\overline{\mathrm{X}}_{\mathrm{j}}\right)=\frac{1}{\mathrm{n}-1} \sum_{\mathrm{k}=1}^{\mathrm{n}} \mathrm{X}_{\mathrm{ki}} \mathrm{X}_{\mathrm{kj}}$

3. 计算协方差矩阵的特征值和特征向量

协方差矩阵R具有n个特征值 $\lambda _{1}$ ， $\lambda _{2}$ ，…， $\lambda _{p}$ ；及n个对应的特征向量，求解特征向量：

$\mathrm{a}_{1}=\left[\begin{array}{c} \mathrm{a}_{11} \\ \mathrm{a}_{21} \\ \vdots \\ \mathrm{a}_{\mathrm{p} 1} \end{array}\right], \mathrm{a}_{2}=\left[\begin{array}{c} \mathrm{a}_{12} \\ \mathrm{a}_{22} \\ \vdots \\ \mathrm{a}_{\mathrm{p} 2} \end{array}\right], \cdots, \mathrm{a}_{\mathrm{p}}=\left[\begin{array}{c} \mathrm{a}_{1 \mathrm{p}} \\ \mathrm{a}_{2 \mathrm{p}} \\ \vdots \\ \mathrm{a}_{\mathrm{pp}} \end{array}\right]$

4.计算主成分贡献率和累计贡献率

每个特征的贡献率为其特征值所占的比值，即：

$rate= \frac{\lambda _{i}}{\sum_{k=1}^{p}\lambda_{k}}$

累计贡献率为：

$total\_rate=\frac{\sum_{\mathrm{k}=1}^{\mathrm{i}} \lambda_{\mathrm{k}}}{\sum_{\mathrm{k}=1}^{\mathrm{p}} \lambda_{\mathrm{k}}}$

5. 找出主成分

统计累计贡献率，只保留贡献率较高的前几个变量值，一般贡献率取85%，写出对应的前m个主成分，其中m肯定是一个小于p的值。

第i个主成分：

$\mathrm{F}_{\mathrm{i}}=\mathrm{a}_{1 \mathrm{i}} \mathrm{X}_{1}+\mathrm{a}_{2 \mathrm{i}} \mathrm{X}_{2}+\cdots+\mathrm{a}_{\mathrm{pi}} \mathrm{X}_{\mathrm{p}}$

实例分析

对某磨机数据进行处理后得到如下的数据：

时间	给矿量x1	返砂量x2	磨机总功率x3	给矿水流量x4	排矿水流量x5	冲筛水流量x6	返砂比
00:00	960.96	90.34	6763.89	217.63	513.99	54.59	0.088
00:40	960.27	85.29	6623.43	214.22	340.37	54.60	0.089
01:20	959.87	85.74	6595.99	213.99	332.08	54.64	0.090
02:00	960.37	86.96	6622.91	214.07	343.91	54.70	0.097
02:40	960.24	94.04	6654.79	214.14	339.09	54.790	0.105
03:20	955.26	101.05	6768.13	214.11	369.28	54.84	0.119
04:00	963.32	113.91	6859.26	217.57	304.72	54.93	0.145
04:40	961.68	140.09	7129.90	219.97	298.97	54.91	0.137
05:20	960.08	132.70	7409.50	225.64	306.73	54.86	0.122
06:00	961.26	117.86	7254.85	223.43	300.85	54.88	0.111
06:40	960.28	107.45	7052.22	220.81	335.64	54.89	0.112
07:20	959.52	107.78	6974.43	217.88	331.53	54.96	0.097
08:00	961.40	93.65	6899.92	218.64	320.95	54.94	0.087
08:40	960.10	83.74	6670.77	215.38	304.63	55.07	0.083

处理后的数据为每隔40分钟某工厂半自磨机（SAG）运行参数整合，数据量n为2907条，上述表格只展示部分数据，特征向量有给矿量、返砂量、磨机总功率、给矿水流量、排矿水流量、冲筛水流量6个指标，现对这六个指标进行主成分分析，对指标进行降维。

Matlab代码

主程序：

Data = xlsread('data.xlsx');
data = Data(:,3:end-1);    %数据导入
data2 = PCA(data,0.85);    %运行主成分分析子程序

导入数据data，为已处理好的数据，PCA为子程序

主成分分析子程序：

function Data_var = PCA(data,n)
%主成分分析，输入数组和累计贡献率，得到处理后的变量数组。
%输入数组data;行为数据条数，列为指标
%输入n,为可接受累计贡献率范围，一般取0.85
%输出Data_var为处理后的变量数组
%程序会输出主成分个数和各个主成分变量

    x = zscore(data);    %数据标准化;    
    std=corrcoef(x);   %计算相关系数矩阵
    [vec,val]=eig(std);   %求特征值(val)及特征向量(vec)
    newval=diag(val) ;  %将特征值作成一个新向量
    [y,ii] = sort(newval);   %对特征根进行排序，y为排序结果，ii为索引
    rate = y/sum(y);   %计算贡献率
    sumrate=0;    newi=[];
    for k=length(y):-1:1
        sumrate=sumrate+rate(k);
        newi(length(y)+1-k)=ii(k);
        if sumrate>n
            break;
        end
    end
    fprintf('主成分个数：%g\n\n',length(newi));    %输出主成分个数
    for i=1:1:length(newi)      %计算载荷aa
        for j=1:1:length(y)
            aa(i,j)=sqrt(newval(newi(i)))*vec(j,newi(i));
        end
    end
    aaa=aa.*aa;          %主成分载荷归一化zcfhz
    for i=1:1:length(newi)
        for j=1:1:length(y)
            zcfhz(i,j)=aa(i,j)/sqrt(sum(aaa(i,:)));
        end
    end
    disp('主成分为：')
    for i = 1:length(newi)
        Data_var(:,i) = data(:,1).*zcfhz(i,1) + data(:,2).*zcfhz(i,2) + data(:,3).*zcfhz(i,3) + data(:,4).*zcfhz(i,4) + data(:,5).*zcfhz(i,5) + data(:,6).*zcfhz(i,6);
        fprintf('F%d = %4.2f * x1 + %4.2f * x2 + %4.2f * x3 + %4.2f * x4 + %4.2f * x5 + %4.2f * x6\n',i,zcfhz(i,1),zcfhz(i,2),zcfhz(i,3),zcfhz(i,4),zcfhz(i,5),zcfhz(i,6));
    end

end

输入数组data;行为数据条数，列为指标；输入n,为可接受累计贡献率范围，一般取0.85
输出Data_var为处理后的变量数组；程序会输出主成分个数和各个主成分变量

运行结果：

主成分个数：5

5个主成分变量为：
F1 = -0.08 * x1 + -0.31 * x2 + 0.59 * x3 + 0.29 * x4 + -0.38 * x5 + 0.57 * x6
F2 = -0.55 * x1 + -0.51 * x2 + -0.15 * x3 + -0.64 * x4 + -0.06 * x5 + 0.07 * x6
F3 = 0.06 * x1 + -0.45 * x2 + 0.42 * x3 + 0.11 * x4 + 0.73 * x5 + -0.25 * x6
F4 = 0.81 * x1 + -0.40 * x2 + -0.05 * x3 + -0.35 * x4 + -0.25 * x5 + -0.04 * x6
F5 = 0.18 * x1 + 0.27 * x2 + -0.08 * x3 + -0.31 * x4 + 0.49 * x5 + 0.74 * x6

处理后的变量data_var：（节选）

时间	F1	F2	F3	F4	F5
00:00	3753.24	-1738.50	3269.51	229.46	-124.01
00:40	3737.15	-1701.82	3084.63	281.72	-198.63
01:20	3724.05	-1697.11	3066.68	284.61	-200.35
02:00	3734.99	-1702.77	3086.21	280.34	-196.26
02:40	3753.36	-1710.75	3092.93	277.15	-199.29