偏最小二乘法算法.docx

1、偏最小二乘法算法偏最小二乘法1.1基本原理偏最小二乘法(PLS)是基于因子分析的多变量校正方法，其数学基础为主成分分析。 但它相对于主成分回归(PCR)更进了一步，两者的区别在于PLS法将浓度矩阵Y和相应 的量测响应矩阵X同时进行主成分分解：X=TP+EY=UQ+F式中T和U分别为X和Y的得分矩阵，而P和Q分别为X和Y的载荷矩阵，E和F 分别为运用偏最小二乘法去拟合矩阵X和Y时所引进的误差。偏最小二乘法和主成分回归很相似，其差别在于用于描述变量Y中因子的同时也用于 描述变量X。为了实现这一点，数学中是以矩阵Y的列去计算矩阵X的因子。同时，矩阵 Y的因子则由矩阵X的列去预测。分解得到的T和U矩阵

2、分别是除去了人部分测量误差的 响应和浓度的信息。偏最小二乘法就是利用各列向量相互正交的特征响应矩阵T和特征浓 度矩阵U进行回归：U=TB得到回归系数矩阵，又称关联矩阵E：B=(TT )F U因此，偏最小二乘法的校正步骤包括对矩阵Y和矩阵X的主成分分解以及对关联矩阵 B的计算。1.2主成分分析主成分分析的中心目的是将数据降维，以排除众多化学信息共存中相互重叠的信息。他 是将原变量进行转换，即把原变量的线性组合成几个新变量。同时这些新变量要尽可能多的 表征原变量的数据结构特征而不丢失信息。新变量是一组正交的，即互不相关的变量。这种 新变量又称为主成分。如何寻找主成分，在数学上讲，求数据矩阵的主成分

3、就是求解该矩阵的特征值和特征矢 量问题。卞面以多组分混合物的量测光谱来加以说明。假设有n个样本包含p个组分，在m 个波长下测定其光谱数据，根据比尔定律和加和定理有：如果混合物只有一种组分，则该光谱矢量与纯光谱矢量应该是方向一致，而人小不同。 换句话说，光谱A表示在由p个波长构成的p维变量空间的一组点(n个)，而这一组点一 定在一条通过坐标原点的直线上。这条直线其实就是纯光谱b。因此由m个波长描述的原始 数据可以用一条直线，即一个新坐标或新变量来表示。如果一个混合物由2个组分组成，各 组分的纯光谱用bl, b2表示，则有：,所以1=1 （=1= Xh.ll2-E（vi）21=1 1=1=勿兀一

4、（叮hj1=1 1=1=ziixiii2-vrx7xvi nun1=1上式等价于vXTXvl max （最大特征值入）上式中5表示第一个主成分轴矢量，即第一个特征矢量，所对应的最大值称为特征值, 用心表示。从上面推导看出，寻找主成分轴就是求X矩阵的协方差矩阵X?X中的最大特征值(Xi)和特征向量(Vi)o下面考虑变量数为m的一般情况。 在m为空间中新变量可以表示为：厂叫】= %】 +片討总+ +叫 d初叫厂冬莎】+冬站2 + +%心I % = %兀】+ %兀+ %為其中系数矩阵V为V11用U和X分别表示新变量和原始矢量，则 If =,x =人2%、u = V x上述m维主成分系数必须满足下面两

5、个条件(1)正交条件：任意两个主成分协、心，其系数的乘积之和为0。(2)归一化条件：对于任一主成分系数的平方和等于1。喑+唁+ +唁=1满足这两个条件的矩阵，称之为正交矩阵。正交矩阵具有如卞性质:V7V = /,V_I = Vr13矩阵的主成分分解根据特征向量和特征值的定义叭=&= 1,2,，加(*)同时令x的协方差矩阵为Z = XTX(*)式两边同时左乘Vi，有主成分系数矩阵V也可写为因此可得ZV = V-diagAi其中diag表示一个对角矩阵，即对角线元素为非对角线元素为0的矩阵。 上式两边同时左乘VL得vrzv =她VTZV = VrXrXV = (XV)TXV =畑令T = XV，则

6、上式变为TtT = 畑人将式T=XV右乘厂1得X = TVt上式是矩阵X的主成分分解的一种表达式，由上式得求解T和V的方法VT = (TtT)1TtXT = XV(VtV)依据矩阵乘法规则即可获得矩阵V和T中每一个矢量的计算公式：vr.=tT.X/tTt,tj = Xvj/v,jtj根据上面两个公式可以设计主成分分解的迭代法算法如卞：(1)取X中任意一列作为起始的to(2)由此 t 计算:vT = tTX/tTt将 Q 归一化：v；H. = v/|v|(4)计算新的 t： t=Xv/vTv(5)比较步骤4所得的t和上一步的t。若二者相等(在给定的误差范围内)，则按(S)计算特征值，转第六步继续

7、进行；否则返回第二步继续迭代。(6)从Y中减去/的贡献：X = X-t-vTo返回1,继续运行，直到最后Y趋近于零。从理论上讲，在m空间中，可以获得m个主成分。但是在实际应用中一般只取前几个 对方差贡献最大的主成分，这样就使高维空间的数据降到低维，如二维或三维空间，非常有 益于数据的观察，同时损失的信息量还不会太大。取前p个主成分的依据为P in比率(尸工工&(=1 1=1般推荐，比率() $80%1.4偏最小二乘法算法(1)矩阵X和Y的标准化处理(2)取Y中任意一列赋给作为起始的u对于X矩阵(3)w,=urX/uIu(4)归一化：W： = w：/|w；|(5)计算新的 t： t=Xw/wrw

8、对于Y矩阵(6)qT=trY/tTt(7)归化：g二.=q；d /1| q； II(8)u=Yq/qTq收敛判据：(9)将步骤8所得的u与前一次迭代的结果相比较，若等于(在允许误差范围内)，到步骤10,否则返回3。(10)pT=tTX/tTt(11)归一化：P二=Pold / | Pold II(12)tnew = told I Pold I(13)w二.=町/|必11计算回归系数b以用于内部关联：(14)b=uTt/tTt对于主成分h计算残差：(15)EiUhP；(16)Fh = Ei-bhZd之后回到步骤(2),去进行下一主成分的运算，直到残差趋近于零。未知样品预测(17)如校正部分，将X

9、矩阵标准化(18)肛0, y=0Q9)h=h+L,q = XW ,y = y + bhgwq； ,x = x- gp；(20)若ha(主成分数)，到步骤(21)o否则返回步骤(19)(21)得到的Y已经标准化，因此需要按标准化步骤的相反操作，将之复原到原坐标注意的 是对预测集进行标准化处理的时，使用的是训练集的均值和标准偏差。因此，在进行反标准 化操作时，使用的也应该是训练集的均值和标准偏差。1.5程序框图与程序代码 程序框图：用C语言编制的PLS程序源代码： C语言的源程序的各函数和功能如下:data_mput();/*数据的输入*/data.standardization Q; pctes

10、tQ principalcomponentQ; factornumQ;nornialeqO；/*数据的标准化*/*在主成分分解时判断t的收敛*/*主成分分解*/*确定必要的主成分分数*/*正规方程的建立*/test(mt k)/*迭代求解系数矩阵的收敛检验*/lterationQ;/*迭代法求解*/predictQ;/*未知样品预测*/biasQ；/*计算预测标准偏差*/repoilQ;/*输出结呆*/inam(argc?aigv)/*主函数*/PLS程序源代码#iiiclude#mclude#iiiclude#mcludedefine N 25define Mil#define L 4def

11、ine TN 8define H 5FILE *mfp,*outfp;double tiain_xNM, tiam_vNL;double test_xTNM, test_yTNL;double avg_xM,std_xM:double avg_yL,std_ yL；double error_valueTN Ldouble uN.new_uN;doublesave_p H N,save_qH L,save_wH M,save_d H:double predictvalueTN L,biasL;void data_mput() /*数据的输入*7 intij;血(j=0；jvMj+)fbr(i=O

12、;iN;i+)fscanf(infA &tiaiii_xij);fbr(i=O;iN;i-H-) foi(j=O;jLj+) fscanf(infA &tiaiii_yij);foi(j=O;jLj+)new_u i+=tmin_y i j *qj; new_ui/=sq; inask=pctest(); while(mask); foi(j=O;jLj+) save_qililj=qlj; fbi(st=O,i=O;iN;i+) st+=ti*ti;pj=O.O;fbr(i=O;iN;i-H-)Pj+=ti*gin_xiIj：pj/=st；fdi(sp=Oj=OjMj+) sp+=pj*pj

13、； sp=sqrt(sp); foi(j=O;jM;j-H-) pj/=sp;save_piliIj=pj; fbr(i=O;iN;i-H-) ti*=sp； fdi(i=O;iM;i+) wi*=sp;save_wih i=wi;fbr(i=O;iTN;i+)fscaiif(mfp, &test_xij);fbr(i=O;iTN;i+) for(j=OjLj+) fscaiif(uifp. u%lfv, &test_yij);fclose(mfp);void data_standardization () /*数据的标准化*/mt ij； av jxj=std_xj=OO; for(i=O;

14、iN;i-H-) av jxj+=tf ain_xi j; avjxj/=N；for(i=O;iN;i-H-)std_xj+=pow(Mam_xiJHvjxj)2); std_xj=sqrt(std_xj)/(N)；for(i=O;iN;i-H-)tram_xij= (tram_xij-avg_xlj)/std-xj;for(i=O;iTN;i+)tram_xij= (tram_xij-avg_xlj)/std-xj; for(j=O;jLj+) avg-yj=std_y|j=OO;fbr(i=O;iN;i-H-)avjyj+=tniin_yi j； avjyj/=N;fbr(i=O;iN;i

15、-H-)std_yj+=pow(tiam_yijMvjyj),2); std_yj=sqrt(std_yj)/(Nl);fbr(i=O;iN;i-H-)for(j=O;jLj+)trauvyi j= (tiain_y 1 (j-avg_yj )/std_y j ; pctestQ /*检验收敛*/ intij;double count=0.0;fbi(i=0;iN;i+) count+=pow(new_u i-u i),2);foi(st=0j=0;iN;i+)st+=ti*ti;foi(b=0a=0:iN;i+)b+=new_ui*ti；b/=st;save_dili=b;for(i=O;i

16、N;i-H-)foi(j=0;jM;j+)elsereturn 1; calibiation() /* 模型的建立*/ mt ij.ilijnask;double su,sw、st，sq.sp,b;double tN.wM.pM.qL;111=0:do for(i=0;iN;i-H-) ui=tram_yil;do for(su=04=0;iN;i-H-) su+=ui*ui; wj=0.0;fbi(i=0;iN;i-H-)w(j+=u 1 *gin_xi j ； wj/=su; fbi(sT=OJ=OjMj+)sw-r=wj*wj； sw=sqrt(sw);foi(j=O;jM;j-H-)w

17、j/=sw;fbi(sT=OJ=OjMj+) sw-r=wj*wj；ti=o.o； for(j=O;jM;j-H-) ti+=train_xi j*w(j； ti/=sw;fbi(st=04=0;iN;i+)st+=ti*ti; foi(j=O;jLj+) qj=O.O；fbr(i=0;iN;i+) qj+=ti*gin_yi|j; intij;fbr(i=0;iTN;i-H-) foi(j=O;jLj+) enoi_valuei j=tesCyi j -predictx-aluei j ; foi(j=O;jLj+) biasj=O.O;fdi(i=0;iTN;i+)biasj/=TN;bi

18、aslj=sqit(biasj);repoitQ /*输出结果*/ mt i,j；fprmtf(outfp,预测的浓度值为：n”);fbr(i=O;iTN;i+)foi(j=O;jLj+)fpiiiitf(outfp, %941f , predictraluei|j);fpriiitf(outfp, );fprmtf(outfp,预测结果的误差为：W ); fbi(i=O;iTN;i+) foi(j=0;jL;j+)fpriiitf(outfp, %94f , enoi_valueij);fpriiitf(outfp, );fpimtf(outfp,相对于各组分的校正标准偏差为n” ); fd

19、i(i=0;iL;i+)rintf(outfp, %1041f” , biasi);void main(aigc,aigr)mt argc；char*argv; if(argc!=3) prmtf(tv forgot enter a filename.Xii); exit(l);if(uifp=fopen(argv 1 , u r )=0 priiitfcant open data filenn);qi=/st;for (sq=0, j=0 JL j+)sq+=qj*qj;sq=sqrt(sq);foi(j=O;jLj+)qj=/sq；for (sq=0, j=0 JL j+)sq+=qj*qlj；for(i=0;iN;i-H-) new_uj=0.0;exit(l);if(outfp=fopen(argv2, uw,T ) = = 0) piintf(caif t open output filenn); exit(l);data_mputQ;data_standaidization(); calibrationQ;predictQ;statistics();repoitQ;