PCA与PLS仿真报告
1.实验目的
1) 了解PCA和PLS的基本原理。
2) 利用Matlab程序语言建立真实线性和非线性数学模型，并利用PCR和PLS预估
该模型，并计算对应的校正标准误差（SEC）、交叉检验标准误差（SECV）和预测标准误差（SEP），判断预估模型的准确性。
3) 分析噪声大小和样本矩阵线性关联对模型稳健性的影响。
2.算法介绍
2.1 主成分分析（Principal Component Regression）
主成分回归法（PCR）是采用多元统计中的主成分分析方法（Principal Component Analysis），先对输入矩阵X进行分解，然后选取其中的主成分来进行多元线性回归分析，故称之为主成分回归。
2.1.1主成分分析原理
PCA将矩阵X（n×k）分解为k个向量的外积之和，即：
TTT （2-1） X?t1p1?t2p2???tkpk式中t为得分向量（Score Vector）；p为载荷向量（Loading Vector）或称为主成分。上式也可写成下列矩阵形式：
X?TPT （2-2）
式中T??t1,t2,?,tn?称为得分矩阵；P??p1,p2,?,pn?称为载荷矩阵。各个得分向量之间是正交的，各个载荷向量之间也是正交的，且每个载荷的向量长度都为1。不难得出，ti?Xpi。此式说明了主成分分析的数学意义，即每个得分向量实际上是矩阵X在其
对应载荷向量方向上的投影。向量ti反映了矩阵X在pi方向的覆盖程度。
如果将得分向量按其长度做以下排列：
t1?t2???tk
则载荷向量p1代表矩阵X变化（方差）最大的方向，p2与p1垂直并代表X变化的第二大方向，对X进行主成分分析实际上等效于对Xpk代表X变化的最小的方向。的协方差矩阵XTX进行特征向量分析。矩阵X的载荷向量实际上是矩阵XTX的特征向量。
主成分分析与矩阵的奇异值分解（SVD）也是密切相关的。奇异值分解法可将任意阶实数矩阵分解为三个矩阵的乘积，即：
X?USVT （2-3）
式中S为对角矩阵，是协方差矩阵XTX特征值的平方根；U和VT分别为标准列正交和标准行正交矩阵。实际上，矩阵U与矩阵S的乘积等于主成分分析中的得分矩阵T，矩阵V则等于载荷矩阵P。矩阵X的协方差矩阵的前f个特征值的和除以它的所有特征值之和，称为前f个主成分的累积贡献率，表示前f个主成分解释的数据变化占全部数据变化的比例。选取主成分的个数取决于主成分的累计方差贡献率，通常使累计方差贡献率大于85%所需的主成分数能够代表原始变量所能提高的绝大部分信息。
2.1.2 主成分回归原理
用矩阵X主成分分析得到的前f个得分向量组成矩阵T?[t1,t2,?,tf],代替原始输入矩阵X进行多元线性回归（MLR），便得到了主成分回归模型：
?E Y?TB （2-4）
B的最小二乘法解为：
B??TTT?TTY （2-5）
?1在实际仿真中我们采用NIPALS迭代算法求取主成分，具体算法如下： （1）取X中某列矢量x作为t的起始值：t?x； （2）计算PT?tTX/tTt；
（3）将PT归一化，PT?PT/P；
（4）计算t，t?Xp
（5）比较新的t与旧的t，看是否满足收敛条件。若满足收敛条件（tnew-t?10-3），继续步骤（6），否则跳回步骤（2）；
（6）若已完成计算所需的主成分，则停止计算；否则计算残差阵E：
E?X?tpT
（7）用E代替X，返回步骤（1），求下一个主成分。直到求取到所要求的主成分数时，算法停止。
利用NIPALS迭代算法对训练集矩阵Xtrain训练集矩阵进行主元分析，得到载荷矩阵P和得分矩阵T，再根据式（2-5）求得B，由此可求得预测集的预测Y值，具体公式如下：
Ypred?ictX
至此，主成分回归完成。
predi （2-6） PB2.2 部分最小二乘法回归（Partial least squares regression）
在PCR中，只对输入矩阵X进行分解，消除无用的噪音信息。同样，输出矩阵Y也包含无用信息，应对其作同样的处理，且在分解矩阵X时应考虑矩阵Y的影响。偏最小二乘法PLS就是基于以上思想提出的多元回归方法。
PLS的首先对输入矩阵X和输出矩阵Y进行分解，其模型为：
?E X?TP （2-7）
Y?UQ?F (2-8)
式中T和U分别为X和Y矩阵的得分矩阵；P和Q分别为X和Y矩阵的载荷矩阵；E和F分别为X和Y矩阵的PLS拟合残差矩阵。PLS的第二步是将T和U作线性回归：
U?TB （2-9） B??TTT??1T U （2-10）
在实际仿真中我们采用迭代算法求取主成分，具体算法如下：
PLS2基本算式： （1）u?Y中的第一列； （2）w?XTu/?uTu?； （3）w?w/w （4）t?Xw； （5）c?YTt/?tTt?； （6）c?c/c； （7）u?Yc/?cTc?；
（8）如果收敛（收敛条件：t?told/t?10?8）则转到（9）否则转到（2）； （9）X负载矢量：p?XTt/?tTt?； （10）Y负载矢量：q?YTu/?uTu?； （11）回归（t对u）：v??uTt?/?tTt?； （12）残差：X?X?tpT，Y?Y?vtcT；
（13）计算下一维：转到（1），直到求得所需的维数。
在PLS2的计算公式中若矩阵Y中只有一个目标变量，则PLS方法的整个计算将
会简化，则算式为PLS1算法。 PLS1基本算式： （1）w?XTy/?yTy?； （2）w?w/w； （3）t?Xw；
（4）回归（t对y）：v??tTy?/?tTt?； （5）X负载矢量：p?XTt/?tTt?； （6）残差：X?X?tpT，y?y?vt；
（7）返回（1），开始下一维计算，直到求得所需维数。
不论是PLS1算法还是PLS2算法，对Xtrain训练集矩阵和训练集Y值矩阵Ytrain进
行计算即可得到Xtrain对应载荷矩阵P和X与Y的相关矩阵W，再根据式（2-5）求得
B1，由此可求得预测集的预测Y值，具体公式如下：
Ypredict?XpredictRB1?XpredictW(PTW)?1B1 （2-11）
至此，PLS回归完成。
三．输入共线性分析(岭回归)
设定输入数据矩阵
中
T?1含有共线性量，则由多元线
性回归求回归系数矩阵B时，算求取?XtrainXtrain?的对应行列式XtrainTXtrain的值接近与0，?XtrainTXtrain?对角线上的元素就会很大，因为有
?1这样求得的回归系数?train的值的偏差较大，甚至无意义，因为可以看出当特征值有一个至少趋于0时，B与b的偏差很大。解决上述问题，需用岭回归思想：当
TXtrainXtrain?0时，设想给XtrainTXtrain加上一个正常数矩阵kI(k>0), I为单位矩阵，
?1T?1那么?XtrainXtrain?kI?接近奇异的可能性就会比?XtrainXtrain?接近奇异的可能性小
T得多，因为矩阵XtrainTXtrain?kI按大小排列的特征根为?1?k,?2?k,?,?n?k， 其中
?1,?2,?,?n为矩阵XtrainTXtrain按大小排列的特征根，当?n?0时，?n?k?0，因此用
?train?k???XtrainTXtrain?kI?XtrainTYtrain （3-1）
?1为回归系数的估计应比最小二乘估计?train稳定。常数k?0，且?train?0时，
?trai?nk???对
train，?train?k?可看作对?train的某种压缩。具体算法实现如下:
中的数据矩阵，
均匀随机分布，用rand函数产生。
Xx
为100*1维的服从[0,1]上的，
，