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
2、 的因子则由矩阵 X 的列去预测。分解得到的 T 和 U 矩阵分别是除去了大部分测量误差的响应和浓度的信息。偏最小二乘法就是利用各列向量相互正交的特征响应矩阵 T 和特征浓度矩阵 U 进行回归：U=TB得到回归系数矩阵，又称关联矩阵 B：B=(TTT-1)TTU因此，偏最小二乘法的校正步骤包括对矩阵 Y 和矩阵 X 的主成分分解以及对关联矩阵B 的计算。1.2 主成分分析主成分分析的中心目的是将数据降维，以排除众多化学信息共存中相互重叠的信息。他是将原变量进行转换，即把原变量的线性组合成几个新变量。同时这些新变量要尽可能多的表征原变量的数据结构特征而不丢失信息。新变量是一组正交的，即互不相关的
3、变量。这种新变量又称为主成分。如何寻找主成分，在数学上讲，求数据矩阵的主成分就是求解该矩阵的特征值和特征矢量问题。下面以多组分混合物的量测光谱来加以说明。假设有 n 个样本包含 p 个组分，在 m 个波长下测定其光谱数据，根据比尔定律和加和定理有：Anm=CnpBpm如果混合物只有一种组分，则该光谱矢量与纯光谱矢量应该是方向一致，而大小不同。换句话说，光谱 A 表示在由 p 个波长构成的 p 维变量空间的一组点（n 个） ，而这一组点一定在一条通过坐标原点的直线上。这条直线其实就是纯光谱 b。因此由 m 个波长描述的原始数据可以用一条直线，即一个新坐标或新变量来表示。如果一个混合物由 2 个组
4、分组成，各组分的纯光谱用 b1，b2 表示，则有： 12TTiiiacb有上式看出，不管混合物如何变化，其光谱总可以用两个新坐标轴 b1,b2 来表示。因此可以推出，如果混合物由 p 个组分组成，那么混合物的光谱就可由 p 个主成分轴的线性组合表示。因而现在的问题就变成了如何求解这些主成分轴。而寻找这些坐标轴的基本原则是使新坐标轴包含原数据的最大方差。即沿着新坐标轴的方向，使方差达到最大。而其他方向，使方差达到最小。从几何角度看，就是变量空间中所有的点到这个新坐标轴的距离最短。以二维空间的为例说明如何寻找主成分坐标轴。变量空间的每一个数据点（一个样本）都可以用通过该点与坐标原点的一个矢量 xi
5、 表征。 主 成 分 轴2x1v1 byic a0上图中直角三角形的三个边长分别以 a,b,c 表示，那么这 n 个点到第一个主成分轴 v1距离的平方和可以通过勾股定理与矢量点积得出： 2211()nniiiiibca因为 与 ，所以2|iicx1|cosTiivx222111(|)nniiib2211|()nnTiiiixv2111|()nnTiiiixmin211|nTiixvX上式等价于max （最大特征值 ）1Tv上式中 v1 表示第一个主成分轴矢量，即第一个特征矢量，所对应的最大值称为特征值，用 1 表示。从上面推导看出，寻找主成分轴就是求 X 矩阵的协方差矩阵 XTX 中的最大特征
6、值（ i）和特征向量（v i） 。下面考虑变量数为 m 的一般情况。在 m 为空间中新变量可以表示为： 11212 212iiimiimimimiuvxvx 其中系数矩阵 V 为V=12112mmmvv 用 u 和 x 分别表示新变量和原始矢量，则，12iimu12iimxuVx上述 m 维主成分系数必须满足下面两个条件（1） 正交条件：任意两个主成分 uk、 ur，其系数的乘积之和为 0。12rkmrvv（2）归一化条件：对于任一主成分系数的平方和等于 1。221kk满足这两个条件的矩阵，称之为正交矩阵。正交矩阵具有如下性质： 1,TTVI1.3 矩阵的主成分分解根据特征向量和特征值的定义(
7、*),12,TiiivXm同时令 X 的协方差矩阵为 TZX(*)式两边同时左乘 vi，有,12,iiZvm主成分系数矩阵 V 也可写为 12(,)Vv因此可得 iZdag其中 表示一个对角矩阵，即对角线元素为 ，非对角线元素为 0 的矩阵。idagi上式两边同时左乘 VT，得 TiVZdag()TTiXV令 ,则上式变为 将式 右乘 得TXVidag1T上式是矩阵 X 的主成分分解的一种表达式，由上式得求解 T 和 V 的方法1()TVXT依据矩阵乘法规则即可获得矩阵 V 和 T 中每一个矢量的计算公式： /,/Tjjjjjvttvt根据上面两个公式可以设计主成分分解的迭代法算法如下：(1)
8、 取 X 中任意一列作为起始的 t。(2) 由此 t 计算： /TTvtX(3) 将 vT 归一化： newoldlv(4) 计算新的 t： /T(5) 比较步骤 4 所得的 t 和上一步的 t。若二者相等（在给定的误差范围内） ，则按（）计算特征值，转第六步继续进行；否则返回第二步继续迭代。Tt(6) 从 Y 中减去 的贡献： 。返回 1，继续运行，直到最后 Y 趋近于零。Ttv TXtv从理论上讲，在 m 空间中，可以获得 m 个主成分。但是在实际应用中一般只取前几个对方差贡献最大的主成分，这样就使高维空间的数据降到低维，如二维或三维空间，非常有益于数据的观察，同时损失的信息量还不会太大。
9、取前 p 个主成分的依据为比率(%) 1/pmii一般推荐，比率(%)80%1.4 偏最小二乘法算法(1) 矩阵 X 和 Y 的标准化处理(2) 取 Y 中任意一列赋给作为起始的 u对于 X 矩阵(3) wT=uTX/uTu(4) 归一化： /|Tnewoldl(5) 计算新的 t：t=Xw/w Tw对于 Y 矩阵(6) qT=tTY/tTt(7) 归一化： /|Tnewoldlq(8) u=Yq/qTq收敛判据：(9) 将步骤 8 所得的 u 与前一次迭代的结果相比较，若等于（在允许误差范围内）,到步骤10，否则返回 3。(10) pT=tTX/tTt(11) 归一化： /|Tnewoldl
10、p(12) tnew = told | pold|(13) /|TTneoldl计算回归系数 b 以用于内部关联：(14) b=uTt/tTt 对于主成分 h 计算残差：(15) 1hhEp(16) TFbtwq之后回到步骤(2),去进行下一主成分的运算，直到残差趋近于零。未知样品预测(17) 如校正部分，将 X 矩阵标准化(18) h=0，y=0(19) h=h+1, , ,ThtWThybtwqThxtp(20) 若 ha(主成分数)，到步骤(21)。否则返回步骤(19)(21) 得到的 Y 已经标准化，因此需要按标准化步骤的相反操作，将之复原到原坐标注意的是对预测集进行标准化处理的时，使
11、用的是训练集的均值和标准偏差。因此，在进行反标准化操作时，使用的也应该是训练集的均值和标准偏差。1.5 程序框图与程序代码程序框图：输入数据矩阵：训练集 X，Y 和预测集 X数据标准化：Xij=(Xij-Xj)/Sj从 Y 中任选一列作为初始 u， h=h+1wT=uTX/uTu, /|Tnewoldlt=Xw/wTw,qT=tTY/tTt, , u=Yq/qTq/|neoldlu=Yq/qTqqT=tTY/tTtt=Xw/wTw两次迭代 t 相等吗pT=tTX/tTt， ，/|Tneoldlptnew = told | pold|， ， /|Tnewoldlpb=uTt/tTt， ，1Thh
12、Et1hFbqh#include#include#include#define N 25#define M 11#define L 4#define TN 8#define H 5FILE *infp,*outfp;double train_xNM, train_yNL;double test_xTNM, test_yTNL;double avg_xM,std_xM；double avg_yL,std_ yL；double error_valueTNLdouble uN,new_uN;double save_pHN,save_qHL,save_wHM,save_dH;double predic
13、tvalueTNL,biasL;void data_input() /*数据的输入*/ int i,j;for(j=0;jM;j+)for(i=0;iN;i+)fscanf(infp,“%lf”，for(i=0;iN;i+)for(j=0;jL;j+)fscanf(infp,“%lf”，for(j=0;jM;j+)for(j=0;jL;j+)new_ui+=train_yij*qj;new_ui/=sq;mask=pctest();while(mask);for(j=0;jL;j+)save_qihj=qj;for(st=0,i=0;iN;i+)st+=ti*ti;for(j=0;jM;j+)
14、pj=0.0;for(i=0;iN;i+)pj+=ti*train_xij；pj/=st;for(sp=0,j=0;jM;j+)sp+=pj* pj;sp=sqrt(sp);for(j=0;jM;j+) pj/=sp;save_pihj=pj;for(i=0;iN;i+)ti*=sp;for(i=0;iM;i+) wi*=sp;save_wihi=wi;for(i=0;iTN;i+)fscanf(infp,“%lf”，for(i=0;iTN;i+)for(j=0;jL;j+)fscanf(infp,“%lf”，fclose(infp);void data_standardization ()
15、/*数据的标准化*/int i,j;for(j=0;jM;j+) avg_xj=std_xj=0.0;for(i=0;iN;i+)avg_xj+=train_xij;avg_xj/=N；for(i=0;iN;i+)std_xj+=pow(train_xij-avg_xj),2);std_xj=sqrt(std_xj)/(N-1);for(i=0;iN;i+)for(j=0;jM;j+)train_xij= (train_xij-avg_xj)/std_xj;for(i=0;iTN;i+)for(j=0;jM;j+)train_xij= (train_xij-avg_xj)/std_xj;for
16、(j=0;jL;j+) avg_yj=std_yj=0.0;for(i=0;iN;i+)avg_yj+=train_yij;avg_yj/=N;for(i=0;iN;i+)std_yj+=pow(train_yij-avg_yj),2);std_yj=sqrt(std_yj)/(N-1);for(i=0;iN;i+)for(j=0;jL;j+)train_yij= (train_yij-avg_yj)/std_yj;pctest() /*检验收敛*/ int i,j;double count=0.0;for(i=0;iN;i+) count+=pow(new_ui-ui),2);for(st=
17、0,i=0;iN;i+)st+=ti*ti;for(b=0,i=0;iN;i+)b+=new_ui*ti；b/=st;save_dih=b;for(i=0;iN;i+)for(j=0;jM;j+)train_xij-=ti*pj;for(i=0;iN;i+)for(j=0;jL;j+)train_yij-=b*ti*qj;ih+;while(ihH);predict() /*计算校正集的浓度*/int i,j,k;double ptTNfor(i=0;iTN;i+)for(j=0;jL;j+)predictvalueij=0.0;for(k=0;kH;k+) for(i=0;iTN;i+) p
18、ti=0.0;for(j=0;jM;j+)pti+=test_xij*save_wkj;for(i=0;iTN;i+)for(j=0;jL;j+)predictvalueij+=save_dk*pti*save_qkj;for(i=0;iTN;i+) for(j=0;jM;j+)test_xij-=pti*save_pkj;for(i=0;iTN;i+) for(j=0;jM;j+)test_xij-=pti*save_pkj;for(i=0;iTN;i+) for(j=0;jL;j+) predictvalueij*=std_yi;predictvalueij+=avg_yj;statist
19、ics() /*计算结果评估 */ui=new_ui;count=sqrt(count);if(count1e-12)return 0;else return 1;calibration() /*模型的建立*/ int i,j,ih,mask;double su,sw,st,sq,sp,b;double tN,wM,pM,qL;ih=0;do for(i=0;iN;i+)ui=train_yil;do for(su=0,i=0;iN;i+)su+=ui*ui;for(j=0;jM;j+) wj=0.0;for(i=0;iN;i+)wj+=ui*train_xij;wj/=su;for(sw=0
20、,j=0;jM;j+)sw+=wj*wj;sw=sqrt(sw);for(j=0;jM;j+)wj/=sw;for(sw=0,j=0;jM;j+)sw+=wj*wj;for(i=0;iN;i+) ti=0.0;for(j=0;jM;j+)ti+=train_xij*wj;ti/=sw;for(st=0,i=0;iN;i+)st+=ti*ti;for(j=0;jL;j+) qj=0.0;for(i=0;iN;i+)qj+=ti*train_yij; int i,j;for(i=0;iTN;i+)for(j=0;jL;j+)error_valueij=test_yij-predictvalueij
21、;for(j=0;jL;j+) biasj=0.0;for(i=0;iTN;i+)biasj/=TN;biasj=sqrt(biasj);report() /*输出结果*/ int i,j;fprintf(outfp,“预测的浓度值为：n ”);for(i=0;iTN;i+)for(j=0;jL;j+)fprintf(outfp, “%9.4lf”, predictvalueij);fprintf(outfp,“n”);fprintf(outfp,“预测结果的误差为：n ”);for(i=0;iTN;i+)for(j=0;jL;j+)fprintf(outfp, “%9.4lf”, error
22、_valueij);fprintf(outfp,“n”);fprintf(outfp,“相对于各组分的校正标准偏差为n ”);for(i=0;iL;i+)fprintf(outfp, “%10.4lf”, biasi);void main(argc,argv) int argc；char*argv; if(argc!=3) printf(“youv forgot enter a filename.n”);exit(1);if(infp=fopen(argv1,“r”)=0 printf(“cant open data filen”);qi=/st;for (sq=0, j=0;jL;j+)sq+=qj*qj;sq=sqrt(sq);for(j=0;jL;j+)qj=/sq;for (sq=0, j=0;jL;j+)sq+=qj*qj;for(i=0;iN;i+) new_uj=0.0;exit(1);if(outfp=fopen(argv2,“w”) = = 0) printf(“cant open output filen” );exit(1);data_input();data_standardization();calibration();predict(); statistics(); report();
