算法模型 | 主成成分分析PCA+实例详解

主成成分分析，顾名思义就是要对描述一件事物的多个主成分（指标）进行分析，然后对这多个指标进行压缩简化（降维），中间会造成数据的损失。

一、介绍

假设想描述如下各地区的居民幸福感，不同地方就是样本个数，不同的支出方面就是指标

	食品	教育	交通	住房	其他
北京	100	200	100	100	200
深圳	200	200	100	100	200
广州	100	200	400	100	200
上海	300	200	100	200	200

所以有如下定义

这里的Zm就是经过PCA之后的新生成的综合性指标，是对前面p个指标的综合性描述，m的大小一定小于p的，因为实现了指标压缩（降维）。所以我们的目标就是来找出若干个Z的表达式。

具体的推导流程可以参考教材：《应用多元统计分析》王学民

二、代码流程

2.1 标准化

这一步是为了去量纲，也就是消除单位的影响，使它们在计算上变得有意义

X1=zscore(x);   % matlab内置的标准化函数（x-mean(x)）/std(x)

2.2 计算协方差矩阵

R = cov(X1);

简单解释，x1和x1的相关性是1，x1和x2的相关性是0.79，所以会按对角线对称。

2.3 计算特征值特征向量

这里的eig返回两个参数

[V,D] = eig(R);  % V 特征向量矩阵  D 特征值构成的对角矩阵

2.4 贡献率累计贡献率

算出来的a1 a2 就是对应的特征向量，按贡献率从大到小一次排列（才能计算出最优解），特征值累加起来一定等于指标数。而特征值除以指标数就等于贡献率，累加贡献率是对前面到现在这个位置的贡献率的累加，超过80%就可以保留，所以如下图保留的就是3个，这三个综合性指标就能很好的概括整个数据的情况。在这一部分可以描述多一点，并撰写到论文中。

lambda = diag(D);  % diag函数用于得到一个矩阵的主对角线元素值(返回的是列向量)
lambda = lambda(end:-1:1);  % 因为lambda向量是从小大到排序的，我们将其调个头
contribution_rate = lambda / sum(lambda);  % 计算贡献率
cum_contribution_rate = cumsum(lambda)/ sum(lambda);   % 计算累计贡献率  cumsum是求累加值的函数
disp('特征值为：')
disp(lambda')  % 转置为行向量，方便展示
disp('贡献率为：')
disp(contribution_rate')
disp('累计贡献率为：')
disp(cum_contribution_rate')
disp('与特征值对应的特征向量矩阵为：')

2.5 保留主成分的值

根据上面的累计贡献率，超过80%的一般需要多少个我们就选择多少个，然后算出的F（也就是第一点提到的Z1…Zm）就可以用spss或stata进行聚类或者回归

m =input('请输入需要保存的主成分的个数:  ');
F = zeros(n,m);  %初始化保存主成分的矩阵（每一列是一个主成分）
for i = 1:m
    ai = V(:,i)';   % 将第i个特征向量取出，并转置为行向量 
    Ai = repmat(ai,n,1);   % 将这个行向量重复n次，构成一个n*p的矩阵
    F(:, i) = sum(Ai .* X1, 2);  % 注意，对标准化的数据求了权重后要计算每一行的和
end

三、实例

接到一个项目，是要用Olivetti压缩图片，其实和matlab代码一样，我们就用matlab的来，委托方的老师发了一个mat文件，我们只用里面的X这个举证，大致意思是X是400x10304个数据，每一个行是一张图片一共有400张，我们用10个或是100个综合性指标来描述这个人的轮廓，此时的均方误差是多少。

这里需要特别注意的是，我们是要压缩图片的样本数，而不是压缩每一张图，也就是要把X进行转置，如果没有理解好可能就会导致直接matlab跑个十分钟才能出结果。

这个问题和上面统计学的可能不太一样，理解起来也需要花一些时间，但没有要求说是计算出至少需要多少个指标并进行描述，使代码也简单很多，不需要上面的2.3 2.4代码部分，在最后加一个计算矩阵的均方误差就好。

A=std2(F);%均方误差
disp('均方误差为：')
disp(A);