(真正的好东西)偏最小二乘回归=多元线性回归分析+典型相关分析+主成分分析教学内容

1、（真正的好东西）偏最 二乘回归=多元线性 归分析+典型相关分 析+主成分分析偏最小二乘回归是一种新型的 多元统计数据分析方法，它与1983年由伍德 和阿巴诺等人首次提出。近十年来，它在理论、方法和应用方面都得到了迅速 的发展。密西根大学的弗耐尔教授称偏最小二乘回归为第二代回归分析方法。偏最小二乘回归方法在统计应用中的重要性主要的有以下几个方面：(1)偏最小二乘回归是一种 多因变量对多自变量 的回归建模方法。(2)偏最小二乘回归可以较好地解决许多以往用普通多元回归无法解决的 问题。在普通多元线形回归的应用中，我们常受到许多限制。最典型的问题就 是自变量之间的多重相关性。如果采用普通的最小二乘方法

2、，这种变量多重相 关性就会严重危害参数估计，扩大模型误差，并破坏模型的稳定性。变量多重 相关问题十分复杂，长期以来在理论和方法上都未给出满意的答案，这一直困 扰着从事实际系统分析的工作人员。在偏最小二乘回归中开辟了一种有效的技 术途径，它利用对系统中的数据信息进行分解和筛选的方式，提取对因变量的 解释性最强的综合变量，辨识系统中的信息与噪声，从而更好地克服变量多重 相关性在系统建模中的不良作用。(3)偏最小二乘回归之所以被称为第二代回归方法，还由于它可以实现多 种数据分析方法的综合应用。偏最小二乘回归二多元线性回归分析+典型相关分析+主成分分析由于偏最小二乘回归在建模的同时实现了数据结构的简化

3、，因此，可以在 二维平面图上对多维数据的特性进行观察，这使得偏最小二乘回归分析的图形 功能十分强大。在一次偏最小二乘回归分析计算后，不但可以得到多因变量对多自变量的回归模型，而且可以在平面图上直接观察两组变量之间的相关关 系，以及观察样本点间的相似性结构。这种高维数据多个层面的可视见性，可 以使数据系统的分析内容更加丰富，同时又可以对所建立的回归模型给予许多 更详细深入的实际解释。一、 偏最小二乘回归的建模策略原理 方法1.1 建模原理设有q个因变量 yi,yq和p自变量 x1,Xp。为了研究因变量和自变量的统计关系,我们观测了n 个样本点,由此构成了自变量与因变量的数据表X= %,Xp和.Y

4、= yi,yq。偏最小二乘回归分别在X与Y中提取出成分ti和Ui (也就是说，ti是Xi，.,Xp的线形组合，Ui是y1,.,yq的线形组合).在提取这两 个成分时,为了回归分析的需要,有下列两个要求:(1) ti和Ui应尽可能大地携带他们各自数据表中的变异信息；(2) ti 与 Ui 的相关程度能够达到最大。这两个要求表明，ti和ui应尽可能好的代表数据表 X和Y,同时自变量的成分 ti 对因变量的成分Ui 又有最强的解释能力。在第一个成分ti 和 Ui 被提取后，偏最小二乘回归分别实施X 对 ti 的回归以及 Y 对 Ui 的回归。如果回归方程已经达到满意的精度，则算法终止；否则,将利用X

5、被ti解释后的残余信息以及Y被ti解释后的残余信息进行第二轮的成分 提取。如此往复，直到能达到一个较满意的精度为止。若最终对X 共提取了 m个成分ti tm偏最小二乘回归将通过实施yk对ti tm 的回归，然，后再表达成y k关于原变量x 1 x的回归方程，k=i,2,q。X i X m，1.2计算方法推导为了数学推导方便起见，首先将数据做标准化处理。X经标准化处理后的数据矩阵记为Eo=(Eo1 ,，E0P)np , Yj经标准化处理后的数据矩阵记为F0 =( F01 ,，F0q ) n p。第一步 记ti是Eo的第一个成分，Wi是Eo的第一个轴，它是一个单位向量，既 l|w1H=1。记Ui是

6、Fo的第一个成分，U1 = F0 clo是Fo的第一个轴，并且|。|二1。如果要t1 U1能分别很好的代表X与Y中的数据变异信息，根据主成分分析原理，应该有Var(U1) maxVar( t1) max另一方面，由于回归建模的需要，又要求 G对5有很大的解释能力，有典型相关分析的思路，t1与5的相关度应达到最大值，既r (t1 U1) max因此，综合起来，在偏最小二乘回归中，我们要求 t1与U1的协方差达到最大，u1) maxCov(t1U1)= , Var(t)Var(5)r(t1正规的数学表述应该是求解下列优化问题，既maxW1，c1E o W1, Fo c1s.tW1 W11'

7、c1 c11值。因此，将在|wi|2二i和11cl |2二1的约束条件下，去求(w1 E0 F0 c1)的最大如果采用拉格朗日算法，记s=W Eo % C1 (Wi1 1)对s分别求关于W1c11和2的偏导并令之为零,2 (C1 C1-1)有wi=EoF0 C1 _ 2 1 W1=0(1-2)由式(1-2)(1-5)，可以推出s '7 = F0E0W1 _2 2c1=0 c1s , (W1 w1 - 1)=01S .'一一二一(C1 C1 -1)=02''2 12 2 w 1E 0F0c1E0W1,F0C1(1-3)(1-4)(1-5)记12 12 2 W1E0

8、F0C1,所以，1正是优化问题的目标函数值.把式(1-2)和式(1-3)写成'E 0F0c11W1(1-6)一 F 0E0W11a(1-7)将式(1-7)代入式(1-6)，有2E 0 F0 F 0 E0W11 W1(1-8)同理，可得l' L l' L2F 0 E0E 0 F0G1 G(1-9). t. »» 9 、 可见,Wi是矩阵E 0F0F 0E0的特征向量，对应的特征值为1 . 1是目标函数值，它要求取最大值，所以，wi是对应于E0FoF '0E0矩阵最大特征值的单位特征向量.而另一方面，Ci是对应于矩阵F'oEoE'

9、oF。最大特征值12的单位特征向量.求得轴Wi和Ci后,即可得到成分tiEoWiuiFoci然后，分别求Eo和Fo对ti, Ui的三个回归方程EotiPiEi(i-i0)Fouiqi(i-ii)FotiriFi(i-i2)式中，回归系数向量是E otiPi211ti II2(i-i3)F oui qi 2lluill2(i-i4)F oti2|ti |2(i-i5)而Ei, F i,Fi分别是三个回归方程的残差矩阵.第二步 用残差矩阵Ei和Fi取代Eo和Fo然后，求第二个轴W2和C2以及第个成分t2, U2,有t2= Ei W2U2 = Fi C2t2 ,u2w 2E iF1c2W2是对应于矩

10、阵E'lFiF'lEi最大特征值2 2的特征值，C2是对应于矩阵F'iEiE'iFi最大特征值的特征向量.计算回归系数,P2E it2211t2 if，F it2因此，有回归方程El12 P 2 E2F1 t2r 2 F2如此计算下去，如果X的秩是A,则会有一 ''Eo tiP1 tAP A(1-16). '. '，一 一.Fo tir itAAFa(1-17)由于,ti, ,tA均可以表示成E01,Eop的线性组合，因此，式(1-17)还可以还原成ykFok关于Xj* Eok的回归方程形式，即*yk*k1X1kpX pFAkk=

11、1,2,. ；qFAk是残差距阵Fa的第k列。1.3交叉有效性下面要讨论的问题是在现有的数据表下，如何确定更好的回归方程。在许多情形下，偏最小二乘回归方程并不需要选用全部的成分 t1, , tA进行回3建模，而是可以象在主成分分析一样，采用截尾的方式选择前 m个成分(m A,A 秩(X),仅用这m个后续的成分t1, ,tm就可以得到一个预测性较好的模型。事实上，如 果后续的成分已经不能为解释F。提供更有意义的信息时，采用过多的成分只会破坏对统计趋势的认识,引导错误的预测结论。在多元回归分析一章中,我们曾在调整复测定系数的内容中讨论过这一观点。下面的问题是怎样来确定所应提取的成分个数。在多元回归

12、分析中,曾介绍过用抽样测试法来确定回归模型是否适于预测应用。我们把手中的数据分成两部分:第一部分用于建立回归方程,求出回归系数估计量bB,拟合值兔以及残差均方和?B ；再用第二部分数据作为实验点，代入刚才 所求得的回归方程，由此求出外和?t o 一般地，若有?T ?B,则回归方程会有更 好的预测效果。若?T2?B2 ,则回归方程不宜用于预测。在偏最小二乘回归建模中 ,究竟应该选取多少个成分为宜,这可通过考察增加一个新的成分后,能否对模型的预测功能有明显的改进来考虑。采用类似于抽样测试法的工作方式，把所有n个样本点分成两部分：第一部分除去某个样本点i的所有样本点集合(共含n-1个样本点)，用这部

13、分样本点并使用h个成分拟合一个 回归方程；第二部分是把刚才被排除的样本点i代入前面拟合的回归方程，得到yj 在样本点i上的拟合值 商(”对于每一个i =1,2,n,重复上述测试，则可以定义 yj 的预测误差平方和为PRESShj ,有n2PRESShj(yijy?hj( i)2(1-18)i1定义 Y 的预测误差平方和为PRESSh ,有pPRESShPRESShj(1-19)j1显然 ,如果回归方程的稳健性不好,误差就很大,它对样本点的变动就会十分敏感,这种扰动误差的作用,就会加大PRESSh 的值。另外,再采用所有的样本点，拟合含h个成分的回归方程。这是，记第i个样本点的预测值为 外，则可

14、以记yj的误差平方和为SSj,有 nsshj(yj yhji )2(1-20)i 1定义Y的误差平方和为SSh，有pSShSShj(1-21)j 1一般说来，总是有PRESSh大于SSh，而SSh则总是小于SSh 1。下面比较SSh 1和PRESSh。SSh 1是用全部样本点拟合的具有h-1个成分的方程的拟合误差；PRES0增加了一个成分th,但却含有样本点的扰动误差。如果 h个成分的回归方程的含扰动误差能在一定程度上小于(h-1)个成分回归方程的拟合误差，则认为增加一个成分th,会使预测结果明显提高。因此我们希望(PRESS /SSh 1)的比值能越小越好。在SIMCA-P软件中，指定(PR

15、ESSh / SSh 1)0.952即JPRESS 0.95南二时，增加成分th就是有益的；或者反过来说，当JPRESS 0.95、SSh 1时,就认为增加新的成分th,对减少方程的预测误差无明显 的改善作用.另有一种等价的定义称为交叉有效性。对每一个变量 y定义Qhk 1PRES8SSh 1)k(1-22)PRESSh对于全部因变量Y,成分t h交叉有效性定义为QhqPRESShkk 1SS；h 1)kSS；h 1)(1-23)用交叉有效性测量成分th对预测模型精度的边际贡献有如下两个尺度。(1) 当Q2 (1 0.952) 0.0975时，th成分的边际贡献是显著的。显而易见，Q； 0.0

16、975与(PRESSh/S&i) 0.952是完全等价的决策原则。 对于k=1,2,q,至少有一个k,使得Qh20.0975这时增加成分th ,至少使一个因变量yk 的预测模型得到显著的改善, 因此 ,也可以考虑增加成分th 是明显有益的。明确了偏最小二乘回归方法的基本原理、方法及算法步骤后，我们将做实证分析。附录function w=maxdet(A)%求矩阵的最大特征值v,d=eig(A);n,p=size(d);d1=d*ones(p,1);d2=max(d1);i=find(d1=d2);w=v(:,i);%function c,m,v=norm1(C)%对数据进行标准化处理n

17、,s=size(C);for i=1:nfor j=1:sc(i,j)=(C(i,j)-mean(C(:,j)/sqrt(cov(C(:,j);endendm=mean(C);for j=1:sv(1,j)=sqrt(cov(C(:,j);end%function t,q,w,wh,f0,FF=fun717(px,py,C)%px自变量的输入个数%py输入因变量的个数。%C输入的自变量和因变量组成的矩阵%t提取的主成分%q为回归系数。%w最大特征值所对应的特征向量。%wh处理后的特征向量%f0回归的标准化的方程系数%FF原始变量的回归方程的系数c=norm1(C);%norm1为标准化函数截取

18、标准化的因变量求最大特征向量提取主成分%y=c(:,px+1:px+py);%E0=c(:,1:px);F0=c(:,px+1:px+py);A=E0'*F0*F0'*E0;w(:,1)=maxdet(A);%t(:,1)=E0*w(:,1);%E(:,1:px)=E0-t(:,1)*(E0'*t(:,1)/(t(:,1)'*t(:,1)'获得回归系数p(:,1:px)=(E0'*t(:,1)/(t(:,1)'*t(:,1)'for i=0:px-2B(:,px*i+1:px*i+px)=E(:,px*i+1:px*i+px)&#

19、39;*F0*F0'*E(:,px*i+1:px*i+px );w(:,i+2)=maxdet(B(:,px*i+1:px*i+px);% maxdet 为求最大特征值的函数t(:,i+2)=E(:,px*i+1:px*i+px)*w(:,i+2);p(:,px*i+px+1:px*i+2*px)=(E(:,px*i+1:px*i+px)'*t(:,i+2)/(t(:,i+2)'*t( :,i+2)'E(:,px*i+px+1:px*i+2*px)=E(:,px*i+1:px*i+px)- t(:,i+2)*(E(:,px*i+1:px*i+px)'*t

20、(:,i+2)/(t(:,i+2)'*t(:,i+2)'endfor s=1:pxq(:,s)=p(1,px*(s-1)+1:px*s)'endn,d=size(q);for h=1:pxiw=eye(d);for j=1:h-1iw=iw*(eye(d)-w(:,j)*q(:,j)');endwh(:,h)=iw*w(:,h);endfor j=1:pyzr(j,:)=(regress1(y(:,j),t)'% 求回归系数endfor j=1:pxfori=1:py%生成标准化变量的方程的系数矩阵w1=wh(:,1:j);zr1=(zr(i,1:j)&

21、#39;f0(i,:,j)=(w1*zr1)'endnormxy,meanxy,covxy=norm1(C);%normxy 标准化后的数据矩阵%meanxyj一列的均值%covxy每一列的方差ccxx=ones(py,1)*meanxy(1,1:px);ccy=(covxy(1,px+1:px+py)'*ones(1,px);ccx=ones(py,1)*(covxy(1,1:px);ff=ccy.*f0(:,:,j)./ccx;fff=-(sum(ccy.*ccxx.*f0(:,:,j)./ccx)')-meanxy(1,px+1:px+py)'FF(:,:

22、,j)=fff,ff;%生成原始变量方程的常数项和系数矩阵end%function r,Rdyt,RdYt,RdYtt,Rdytt,VIP=fun8y(px,py,c)X=c(:,1:px);Y=c(:,px+1:px+py);x=norm1(X);y=norm1(Y);t,q,w=fun717(px,py,X,Y);r1=corrcoef(y,t);r=r1(py+1:px+py,1:py)'Rdyt=r.A2;RdYt=mean(Rdyt)for m=1:pxRdYtt(1,m)=sum(RdYt(1,1:m)');endfor j=1:pyfor m=1:pyRdytt(

23、j,m)=sum(Rdyt(j,1:m)');endendfor j=1:pxfor m=1:pxRd(j,m)=RdYt(1,1:m)*(w(j,1:m).A2)');endendfor j=1:pxVIP(j,:)=sqrt(px*ones(1,px)./RdYtt).*Rd(j,:);end%function r,Rdxt,RdXt,RdXtt,Rdxtt=fun8x(px,py,c)X=c(:,1:px);Y=c(:,px+1:px+py);x=norm1(X);y=norm1(Y);t,q,w=fun717(px,py,X,Y);r1=corrcoef(x,t);r=

24、r1(px+1:px+px,1:px)'Rdxt=r.A2;RdXt=mean(Rdxt);for m=1:pxRdXtt(1,m)=sum(RdXt(1,1:m)');endfor j=1:pxfor m=1:pxRdxtt(j,m)=sum(Rdxt(j,1:m)');endend% for j=1:px% for m=1:px%Rd(j,m)=RdXt(1,1:m)*(w(j,1:m).A2)');% end% end% for j=1:px% VIP(j,:)=sqrt(px*ones(1,px)./RdYtt).*Rd(j,:);% end% func

25、tion t,u=TU(px,py,C)%t 提取的自变量的主成分% u 提取的因变量的主成分c=norm1(C);y=c(:,px+1:px+py);E0=c(:,1:px);F0=c(:,px+1:px+py);A=E0'*F0*F0'*E0;w(:,1)=maxdet(A);t(:,1)=E0*w(:,1);B=F0'*E0*E0'*F0;cc(:,1)=maxdet(B);u(:,1)=F0*cc(:,1);% function drew(px,py,c)X=c(:,1:px);Y=c(:,px+1:px+py);line,l=size(Y);t,q,w,wh,f0,FF=fun717(px,py,c);YY=X*FF(:,2:px+1,3)'+ones(line,1)*FF(:,1,3)'subplot(1,1,1,1)bar(f0(:,:,3)title(' 直方图 ')legend('SG','TZBFB','FHL','