pls代码偏最小二乘

PLS算法1偏最小二乘回归是一种新型的多元统计数据分析方法，它与1983年由伍德和阿巴诺等人首次提出。近十年来，它在理论、方法和应用方面都得到了迅速的发展。密西根大学的弗耐尔教授称偏最小二乘回归为第二代回归分析方法。偏最小二乘回归方法在统计应用中的重要性主要的有以下几个方面：（1）偏最小二乘回归是一种多因变量对多自变量的回归建模方法。（2）偏最小二乘回归可以较好地解决许多以往用普通多元回归无法解决的问题。在普通多元线形回归的应用中，我们常受到许多限制。最典型的问题就是自变量之间的多重相关性。如果采用普通的最小二乘方法，这种变量多重相关性就会严重危害参数估计，扩大模型误差，并破坏模型的稳定性。变量多重相关问题十分复杂，长期以来在理论和方法上都未给出满意的答案，这一直困扰着从事实际系统分析的工作人员。在偏最小二乘回归中开辟了一种有效的技术途径，它利用对系统中的数据信息进行分解和筛选的方式，提取对因变量的解释性最强的综合变量，辨识系统中的信息与噪声，从而更好地克服变量多重相关性在系统建模中的不良作用。（3）偏最小二乘回归之所以被称为第二代回归方法，还由于它可以实现多种数据分析方法的综合应用。偏最小二乘回归=多元线性回归分析+典型相关分析+主成分分析由于偏最小二乘回归在建模的同时实现了数据结构的简化，因此，可以在二维平面图上对多维数据的特性进行观察，这使得偏最小二乘回归分析的图形功能十分强大。在一次偏最小二乘回归分析计算后，不但可以得到多因变量对多自变量的回归模型，而且可以在平面图上直接观察两组变量之间的相关关系，以及观察样本点间的相似性结构。这种高维数据多个层面的可视见性，可以使数据系统的分析内容更加丰富，同时又可以对所建立的回归模型给予许多更详细深入的实际解释。一、偏最小二乘回归的建模策略\原理\方法1.1建模原理PLS算法2设有q个因变量{qyy,...,1}和p自变量{pxx,...,1}。为了研究因变量和自变量的统计关系,我们观测了n个样本点,由此构成了自变量与因变量的数据表X={pxx,...,1}和.Y={qyy,...,1}。偏最小二乘回归分别在X与Y中提取出成分1t和1u(也就是说,1t是pxx,...,1的线形组合,1u是qyy,...,1的线形组合).在提取这两个成分时,为了回归分析的需要,有下列两个要求:(1)1t和1u应尽可能大地携带他们各自数据表中的变异信息;(2)1t与1u的相关程度能够达到最大。这两个要求表明，1t和1u应尽可能好的代表数据表X和Y,同时自变量的成分1t对因变量的成分1u又有最强的解释能力。在第一个成分1t和1u被提取后，偏最小二乘回归分别实施X对1t的回归以及Y对1u的回归。如果回归方程已经达到满意的精度，则算法终止；否则,将利用X被1t解释后的残余信息以及Y被1t解释后的残余信息进行第二轮的成分提取。如此往复，直到能达到一个较满意的精度为止。若最终对X共提取了m个成分1t，…，mt，偏最小二乘回归将通过实施ky对1t，…，mt，的回归,然后再表达成ky关于原变量x1，…，xm，的回归方程,k=1,2,…,q。1.2计算方法推导为了数学推导方便起见,首先将数据做标准化处理。X经标准化处理后的数据矩阵记为0E=(E01，…，Ep0)pn，jY经标准化处理后的数据矩阵记为0F=(01F，…，qF0)pn。第一步记1t是0E的第一个成分，1w是0E的第一个轴，它是一个单位向量，既||1w||=1。记1u是0F的第一个成分，1u=0Fc1。c1是0F的第一个轴，并且||c1||=1。如果要1t，1u能分别很好的代表X与Y中的数据变异信息，根据主成分分PLS算法3析原理，应该有Var(1u)maxVar(1t)max另一方面，由于回归建模的需要，又要求1t对1u有很大的解释能力，有典型相关分析的思路，1t与1u的相关度应达到最大值，既r（1t，1u）max因此，综合起来，在偏最小二乘回归中，我们要求1t与1u的协方差达到最大，既Cov(1t，1u)=)()(11utVarVarr(1t，1u)max正规的数学表述应该是求解下列优化问题，既cw11,maxcFwE1010,s.t111'11'1ccww因此，将在||1w||2=1和||c1||2=1的约束条件下，去求(w'1E'0F0c1)的最大值。如果采用拉格朗日算法，记s=w'1E'00Fc1－1(w'11w－1)－2(c'1c1－1)对s分别求关于1w，c1，1和2的偏导并令之为零，有1ws=E'00Fc1－121w=0(1-2)1cs=F'00E1w－22c1=0(1-3)1s=－(w'11w－1)=0(1-4)2s=－(c'1c1－1)=0(1-5)PLS算法4由式(1-2)~(1-5),可以推出1010100'1'21,22cFwEcFEw记100'1'21122cFEw,所以,1正是优化问题的目标函数值.把式(1-2)和式(1-3)写成11100'wcFE(1-6)11100'cwEF(1-7)将式(1-7)代入式(1-6),有121100'00'wwEFFE(1-8)同理,可得121100'00'ccFEEF(1-9)可见,1w是矩阵00'00'EFFE的特征向量,对应的特征值为21.1是目标函数值,它要求取最大值,所以,1w是对应于00'00'EFFE矩阵最大特征值的单位特征向量.而另一方面,1c是对应于矩阵00'00'FEEF最大特征值21的单位特征向量.求得轴1w和1c后,即可得到成分101wEt101cFu然后,分别求0E和0F对1t,1u的三个回归方程11'10EptE(1-10)11'10FquF(1-11)11'10FrtF(1-12)式中,回归系数向量是2110'1||||ttEp(1-13)PLS算法52110'1||||uuFq(1-14)2110'1||||ttFr(1-15)而1E,1F,1F分别是三个回归方程的残差矩阵.第二步用残差矩阵1E和1F取代0E和0F,然后,求第二个轴2w和2c以及第二个成分2t,2u,有2t=1E2w2u=1F2c211'2'222,cFEwut2w是对应于矩阵11'11'EFFE最大特征值22的特征值,2c是对应于矩阵11'11'FEEF最大特征值的特征向量.计算回归系数2221'2||||ttEp2221'2||||ttFr因此,有回归方程22'21EptE22'21FrtF如此计算下去,如果X的秩是A,则会有AAptptE'1'10（1-16）AAAFrtrtF'1'10（1-17）由于,Att,,1均可以表示成pEE001,,的线性组合,因此,式(1-17)还可以还原成kkFy0*关于kjEx0*的回归方程形式，即PLS算法6AkpkpkkFxxy**11*k=1,2,…,qAkF是残差距阵AF的第k列。1.3交叉有效性下面要讨论的问题是在现有的数据表下,如何确定更好的回归方程。在许多情形下,偏最小二乘回归方程并不需要选用全部的成分Att,,1进行回归建模,而是可以象在主成分分析一样,采用截尾的方式选择前m个成分))(,(XAAm秩,仅用这m个后续的成分mtt,,1就可以得到一个预测性较好的模型。事实上,如果后续的成分已经不能为解释0F提供更有意义的信息时,采用过多的成分只会破坏对统计趋势的认识,引导错误的预测结论。在多元回归分析一章中,我们曾在调整复测定系数的内容中讨论过这一观点。下面的问题是怎样来确定所应提取的成分个数。在多元回归分析中,曾介绍过用抽样测试法来确定回归模型是否适于预测应用。我们把手中的数据分成两部分:第一部分用于建立回归方程,求出回归系数估计量Bb,拟合值Byˆ以及残差均方和2ˆB;再用第二部分数据作为实验点,代入刚才所求得的回归方程,由此求出2ˆˆTTy和。一般地,若有2ˆT2ˆB,则回归方程会有更好的预测效果。若2ˆT2ˆB,则回归方程不宜用于预测。在偏最小二乘回归建模中,究竟应该选取多少个成分为宜,这可通过考察增加一个新的成分后,能否对模型的预测功能有明显的改进来考虑。采用类似于抽样测试法的工作方式,把所有n个样本点分成两部分:第一部分除去某个样本点i的所有样本点集合(共含n-1个样本点),用这部分样本点并使用h个成分拟合一个回归方程;第二部分是把刚才被排除的样本点i代入前面拟合的回归方程,得到jy在样本点i上的拟合值)(ˆihjy。对于每一个i=1,2,…,n,重复上述测试,则可以定义jy的预测误差平方和为hjPRESS,有niihjijhjyyPRESS12)()ˆ((1-18)定义Y的预测误差平方和为hPRESS,有PLS算法7pjhjhPRESSPRESS1(1-19)显然,如果回归方程的稳健性不好,误差就很大,它对样本点的变动就会十分敏感,这种扰动误差的作用,就会加大hPRESS的值。另外,再采用所有的样本点,拟合含h个成分的回归方程。这是,记第i个样本点的预测值为hjiyˆ,则可以记jy的误差平方和为hjSS,有nihjiijhjyySS12)ˆ((1-20)定义Y的误差平方和为hSS,有pjhjhSSSS1(1-21)一般说来,总是有hPRESS大于hSS,而hSS则总是小于1hSS。下面比较1hSS和hPRESS。1hSS是用全部样本点拟合的具有h-1个成分的方程的拟合误差;hPRESS增加了一个成分ht,但却含有样本点的扰动误差。如果h个成分的回归方程的含扰动误差能在一定程度上小于(h-1)个成分回归方程的拟合误差,则认为增加一个成分ht,会使预测结果明显提高。因此我们希望)/(1hhSSPRESS的比值能越小越好。在SIMCA-P软件中,指定2195.0)/(hhSSPRESS即195.0hhSSPRESS时,增加成分ht就是有益的;或者反过来说,当195.0hhSSPRESS时,就认为增加新的成分ht,对减少方程的预测误差无明显的改善作用.另有一种等价的定义称为交叉有效性。对每一个变量ky,定义khhkhkSSPRESSQ)1(21(1-22)对于全部因变量Y,成分ht交叉有效性定义为PLS算法8)1()1(1211hhkhqkhkhSSPRESSSSPRESSQ(1-23)用交叉有效性测量成分ht对预测模型精度的边际贡献有如下两个尺度。(1)当0975.0)95.01(22hQ时,ht成分的边际贡献是显著的。显而易见,0975.02hQ与2195.0)/(hhSSPRESS是完全等价的决策原则。(2)对于k=1,2,…,q,至少有一个k,使得0975.02hQ这时增加成分ht,至少使一个因变量ky的预测模型得到显著的改善,因此,也可以考虑增加成分ht是明显有益的。明确了偏最小二乘回归方法的基本原理、方法及算法步骤后，我们将做实证分析。附录functionw=maxdet(A)%求矩阵的最大特征值[v,d]=eig(A);[n,p]=size(d);d1=d*ones(p,1);d2=max(d1);i=find(d1==d2);PLS算法9w=v(:,i);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[c,m,v]=norm1(C)%对数据进行标准化处理[n,s]=size(C);fori=1:nforj=1:sc(i,j)=(C(i,j)-mean(C(:,j)))/sqrt(cov(C(:,j)));endendm=mean(C);forj=1:sv(1,j)=sqrt(cov(C(:,j)));end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[t,q,w,wh,f0,FF]=fun717(px,p