【数学建模】day10-主成分分析

0.

关于主成分分析的详细理解以及理论推导，这篇blog中讲的很清楚。

主成分分析是一种经常使用手段。这应该与因子分析等区别开来，重点在于理解主成分分析的做用以及什么状况下使用主成分分析，本文重点讲解如何使用PCA。

 

1.

主成分分析是一种降维方法。

实际上这个降维是这样作的：原始变量有m维，PCA主成分变量有t维（t<m），那么就至关于把这m维分别往t维上投影。

 

例如咱们要作回归分析，若是自变量众多，彼此之间又具备复杂的相关性，那么咱们考虑对自变量个数进行“减小”。而这个减小不可以丢失有效信息，由此可使用主成分分析。

 

主成分分析的主要思想是，对原始众多自变量进行【线性组合】（这实际上是向PCA方向投影的过程），从而获得新的PCA变量(个数一般比原始变量少)。对于每个新变量，其线性组合的系数向量叫该主成分的方向（重点在于如何获得这个PCA方向）；不一样主成分的方向是正交的，从而保证新的变量彼此不相关，消除了糅合性。

具体来讲，是：

原始变量Xi(i = 1,2,…p)，主成分变量Z(z = 1,2…p)，则：

其中，要求知足：

 

注意：

　　1）主成分分析的结果受量纲的影响，因为各变量的单位可能不同，若是各自改 变量纲，结果会不同，这是主成分分析的大问题，回归分析是不存在这种状况的， 因此实际中能够先把各变量的数据标准化，而后使用协方差矩阵或相关系数矩阵进行分析。（相关系数就是标准化变量的协方差） 

　　2）使方差达到大的主成分分析不用转轴（因为统计软件常把主成分分析和因子 分析放在一块儿，后者每每须要转轴，使用时应注意）。 

　　3）主成分的保留。用相关系数矩阵求主成分时，Kaiser主张将特征值小于1的主成 分予以放弃（这也是SPSS软件的默认值）。 

　　4）在实际研究中，因为主成分的目的是为了降维，减小变量的个数，故通常选取 少许的主成分（不超过5或6个），只要它们能解释变异的70％～80％（称累积贡献率） 就好了。 

 

2.

如何选取主成分？

或者说：重点在于找到主成分方向，由于有了PCA方向，把原始变量向PCA方向上投影就得新的PCA变量。

明确目的在于线性组合原始变量获得：

（z是向量）

在系数平方和为1下，使得z的方差最大。（这样就使得新的变量z1，z2…差别最大，就表明咱们抓住了原始变量的大多数信息）；同时，保证各个主成分变量的系数矩阵c两两正交（这表明咱们使得新的变量彼此不糅合，不相关）。

 

简单推导：

设原始变量数据：

（X称做设计阵）每列表明一个变量指标，每行是一组数据。

 

首先将X标准化。

 

设新变量z = c1x1 + c2x2 +… + cpxp；

系数向量C = (c1,c2,…cp).T

 

系数向量C就是主成分的方向。

 

目的在于使得原始变量X的每一行，在C的方向投影后，在C的方向上达到方差最大。也就是：

　　Z = CX （Z是X向C投影）

　　max E（(Z-E(Z))‘* (Z-E(Z))）

这个求max能够通过变化，转换^{【具体看推导blog，仅仅使用PCA能够略过】}为求max E((X - E(x))' * (X - E(X)))

因为X标准化了（归一化过），E(X) = 0;

从而变为 max E（X.T*X）

而X.T * X 是一个二次型的。

那么，最大值就在：这个矩阵X.T * X 的最大特征值λ对应的特征向量c = [c1,c2...cp]处取得。

c还要通过单位化处理。

也就是说，这就获得了第一主成分的方向c；把X向c上投影获得了第一主成分Z1.

而这个半正定矩阵X.T *X的不一样特征向量是正交的，从而把第二特征向量做为第二主成分的方向。

以此类推。

 

这里，X标准化后，X.T*X / (n-1)其实是原始变量x1,x2…xp的相关系数矩阵。（协方差矩阵）

求解时，其实是计算这个相关系数矩阵。

 

3.

PCA步骤：

　　1. 列出设计阵X，将X标准化处理。

　　2. 求相关系数矩阵R（标准化后，这也就是协方差矩阵）。这里要知道，这个R能够用：

X.T*X / (n-1)求得。也可使用MATLAB求相关系数矩阵命令。

　　3. 求R的特征值（从大到小排列），以及对应的单位标准正交特征向量。

 

　　3. 从特征值大到小选择PCA变量：先选取最大的特征值以及其特征向量，从而构成第一个主成分变量，这个特征向量就是PCA方向。计算累计贡献率：已选取的特征值之和占全部特征值之和的比重。

　　4. 重复步骤3，直到累计贡献率达标或者选取了足够的PCA变量。

　　5. 单纯考虑累积贡献率有时是不够的，还须要考虑选择的主成分对原始变量的贡献 值，咱们用相关系数的平方和来表示.若是选取的主成分为 z1,z2,…zr ，则它们对原 变量 xi 的贡献值为:

　　　　　　　　　　pi = ∑ (r(zj,xi))^2；（注1）

　　　　　（即每个主成分变量与xi的相关系数平方的和。）

　　6. 进而咱们能够用主成分变量对问题作出其余分析（如回归分析等）。

 

注1：两个向量的相关系数就是两向量夹角的余弦，展开来讲就是：

 

4.

一个例子：

 

5.

MATLAB实现

函数使用：

1. 求矩阵X列向量间的相关系数矩阵： r = corrcoef(X) param: 　　X：这个矩阵X的每一列看作一个向量 return： 　　r： 相关系数矩阵。返回每两列之间的相关系数组成的矩阵，对称阵，对角为方差，一个向量的方差就是对向量的每一个数求var。 　　注意使用时没必要对设计阵作标准化处理。由于咱们实际要求的是标准化变量的协方差矩阵，而就是原设计阵的相关系数矩阵。

2. matlab作PCA分析 [vec1,lamda,rate] = pcacov(r) param： 　　r：原始数据的相关系数矩阵 return： 　　vec1: r的特征向量 　　lamda: 对应的特征值 　　rate: 各个主成分的贡献率

3. 累计求和函数cumsum y = cumsum(x,axis) 对x矩阵累计求和 param: 　　axis：轴，axis = 1则对列向量累计求和 　　axis = 2则对行向量累计求和 return: 　　y: 累和矩阵

4. 标准化化处理函数 y = zccore(x) param: 　　x： 矩阵或者向量，标准化处理，如果矩阵，是对列向量标准化处理 return: 　　y： 处理后的矩阵或者向量

5. 主成分分析作线性回归最小二乘估计函数 [c,s,t] = princomp(x) param: 　　x是设计阵 return: 　　c:  对主成分变量作多元线性回归分析，回归方程的系数 　　s:  这个是作主成分分析获得的特征向量矩阵，每一列是一个特征向量（单位化） 　　t：相应的特征值

 

 

 

应用主成分分析+回归分析案例:

分别作主成分回归以及原始变量直接回归。

数据保存：

data.txt

7 26 6 60 78.5 
 1 29 15 52 74.3 
 11 56 8 20 104.3 
 11 31 8 47 87.6  
 7 52 6 33 95.9 
 11 55 9 22 109.2 
 3 71 17 6 102.7  
 1 31 22 44 72.5  
 2 54 18 22 93.1 
 21 47 4 26 115.9 
 1 40 23 34 83.8 
 11 66 9 12 113.3  
 10 68 8 12 109.4

matlab求解以下：

（注：也能够直接用princomp函数，更简单，下面使用标准的pcacov函数）

 1 clc,clear  2  load data.txt %导入数据  3 [m,n] = size(data);  4  x0 = data(:,[1:n-1]);  5  y0 = data(:,n);  6  hg1 = [ones(m,1),x0] \ y0; %普通多元线性回归系数.列向量  7 hg1 = hg1' %显示
 8 fprintf('y = % f',hg1(1));  9  for i = 2:n 10      if hg1(i) > 0
11          fprintf('+ % f*x% d',hg1(i),i-1); 12      else
13          fprintf('% f*x% d',hg1(i),i-1); 14  end 15  end 16 
17 fprintf('\n'); 18  r = corrcoef(x0) %相关系数矩阵 19 xd = zscore(x0); 20  yd = zscore(y0); %标准化处理 21 [vec1,lamda,rate] = pcacov(r) %PCA 22  f = repmat(sign(sum(vec1)),size(vec1,1),1); 23 %产生与vec1同维数的元素为+1/-1的矩阵 24 vec2 = vec1.*f %修改特征值正负号，使得特征值的全部份量和为+
25 contr = cumsum(rate) %计算累计贡献率 26 df = xd * vec2; %计算全部主成分的得分 27 num = input('请输入主成分个数：') 28 hg21 = df(:,[1:num]) \ yd %主成分变量回归系数 29 hg22 = vec2(:,1:num) * hg21 %标准化变量的回归系数 30 hg23 = [mean(y0) - std(y0)*mean(x0)./std(x0)*hg22,std(y0)*hg22'./std(x0)]
31  % 转换，求原始变量的回归方程系数 32 
33 fprintf('y = % f',hg23(1)); 34  for i = 2:n 35      if hg23(i) > 0
36          fprintf('+ % f*x% d',hg23(i),i-1); 37      else
38          fprintf('% f*x% d',hg23(i),i-1); 39  end 40  end 41  fprintf('\n') 42  %下面计算两种回归分析的剩余标准差 43 rmse1=sqrt(sum((x0*hg1(2:end)'+hg1(1)-y0).^2)/(m-n)) 
44  rmse2=sqrt(sum((x0*hg23(2:end)'+hg23(1)-y0).^2)/(m-num)) 

结果以下：

相关系数矩阵：

PCA分析获得的：特征值、特征向量、贡献率等结果：

作代码所示处理后的特征向量：

累计贡献率：

能够看到，最后一个变量几乎没有贡献，使用前三个PCA成分获得:（其实两个也能够）

其中，hg21是主成分变量回归系数、hg22是标准化X的回归系数(没有常数项)，hg23是原始X的回归系数（第一个是常数项）

从而，PCA后的结果为：

而在原始变量上直接作回归分析的结果为：

 

代码中还计算了两个剩余标准差：

可见，PCA后再回归的剩余标准差要小，效果更好，这也证实了使用PCA是有效的。