PCR与PLS实验报告(第二稿)

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]上的，

，

，，，令随机噪声e服

。输出矩阵

从[0,1]上均匀随机分布，噪声强度系数m=10，Y=X*b+m*e。将随机矩阵

样本数为50；另一部分为预测集

进行建模，用预测集

当K分别取1，100，10000，1000000时散点分布图如下：

岭回归K 1 100 10000 1000000 XxxxX?X?b?，Y分成两部分，一部分为训练集

，

，样本数为50。对训练集对模型进行预测。

（校正标准 R2（复相关 系数) （留一法交叉R2（复相关 系数) 误差）SEC 检验误差 ）SECV 0.5611 -37.8523 0.7634 0.9980 0.9987 0.0547 0.0208 0.0063 0.0042 0.7913 0.8246 0.9856 0.9945 0.0532 0.0051 0.0049 讨论：上述岭回归图形和表格可以看出，当K取越大时，岭回归的回归精度越高。实验过程中发现K取越大越好。

3.1 对岭回归模型进行预测

用

点分布图如下：

岭回归因子K 1 100 10000 1000000

讨论:由上述岭回归预测图形和表格可以看出，当K越大时，岭回归预测模型精度越好。

发现的问题：岭回归得出的回归系数矩阵B与模型给定的系数b之间的误差较大。

3.2 随机噪声对岭回归回归模型的影响

对岭回归进行回归噪声影响分析时，先设定岭回归因子K=10000，调节噪声强度系数m的值，观察m分别取5，10，20，50时回归分布图：

X? （预测标准误差）SEP 0.7388 （复相关系数）R2 -53.5189 0.0965 0.9970 0.9982 0.1081 0.0058 0.0051 数据对模型进行预测，当K分别取1，100，10000，1000000时散

讨论：由上图可以看出，噪声强度系数m取较小时，模型回归精确度比较好。当噪声强度系数m越大时，岭回归的回归精度越差。

3.3 随机噪声对岭回归模型预测的影响

对岭回归模型进行预测分析，令岭回归因子K=10000，调节噪声强度系数m的值，观察m分别取5，10，20，50时回归分布图：

噪声强度系数（m） （预测标准误差）SEP 5 10 20 50 0.0067 0.0115 0.0283 0.0706 （复相关系数）R2 0.9951 0.9907 0.9479 0.5745

讨论：由分布散点图和误差分析表格，可以看出，当噪声强度系数m取值较小时，岭回归预测精度较高，但当m取值越来越大时，岭回归预测精度变差。

四．模型估计以及误差分析

4.1 线性模型进行PCA分析

同样运用岭回归的数据，得到样本集训练集对训练集

,并将样本集分为两部分，一部分为

，样本数为50。

对模型进行预测。根据Xtrain，样本数为50；另一部分为预测集进行建模，用预测集

和最初设定的主元数，以及NIPALS算法，可求得得分矩阵T和载荷矩阵P。

4.1.1 在不同主元下回归误差以及交叉检验误差

计算该模型的回归误差和留一法交叉检验误差，并绘制对应散点图。 主元数(m) （校正标准 误差）SEC R2（复相关 系数) （留一法交叉检验误差 ）SECV 1 3 5 8 0.0985 0.0767 0.0501 0.009 0.8467 0.9235 0.9952 0.9966 0.0876 0.0755 0.0449 0.0083 0.8049 0.9150 0.9908 0.9942 R2（复相关 系数) 当主元数分别为1、3、5和8时，模型回归对应的散点图如下：

4.1.2 对PCR模型进行预测

将输入数据矩阵X分为两堆行列维数相同数据，分别为，，用

来对模型进行预测，由NIPALS算法算出的得分向量和载荷向量，再根据式

（2-6），即可求出预测集的值分布图如下：

XXY。主元数m分别为1，3，5，8时，对应的散点

主元数（m） 1 3 5 8 15 （预测标准误差）SEP 0.0237 0.0162 0.0090 0.0076 0.0064 （复相关系数）R2 0.9525 0.9768 0.9897 0.9952 0.9970 讨论：由模型预测分布散点图和预测误差表格可知，利用建模的主成分越多，模型预测精度越高。但当主成分数和建模样本的组分数相近时，模型预测能力变差，基至不能利用其进行预测。且实验发现对样本矩阵进行均值化处理的模型预测效果没有未进行均值化处理的效果好。

4.1.3 分析不同噪声水平对预测的影响

根据Xtrain和最初设定的主元数8，以及NIPALS算法，可求得得分矩阵T和载荷矩阵P。观察不同噪声强度m下的散点分布图：

噪声强度系数m 5 10 20 50 预测标准误差（SEP） 0.0082 0.0112 0.0231 0.0577

讨论：当噪声水平较小时，PCR预测模型较精确。当噪声水平逐渐增大时，模型预测能力减弱，基至不能进行预测。

复相关系数（R2） 0.9947 0.9874 0.9605 0.6570

4.2 用线性模型进行PLS估计

依旧运用上面的数学模型，获得X，Y。将样本集分，一部分为训练集样本数为50。对训练集

预测。由2.2的PLS1算法可知，对Xtrain训练集矩阵和训练集Y值矩阵Ytrain进行计

算即可得到Xtrain对应载荷矩阵P和X与Y的相关矩阵W，再根据式（2-5）求得B1，由此可求得预测集的预测Y值，具体公式如下：

Ypredict?XpredictRB1?XpredictW(PTW)?1B1

4.2.1 在不同的主元数下分析回归误差和交叉检验误差

主元数(m) （校正标准 R2（复相关 系数) （留一法交误差）SEC 1 3 5 10 0.1038 0.0856 0.0243 0.0098 主元数分别为1，3，5，10时对应的分布散点图如下：

X???,分为维数相同的两部

，样本数为50，另一部分为预测集进行建模，用预测集

，

对模型进行

R2（复相关 系数) 叉检验误差 ）SECV 0.8716 0.9428 0.9673 0.9857 0.0956 0.0714 0.0208 0.0057 0.8426 0.9025 0.9472 0.9896

4.2.2在不同主元数下对模型进行预测

将输入数据矩阵X分为两堆行列维数相同数据，分别为

，

，用

来对模型进行预测，由2.2的PLS1算法可知，对Xtrain训练集矩阵和训练集

Y值矩阵Ytrain进行计算即可得到Xtrain对应载荷矩阵P和X与Y的相关矩阵W，再根据式（2-5）求得B1，由此可求得预测集的预测Y值，具体公式如下：

Ypredict?XpredictRB1?XpredictW(PTW)?1B1

主元数分别为1，3，5，10，18，20时对应的预测散点分布图如下：

X

主元数（m） 1 3 5 10 18 （预测标准误差）SEP 0.0223 0.0073 0.0068 0.0064 0.1101 （复相关系数）R2 0.9485 0.9921 0.9953 0.9961 -0.0545

讨论：由模型回归和预测分布散点图以及误差分析表格可知，当主元数较少时，PLS算法可以获得较好的回归精度和预测精度。但当主元数接近样本组分的个数时，则会使模型回归和预测能力都减弱，甚至不能对模型进行预测。实验仿真发现，对样本进行均值化处理模型预测精度没有未进行均值化处理的高。

4.2.3 分析不同噪声水平在PLS下对预测的影响

根据Xtrain和最初设定的主元数8，以及PLS1算法，可求得预测值。观察不同噪声强度m下的散点分布图：

噪声强度系数m 5 10 20 50 预测标准误差（SEP） 0.0066 0.0133 0.0255 0.0547 复相关系数（R2） 0.9954 0.9867 0.9393 0.7504

讨论：由以上散点分布图和误差表格可知，当噪声水平较小时，PLS下模型的预测能力较好。但当噪声强度逐渐增大时，模型预测能力亦越差。当噪声强度系数过大时到和X中的数据大小相当时，PLS算法不能很好地实现模型预测。

本文来源：https://www.bwwdw.com/article/ymxg.html