重庆工商大学数学建模算法讲义第30章 偏最小二乘回归

-531-第三十章偏昀小二乘回归在实际问题中，经常遇到需要研究两组多重相关变量间的相互依赖关系，并研究用一组变量（常称为自变量或预测变量）去预测另一组变量（常称为因变量或响应变量），除了昀小二乘准则下的经典多元线性回归分析（MLR），提取自变量组主成分的主成分回归分析（PCR）等方法外，还有近年发展起来的偏昀小二乘（PLS）回归方法。偏昀小二乘回归提供一种多对多线性回归建模的方法，特别当两组变量的个数很多，且都存在多重相关性，而观测数据的数量（样本量）又较少时，用偏昀小二乘回归建立的模型具有传统的经典回归分析等方法所没有的优点。偏昀小二乘回归分析在建模过程中集中了主成分分析，典型相关分析和线性回归分析方法的特点，因此在分析结果中，除了可以提供一个更为合理的回归模型外，还可以同时完成一些类似于主成分分析和典型相关分析的研究内容，提供更丰富、深入的一些信息。本章介绍偏昀小二乘回归分析的建模方法；通过例子从预测角度对所建立的回归模型进行比较。§1偏昀小二乘回归考虑p个变量pyyy,,,21与m个自变量mxxx,,,21的建模问题。偏昀小二乘回归的基本作法是首先在自变量集中提出第一成分1t（1t是mxx,,1的线性组合，且尽可能多地提取原自变量集中的变异信息）；同时在因变量集中也提取第一成分1u，并要求1t与1u相关程度达到昀大。然后建立因变量pyy,,1与1t的回归，如果回归方程已达到满意的精度，则算法中止。否则继续第二对成分的提取，直到能达到满意的精度为止。若昀终对自变量集提取r个成分rttt,,,21，偏昀小二乘回归将通过建立pyy,,1与rttt,,,21的回归式，然后再表示为pyy,,1与原自变量的回归方程式，即偏昀小二乘回归方程式。为了方便起见，不妨假定p个因变量pyy,,1与m个自变量mxx,,1均为标准化变量。因变量组和自变量组的n次标准化观测数据阵分别记为⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=npnpyyyyF##11110，⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=nmnmxxxxE##11110偏昀小二乘回归分析建模的具体步骤如下：-532-（1）分别提取两变量组的第一对成分，并使之相关性达昀大。假设从两组变量分别提出第一对成分为1t和1u，1t是自变量集TmxxX),,(1=的线性组合：XwxwxwtTmm111111=++=，1u是因变量集TpyyY),,(1=的线性组合：YvyvyvuTpp111111=++=。为了回归分析的需要，要求：①1t和1u各自尽可能多地提取所在变量组的变异信息；②1t和1u的相关程度达到昀大。由两组变量集的标准化观测数据阵0E和0F，可以计算第一对成分的得分向量，记为1ˆt和1ˆu：⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡==1111111111101ˆnmnmnmttwwxxxxwEt####⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡==1111111111101ˆnpnpnpuuvvyyyyvFu####第一对成分1t和1u的协方差),(Cov11ut可用第一对成分的得分向量1ˆt和1ˆu的内积来计算。故而以上两个要求可化为数学上的条件极值问题：⎪⎩⎪⎨⎧====⇒==1,1max,ˆ,ˆ21112111001101011vvv乘数法，问题化为求单位向量1w和1v，使⇒=10011vFEwTTθ昀大。问题的求解只须通过计算mm×矩阵0000EFFEMTT=的特征值和特征向量，且M的昀大特征值为21θ，相应的单位特征向量就是所求的解1w，而1v可由1w计算得到100111wEFvTθ=。（2）建立pyy,,1对1t的回归及mxx,,1对1t的回归。假定回归模型为-533-⎪⎩⎪⎨⎧+=+=11101110ˆˆFuFEtETTβα其中Tm),,(1111ααα=，Tp),,(1111βββ=分别是多对一的回归模型中的参数向量，1E和1F是残差阵。回归系数向量11,βα的昀小二乘估计为⎪⎩⎪⎨⎧==2110121101ˆˆˆˆttFttETTβα，称11,βα为模型效应负荷量。（3）用残差阵1E和1F代替0E和0F重复以上步骤。记TtE110ˆˆα=，TtF110ˆˆβ=，则残差阵001ˆEEE−=，001ˆFFF−=。如果残差阵1F中元素的绝对值近似为0，则认为用第一个成分建立的回归式精度已满足需要了，可以停止抽取成分。否则用残差阵1E和1F代替0E和0F重复以上步骤即得：Tm),,(2212=；Tpvvv),,(2212=分别为第二对成分的权数。而212ˆwEt=，212ˆvFu=为第二对成分的得分向量。22212ˆˆttET=α，22212ˆˆttFT=β分别为YX,的第二对成分的负荷量。这时有⎪⎩⎪⎨⎧++=++=222110222110ˆˆˆˆFttFEttETTTTββαα（4）设mn×数据阵0E的秩为),1min(mnr−≤，则存在r个成分rttt,,,21，使得⎪⎩⎪⎨⎧+++=+++=rTrrTrTrrTFttFEttEββααˆˆˆˆ110110把mkmkkxwxwt++=11（rk,,2,1=），代入rrttYββ++=11，即得p个因变量的偏昀小二乘回归方程式-534-mjmjjxaxay++=11，（mj,,2,1=）（5）交叉有效性检验。一般情况下，偏昀小二乘法并不需要选用存在的r个成分rttt,,,21来建立回归式，而像主成分分析一样，只选用前l个成分（rl≤），即可得到预测能力较好的回归模型。对于建模所需提取的主成分个数l，可以通过交叉有效性检验来确定。每次舍去第i个观测（ni,,2,1=），用余下的1−n个观测值按偏昀小二乘回归方法建模，并考虑抽取h个成分后拟合的回归式，然后把舍去的第i个观测点代入所拟合的回归方程式，得到),,2,1(pjyj=在第i个观测点上的预测值)(ˆ)(hyji。对ni,,2,1=重复以上的验证，即得抽取h个成分时第j个因变量),,2,1(pjyj=的预测误差平方和为∑=−=nijiijjhyyh12)())(ˆ()(PRESS（pj,,2,1=）TpyyY),,(1=的预测误差平方和为∑==pijhh1)(PRESS)(PRESS。另外，再采用所有的样本点，拟合含h个成分的回归方程。这时，记第i个样本点的预测值为)(ˆhyij，则可以定义jy的误差平方和为∑=−=niijijjhyyh12))(ˆ()(SS定义Y的误差平方和为∑==pjjhh1)(SS)(SS当)(PRESSh达到昀小值时，对应的h即为所求的成分个数。通常，总有)(PRESSh大于)(SSh，而)(SSh则小于)1(SS−h。因此，在提取成分时，总希望比值)1(SS)(PRESS−hh越小越好；一般可设定限制值为0.05，即当-535-2295.0)05.01()1(SS)(PRESS=−≤−hh时，增加成分ht有利于模型精度的提高。或者反过来说，当295.0)1(SS)(PRESS−hh时，就认为增加新的成分ht，对减少方程的预测误差无明显的改善作用。为此，定义交叉有效性为)1(SS)(PRESS12−−=hhQh，这样，在建模的每一步计算结束前，均进行交叉有效性检验，如果在第h步有0985.095.0122=−hQ，则模型达到精度要求，可停止提取成分；若0975.02≥hQ，表示第h步提取的ht成分的边际贡献显著，应继续第1+h步计算。§2一种更简洁的计算方法上节介绍的算法原则和推导过程的思路在目前的文献中是昀为常见的。然而，还有一种更为简洁的计算方法，即直接在10,,−rEE矩阵中提取成分rtt,,1(mr≤)。要求ht能尽可能多地携带X中的信息，同时，ht对因变量系统0F有昀大的解释能力。注意，无需在0F中提取成分得分hu，这可以使计算过程大为简化，并且对算法结论的解释也更为方便。偏昀小二乘法的简记算法的步骤如下：（1）求矩阵0000EFFETT昀大特征值所对应的特征向量1w，求得成分XwtT11=，计算成分得分向量101ˆwEt=，和残差矩阵TtEE1101ˆα−=，其中21101ˆˆttET=α。（2）求矩阵1001EFFETT昀大特征值所对应的特征向量2w，求得成分XwtT22=，计算成分得分向量212ˆwEt=，和残差矩阵TtEE2212ˆα−=，其中22212ˆˆttET=α。#（r）至第r步，求矩阵1001−−rTTrEFFE昀大特征值所对应的特征向量rw，求得成分XwtTrr=，计算成分得分向量rrrwEt1ˆ−=。-536-如果根据交叉有效性，确定共抽取r个成分rtt,,1可以得到一个满意的预测模型，则求0F在rttˆ,,ˆ1上的普通昀小二乘回归方程为rTrrTFttF+++=ββˆˆ110把mkmkkxwxwt*1*1++=（rk,,2,1=），代入rrttYββ++=11，即得p个因变量的偏昀小二乘回归方程式mjmjjxaxay++=11，（mj,,2,1=）。这里的*hw满足*0ˆhhwEt=，∏−=−=11*)(hjhTjjhwwIwα。§3案例分析本节采用兰纳胡德（Linnerud）给出的关于体能训练的数据进行偏昀小二乘回归建模。在这个数据系统中被测的样本点，是某健身俱乐部的20位中年男子。被测变量分为两组。第一组是身体特征指标X，包括：体重、腰围、脉搏。第二组变量是训练结果指标Y，包括：单杠、弯曲、跳高。原始数据见表1。表1体能训练数据No体重(x1)腰围(x2)脉搏(x3)单杠(y1)弯曲(y2)跳高(y3)119136505162602189375221106031933858121011014162356212105375189354613155586182365641014272113856810138816734606125409176317415200401015433561725125011169345017120381216633521321011513154346414215105142474650150501519336466703116202376212210120-537-1717637544602518157325211230801915633541522573201383368211043均值标准差178.624.690535.43.20256.17.21049.455.2863145.5562.566670.351.2775表2给出了这6个变量的简单相关系数矩阵。从相关系数矩阵可以看出，体重与腰围是正相关的；体重、腰围与脉搏负相关；而在单杠、弯曲与跳高之间是正相关的。从两组变量间的关系看，单杠、弯曲和跳高的训练成绩与体重、腰围负相关，与脉搏正相关。表2相关系数矩阵10.8702-0.3658-0.3897-0.4931-0.22630.87021-0.3529-0.5522-0.6456-0.1915-0.3658-0.352910.15060.2250.0349-0.3897-0.55220.150610.69570.4958-0.4931-0.64560.2250.695710.6692-0.2263-0.19150.03490.49580.66921利用如下的MATLAB程序：clc,clearloadpz.txt%原始数据存放在纯文本文件pz.txt中mu=mean(pz);sig=std(pz);%求均值和标准差rr=corrcoef(pz);%求相关系数矩阵data=zscore(pz);%数据标准化n=3;m=3;%n是自变量的个数,m是因变量的个数x0=pz(:,1:n);y0=pz(:,n+1:end);e0=data(:,1:n);f0=data(:,n+1:end);num=size(e0,1);%求样本点的个数chg=eye(n);%w到w*变换矩阵的初始化fori=1:n%以下计算w，w*和t的得分向量，matrix=e0'*f0*f0'*e0;[vec,val]=ei