真正的好东西偏最小二乘回归多元线性回归分析+典型相关分析+主成分分析Word格式.docx

1、为了数学推导方便起见 ,首先将数据做标准化处理。op n p据矩阵记为Eo=（Eoi,，Eop）n p，Yj经标准化处理后的数据矩阵记为F0 =（ F0i，F0q ） n P第一步 记ti是Eo的第一个成分，Wi是Eo的第一个轴，它是一个单位向量,既 |w1|=1。记ui是Fo的第一个成分，uFoc，。C1是Fo的第一个轴，并且|ci|=1o如果要ti Ui能分别很好的代表X与丫中的数据变异信息，根据主成分分析原理，应该有Var（u1） maxVar（t1） max另一方面，由于回归建模的需要，又要求ti对ui有很大的解释能力，有典型相关 分析的思路，ti与ui的相关度应达到最大值，既r （t

2、1 u1） max因此，综合起来，在偏最小二乘回归中，我们要求 ti与Ui的协方差达到最大，既Cov（t1U1）=#Var （tJVar （U1）r（t1，uj max正规的数学表述应该是求解下列优化问题，既maxE0W1,F 0C1s.t W1W1 1C1C1 1因此，将在|wi|2=1和|C1|2=1的约束条件下，去求（W E0 F0C1）的最大值。如果采用拉格朗日算法，记2 （C1 C1 1）s=W1 E 0 Fo C1 1 （W1 W1 1）对s分别求关于W1 C1 1和2的偏导并令之为零，有s=-（c1c1 -1）=02由式（1-2）（1-5），可以推出2 1 2 2 wlEoFoq

3、 EoWFoS记1 2 1 2 2 w1EoFoC1 ,所以，1正是优化问题的目标函数值.把式（1-2）和式（1-3）写成将式（1-7）代入式（1-6）,有同理，可得可见,W1是矩阵EoFoFoEo的特征向量，对应的特征值为12. 1是目标函数值，它要求取最大值，所以,W1是对应于EoE。矩阵最大特征值的单位特征向量.而另I I . _ 2 、.一万面，C1是对应于矩阵FoEoEoFo最大特征值1的单位特征向量.求得轴W1和C1后,即可得到成分t1 E0W1U1 F0C1然后,分别求Eo和Fo对t1, U1的三个回归方程Eo1E1（1-10）FoF 1（1-11）t1rF1（1-12）式中，回

4、归系数向量是而E1, F 1, F1分别是三个回归方程的残差矩阵.第二步 用残差矩阵E1和F1取代Eo和Fo，然后，求第二个轴W2和C2以及第 二个成分t2,U2,有t2= E1 w2U2= F1 C22 t2, u2 w 2E 1F1C2C2是对应于矩阵I I QW2是对应于矩阵 E 1F1F 1E1最大特征值 2的特征值,F 1t2|t2|2因此，有回归方程由于,t1, ,tA均可以表示成Eo1, ,Eop的线性组合，因此，式（1-17）还可以还原FAk是残差距阵Fa的第k列。1.3交叉有效性面要讨论的问题是在现有的数据表下 ,如何确定更好的回归方程。在许多情形下 ,偏最小二乘回归方程并不

5、需要选用全部的成分 t1, ,tA 进行回归建模 ,而是 可 以 象 在 主 成 分 分 析 一 样 , 采 用 截 尾 的 方 式 选 择 前 m 个 成 分（m A,A 秩（ X ） ,仅用这 m 个后续的成分 t1, ,tm 就可以得到一个预测性较好 的模型。事实上,如果后续的成分已经不能为解释 F0 提供更有意义的信息时 ,采用过多的成分只会破坏对统计趋势的认识 ,引导错误的预测结论。在多元回归分析一章中 ,我们曾在调整复测定系数的内容中讨论过这一观点。面的问题是怎样来确定所应提取的成分个数。在多元回归分析中 ,曾介绍过用抽样测试法来确定回归模型是否适于预测应用。我们把手中的数据分成两

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

7、代入前面拟合的回归方程，得到yj在 样本点i上的拟合值?hj（ i）。对于每一个i=1,2,n,重复上述测试，则可以定义yj的预测误差平方和为 PRESShj ,有（1-18）n 2PRESShj （yij y?hj（ i）i1定义丫的预测误差平方和为PRESSh，有（1-19）PP RESS P RESSjj 1显然,如果回归方程的稳健性不好，误差就很大,它对样本点的变动就会十分敏感， 这种扰动误差的作用，就会加大PRESSh的值。另外,再采用所有的样本点，拟合含h个成分的回归方程。这是，记第i个样本点的预测值为?hji，则可以记yj的误差平方和为SShj，有SShj（1-20）定义丫的误差

8、平方和为S&，有一般说来,总是有PRESSh大于SSh，而SSh则总是小于SSh 1。下面比较SSh 1和PRESSh 。 SS.1是用全部样本点拟合的具有 h-1个成分的方程的拟合误差PRESSh增加了一个成分th，但却含有样本点的扰动误差。如果h个成分的回归方程的含扰动误差能在一定程度上小于（h-1）个成分回归方程的拟合误差，则认为增 加一个成分th,会使预测结果明显提高。因此我们希望 （PRESSh /SSh 1）的比值能越小越好。在SIMCA-P软件中，指定（P RESSh/SSh 1） 0.952即JP RESS 0.95JSS 1时,增加成分th就是有益的；或者反过来说,当JPRE

9、SS, 0.95JSS7时,就认为增加新的成分th,对减少方程的预测误差无明显的改善作用.另有一种等价的定义称为交叉有效性。对每一个变量 yk,定义对于全部因变量丫,成分th交叉有效性定义为用交叉有效性测量成分th对预测模型精度的边际贡献有如下两个尺度。见，q2 0.0975与（P RESSh/SSh1）0.952是完全等价的决策原则。可以考虑增加成分th是明显有益的。明确了偏最小二乘回归方法的基本原理、 方法及算法步骤后，我们将做 实证分析。fun cti on w=maxdet（A） %求矩阵的最大特征值v,d=eig（A）;n,p =size（d）;d1=d*o nes（ p,1）;d2

10、=max（d1）;i=find（d1=d2）;w=v（:,i）;% function c,m,v=norm1（C） %对数据进行标准化处理n,s=size（C）;for i=1:nfor j=1:sc（i,j）=（C（i,j）-mean（C（:,j）/sqrt（cov（C（:,j）;end endm=mean（C）;v（1,j）=sqrt（cov（C（:end %function t,q,w,wh,f0,FF=fun717（px,py,C）% px自变量的输入个数% py输入因变量的个数。% C输入的自变量和因变量组成的矩阵% t提取的主成分% q为回归系数。% w最大特征值所对应的特征向量。

11、E0=c（:,1:px）;F0=c（:,px+1:px+py）;A=E0*F0*F0*E0;提取主成分t（:,1）=E0*w（:,1）;E（:px）=E0-t（:,1）*（E0*t（:,1）/（t（:,1）,1）;获得回归系数p（:px）=（E0for i=0:px-2B（:,px*i+1:px*i+px）=E（:px*i+px）*E（:px*i+px）w（:,i+2）=maxdet（B（:px*i+px）;% maxdet 为求最大特征值的函数,i+2）=E（:px*i+px）*w（:,i+2）;,px*i+px+1:px*i+2*px）=（E（:,i+2）/（t（:,i+2）*t（ :,i

12、+2）px*i+2*px）=E（:px*i+px）-t（:,i+2）*（E（:px* i+px）endfor s=1:px求回归系数%noq（:,s）=p（1,px*（s-1）+1:px*s）n,d=size（q）;for h=1:iw=eye（d）;h-1iw=iw*（eye（d）-w（:,j）*q（:,j）;end wh（:,h）=iw*w（:,h）;pyzr（j,:）=（regress1（y（:,j）,t）fori=1:生成标准化变量的方程的系数矩阵w1=wh（:j）;zr1=（zr（i,1:j）f0（i,:,j）=（w1*zr1）normxy,meanxy,covxy=norm1（C）

13、;rmxy 标准化后的数据矩阵%mea nx每一列的均值%covxy每一列的方差ccxx=ones（py,1）*meanxy（1,1: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:FF（:,j）=fff,ff;原始变量方程的常数项和系数矩阵%function r,Rdyt,RdYt,RdYtt,Rdytt,VIP=fun8y（px,py,c）X=c（:Y=c（:x

14、=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:RdYtt（1,m）=sum（RdYt（1,1:m）for m=1:Rdytt（j,m）=sum（Rdyt（j,1:Rd（j,m）=RdYt（1,1:m）*（w（j,1:m）A2）end for j=1:VIP（j,:）=sqrt（px*ones（1,px）./RdYtt）.*Rd（j,:）;function r,Rdxt,RdXt,RdXtt,Rdxtt=fun

15、8x（px,py,c）r1=corrcoef（x,t）;r=r1（px+1:px+px,1:px）Rdxt=r.A2;RdXt=mean（Rdxt）;RdXtt（1,m）=sum（RdXt（1,1:Rdxtt（j,m）=sum（Rdxt（j,1:% for j=1:Rd（j,m）=RdXt（1,1:function t,u=TU（px,py,C）%t 提取的自变量的主成分% u 提取的因变量的主成分c=norm1（C）;y=c（:,1）=maxdet（A）;B=F0*E0*E0*F0;cc（:,1）=maxdet（B）;u（:,1）=F0*cc（:function drew（px,py,c）l

16、ine,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）legend（SG,TZBFBFHLJKHPZDJPZDTZZGGPK）grid onplot（YY（:,4）,Y（:,4）,+lslinev=mod（i,4）;d=（i-v）/4;subplot（2,2,v,d+1） plot（YY（:,i）,Y（:,i）,*%function Qhj,Qh,prey=crossval7（px,py,c）%px 是自变量的个数；%py

17、 是因量PRESShj=zeros（px,py）;line,row=size（x）;linenewx=X;newy=Y;newx（j,:）=;newy（j,:t,p0,w,wh,f0,FF=fun717（px,py,newx,newy）;prey（j,:,h）=X（j,:）*FF（:px+1,h）+FF（:,1,h）P RESShj（h,:）=sum（Y- prey（:，:,6）八2）;PRESSh=（sum（PRESShj）t1,p0,w,wh,f0,FF=fun717（px,py,c）;prey2（:,h）=X（:+ones（line,1）*FF（:SShj（h,:）=sum（Y- prey2（:,h）A2）;SSh=（sum（SShjQhj=ones（px-1,py）-PRESShj（2:px,:）./SShj（1:px -1,: % 错位Qh=ones（px-1,1）-PRESSh（2:px,1）./SSh（1:px-1,1）;