1、偏最小二乘方法偏最小二乘方法1第六章第六章 偏最小二乘方法偏最小二乘方法 偏最小二乘方法偏最小二乘方法(PLS-Partial Least Squares)是近年来发展是近年来发展起来的一种新的多元统计分析法起来的一种新的多元统计分析法,现已成功地应用于分析化学现已成功地应用于分析化学,如紫外光谱、气相色谱和电分析化学等等。该种方法,在化合如紫外光谱、气相色谱和电分析化学等等。该种方法,在化合物结构物结构-活性活性/性质相关性研究中是一种非常有用的手段。如美国性质相关性研究中是一种非常有用的手段。如美国TriposCoMFA(Comparative Molecular Field Analys
2、is)方法方法,其其中,数据统计处理部分主要是中,数据统计处理部分主要是PLS。在在PLS方法中用的是替潜变方法中用的是替潜变量,其数学基础是主成分分析。替潜变量的个数一般少于原自量,其数学基础是主成分分析。替潜变量的个数一般少于原自变量的个数,所以变量的个数,所以PLS特别适用于自变量的个数多于试样个数特别适用于自变量的个数多于试样个数的情况。在此种情况下,亦可运用主成分回归方法,但不能够的情况。在此种情况下,亦可运用主成分回归方法,但不能够运用一般的多元回归分析,因为一般多元回归分析要求试样的运用一般的多元回归分析,因为一般多元回归分析要求试样的个数必须多于自变量的个数。个数必须多于自变量
3、的个数。偏最小二乘方法偏最小二乘方法2 6.1 多元线性回归(多元线性回归(MLR)若自变量为若自变量为m个,个,xj(j=1,2,m),因变量为因变量为y,在在y与与xj间,间,我们可以建立一线性模型,即我们可以建立一线性模型,即exbxbxbymm.2211(6.1a)exbyjmjj1(6.1b)ebxy(6.1c)在式中,在式中,bj为回归系数。为回归系数。在式(在式(6.1)中仅有一个试样,若有)中仅有一个试样,若有n个试样,即个试样,即为为yi(i=1,2,n),它的列向量形式为它的列向量形式为y,b与原来相同,与原来相同,矢量矢量xj为矩阵为矩阵X的行,则:的行,则:y=Xb+e
4、偏最小二乘方法偏最小二乘方法3若用图形表示,则为:若用图形表示,则为:y =X B +e 1m11 n n nm在此情况下,在此情况下,n为试样数,为试样数,m为自变量数。有如下三种情况:为自变量数。有如下三种情况:(1)mn,即变量数多于试样数,对于即变量数多于试样数,对于b来说,则有无穷多个解。来说,则有无穷多个解。(2)m=n,变量数与试样数相等,若矩阵变量数与试样数相等,若矩阵X满秩时,则矢量满秩时,则矢量b有有唯一解。但是,在实际工作中,这种情况是极少能碰到的。唯一解。但是,在实际工作中,这种情况是极少能碰到的。此时我们有:此时我们有:e=y Xb=0偏最小二乘方法偏最小二乘方法4
5、(3)m (主成分数)到步主成分数)到步(5)(5),否则到步,否则到步(3)(3)。(5)(5)得到的得到的Y Y为已经标准化,因此需按照标准化步骤的相为已经标准化,因此需按照标准化步骤的相反操作,将之恢复到原始坐标。反操作,将之恢复到原始坐标。4.4.关于主成分数关于主成分数 若若X X和和Y Y间关系符合线性模型,则描述模型的主成分数应间关系符合线性模型,则描述模型的主成分数应与模型的维数相等。主成分数是偏最小二乘模型的重要性质。与模型的维数相等。主成分数是偏最小二乘模型的重要性质。由于测试数据一般隐含噪声,故主成分数通常与由于测试数据一般隐含噪声,故主成分数通常与X X的秩不的秩不相等
6、。如前已述及,在实际问题的处理中,总是要消去一些相等。如前已述及,在实际问题的处理中,总是要消去一些因子(成分),因子(成分),因为这些因子所表征的主要是测试误差、噪因为这些因子所表征的主要是测试误差、噪声及由于变量间相关所引起的共线问题等。声及由于变量间相关所引起的共线问题等。偏最小二乘方法偏最小二乘方法32 确定主成分数的一种方法是以式(确定主成分数的一种方法是以式(6.86.8)中)中F Fh h的模数为判的模数为判据。图据。图6.36.3为模数为模数 对主成分数所得关系曲线,可以选定某对主成分数所得关系曲线,可以选定某值作为门限,当值作为门限,当 小于此值时,则停止迭代。小于此值时,则
7、停止迭代。hFhF图图6.3 6.3 与偏最小二乘中因子书的关系与偏最小二乘中因子书的关系hF 另一种方法是运用另一种方法是运用F F 检验来测试内在相关性检验来测试内在相关性(inner inner relation)relation)以确证所建立的模型。以确证所建立的模型。偏最小二乘方法偏最小二乘方法33 再一种方法为交叉验证法。再一种方法为交叉验证法。在这种方法中计算一统计量在这种方法中计算一统计量PRESSPRESS(prediction residual sum of prediction residual sum of squaressquares),),即预测残差之平方和。即预测
8、残差之平方和。如图如图6.46.4所示,显然,人们总是希所示,显然,人们总是希望采用某一主成分数时所产生的望采用某一主成分数时所产生的PRESSPRESS为最小。但最小的位置常难为最小。但最小的位置常难以准确确定。用这种方法确定主成以准确确定。用这种方法确定主成数非常类似于测定下限的概念。所数非常类似于测定下限的概念。所谓测定下限即在噪声存在下最小可谓测定下限即在噪声存在下最小可以检出的信号。在图以检出的信号。在图6.46.4的情况下,的情况下,因子数可取因子数可取4 48 8。图图6.4 6.4 PRESSPRESS与因子数的关系与因子数的关系偏最小二乘方法偏最小二乘方法34 5.5.应用实
9、例应用实例腐植酸和木质磺酸盐的荧光分光光度分析腐植酸和木质磺酸盐的荧光分光光度分析55 磺酸木质素磺酸木质素(ligninsulfonate)ligninsulfonate)是水中的是水中的一种污染物,可用荧光分光一种污染物,可用荧光分光光度法测定光度法测定.尽管此种方法具尽管此种方法具有高灵敏度和高选择性,但有高灵敏度和高选择性,但在磺酸木质素的测试中腐植在磺酸木质素的测试中腐植酸和去污剂中的光白剂酸和去污剂中的光白剂(optical whiteneroptical whitener)对其对其严重干扰。这三种化合物的严重干扰。这三种化合物的发射光谱重叠非常严重发射光谱重叠非常严重(见图见图6
10、.5).6.5).由图可见,没有一个区由图可见,没有一个区域仅为一种化合物所具有的域仅为一种化合物所具有的发射光谱发射光谱.图图6.5 6.5 腐植酸腐植酸(),),磺酸木质素磺酸木质素(-)(-)和去污剂和去污剂()的发射光的发射光谱谱(均由纯物质测试所得均由纯物质测试所得)偏最小二乘方法偏最小二乘方法35这三种化合物不仅发射光谱严这三种化合物不仅发射光谱严重重叠重重叠,同时在溶液中相互间有同时在溶液中相互间有影响影响,如图如图6.66.6所示所示,三种纯物质三种纯物质的发射光谱加和的发射光谱加和()与其混合与其混合溶液的发射光谱溶液的发射光谱()并不一并不一样样,这就进一步增加了问题的复这
11、就进一步增加了问题的复杂性杂性.但是借助于偏最小二乘法但是借助于偏最小二乘法,可以进行单一成分的测试可以进行单一成分的测试,所得所得结果尚较满意结果尚较满意.图图6.6 6.6 腐植酸腐植酸,磺酸木质素和去污磺酸木质素和去污剂纯溶液发射光谱加剂纯溶液发射光谱加和和()及三物质混合溶液的发射及三物质混合溶液的发射光谱光谱(-)(-)偏最小二乘方法偏最小二乘方法36 首先首先,看一下二组分的情况看一下二组分的情况,表表6.16.1所示为腐植酸和磺酸木所示为腐植酸和磺酸木质素混合样品的浓度测定结果。质素混合样品的浓度测定结果。表表6.1 腐植酸与磺酸木质素混合物溶液测试结果腐植酸与磺酸木质素混合物溶
12、液测试结果(gg/ml)腐 植 酸 磺酸木质素 试样号实际值预测误差实际值预测误差非相似度因子14.020.130.5070.0620.3622.48 0.100.2660.0040.4234.260.120.3640.0050.544 0.9970.110.203 0.0490.7952.96 0.640.1160.0820.9461.140.00.404 0.0870.6575.080.00.2090.1551.03 84.990.340 0.0751.0890 0.300.7230.119 5.30标准偏差 0.27 0.085其中其中,预测误差为预测浓度与实际浓度之差预测误差为预测浓度
13、与实际浓度之差.如对于小组分磺如对于小组分磺酸木质素酸木质素,平均误差为平均误差为-0.024(-0.024(g/ml)g/ml),相应的标准偏差为相应的标准偏差为0.085(0.085(g/ml)g/ml)。标准偏差所用公式为:标准偏差所用公式为:212)(1测试值真值yyn偏最小二乘方法偏最小二乘方法37而非相似度因子而非相似度因子(dissimilarity factor)dissimilarity factor)的表达式为:的表达式为:)(/22Exssdisa式中式中,s sa a2 2(ExEx)为为X X阵的主成分模型所引进的残余标准方差。阵的主成分模型所引进的残余标准方差。而而
14、s s2 2为为s s2 2=/=/(m m a a)e其中,其中,m m为为X X 的维,的维,a a为主成分数,为主成分数,e e为:为:iiiiptxe 运用运用F F显著性检验,其自由度为显著性检验,其自由度为(m m-a a)/2 )/2 和和(m m-a a)()(n n a a 1)/2 1)/2,显著性水平为显著性水平为,若若s s2 2 s sa a2 2 (ExEx)F F ,则计算值可信。则计算值可信。偏最小二乘方法偏最小二乘方法38 若试样增加一组分,即去污剂(含光白剂),其结果若试样增加一组分,即去污剂(含光白剂),其结果示于表示于表6.26.2。由此表可见,对于腐植
15、酸和磺酸木质素来说,。由此表可见,对于腐植酸和磺酸木质素来说,三组分与二组分浓度预测准确性大体上相当。对于去污剂三组分与二组分浓度预测准确性大体上相当。对于去污剂来说,也得到了较好的结果。在表来说,也得到了较好的结果。在表6.26.2的情况下,由于为三的情况下,由于为三组分混合物,所以构造主成分模型时,也相应增加一因子。组分混合物,所以构造主成分模型时,也相应增加一因子。表表6.2 腐植酸,磺酸木质素和去污剂混合溶液测试结果腐植酸,磺酸木质素和去污剂混合溶液测试结果(ggml)腐植酸磺酸木质素去污剂样本号实际值预测误差 实际值预测误差实际值预测误差非 相 似度因子12.440.080.2890
16、.04280.00.91.2924.080.070.3610.05888.514.43.5131.060.090.2340.02469.21.20.7443.320.070.1230.05940.116.31.2950.9980.080.1460.05930.73.00.9262.980.100.4030.0741202.41.3975.130.320.2290.01549.014.82.2285.060.2000.01751.112.41.85900.490.7350.68899.631.921.3标准偏差 0.22 0.234 14.4偏最小二乘方法偏最小二乘方法39若试样仍如表若试样仍如
17、表6.2,即混合物为三组分,但预测为两个组分,也即混合物为三组分,但预测为两个组分,也就是说构造的预测模型为二因子,其结果示于表就是说构造的预测模型为二因子,其结果示于表6.3。由此表。由此表可见,预测误差反而比表可见,预测误差反而比表6.2为小。原因为:模型中少一因子,为小。原因为:模型中少一因子,所以可使结果更稳定。所以可使结果更稳定。表表6.3 三组分混合物,但仅测试腐植酸和磺酸木质素二组分三组分混合物,但仅测试腐植酸和磺酸木质素二组分(gg/ml)腐 植 酸 磺酸木质素 样本号实际值预测误差实际值预测误差非相似度因子12.440.040.289 0.0161.3124.080.120.
18、3610.1003.5131.060.040.234 0.0440.7543.320.070.1230.0241.125 0.9980.110.416 0.0370.9562.980.100.403 0.0101.5575.130.220.2290.1092.20 85.060 00.0171.80900.350.735 0.532 18.5标准偏差0.15 0.185偏最小二乘方法偏最小二乘方法40 6.4 非线性偏最小二乘非线性偏最小二乘 非线性偏最小二乘与线性偏最小二乘的区别仅仅在于非线性偏最小二乘与线性偏最小二乘的区别仅仅在于X与与Y的内在相关性,即后者为一直线,而前者为一曲线的内在相
19、关性,即后者为一直线,而前者为一曲线,如一抛物线。如一抛物线。曲线的表示有多种数学模型,如二次多项式,三次多项式,曲线的表示有多种数学模型,如二次多项式,三次多项式,指数函数和对数函数等。其中,最简单的为二次多项式:指数函数和对数函数等。其中,最简单的为二次多项式:FQUYEPTXaaaaaaahtctccu2210式中,式中,T,U分别为分别为X,Y的得分矩阵,的得分矩阵,p,Q分别为分别为X,Y的装载的装载矩阵,矩阵,a为某一主成分,这种最简单的二次项扩展的偏最小二为某一主成分,这种最简单的二次项扩展的偏最小二乘可简记为乘可简记为QPLS,QPLS的基本思想是:将的基本思想是:将X和和Y分
20、别投影于分别投影于t 和和u:(1)将将X和和Y分别以分别以tp和和uq近似;近似;(2)同时满足同时满足u和和t 内在的内在的二次相关性。二次相关性。偏最小二乘方法偏最小二乘方法41QPLS的算法为:的算法为:(1)进行数据进行数据X和和Y 的标准化处理;的标准化处理;(2)置因子数(即主成分数)置因子数(即主成分数)a的初值为的初值为0;(3)将将a a增增1:a a=a a+1;.将将Y的某一列(具有最大方差)作为的某一列(具有最大方差)作为u的起始矢量;的起始矢量;.由由PLS计算计算和和t的起始矢量:的起始矢量:uuXu/将将归一化:归一化:=1 =1 /Xt偏最小二乘方法偏最小二乘
21、方法42(4)由最小二乘法测定系数由最小二乘法测定系数c:etctccu2210以以r 标记所计算的标记所计算的u(记为记为 )u 2210tctccr(5)计算计算Y的装载的装载q:rrrYq/将将q 归一化:归一化:1qYY)6(计算新的计算新的u:qqYqu/并由最小二乘法,利用新值计算并由最小二乘法,利用新值计算 和和 。10,cc2c偏最小二乘方法偏最小二乘方法43XY)7(校正校正:eXcXccetctccu22102210)(重写此式,以重写此式,以u=F(X,c)表示,并进行未知参数表示,并进行未知参数c和和的的线性化,在线性化,在F00处可表示为:处可表示为:2 /22211
22、0202010000ikikikikikkiiikkkjjjidxtcsumtcdxcsumtcctctccddFsumdcdFcsumFu 合并合并c00和和c0得新得新c0,合并合并c01ti与与c0ti得新得新c1ti等等等等,由此,由此,上式为:上式为:)2(212210iikikiiidxtccsumtctccu 偏最小二乘方法偏最小二乘方法44此式可用于新此式可用于新的计算:的计算:a.由一维线性由一维线性PLS计算:计算:(i)以以 为前为前k列,以列,以1,ti 和和ti2为后为后3列,建造列,建造Z 矩阵。矩阵。ikixtcc)2(21uuZuv/(ii)将将v归一化:归一化
23、:1vvvZvs/(iii)ssusb/(iv)(v)(前前k 列列)bvd 偏最小二乘方法偏最小二乘方法45b.计算新计算新:=+d 将将归一化:归一化:1 XX)8(/Xt(9)检验收敛:检验收敛:若若 ,或者已经迭代到某一最大,或者已经迭代到某一最大次数(如次数(如50次或次或100次),之后到步次),之后到步(9),否则到步,否则到步(4)。1010/oldoldnewuuu(10)由最后由最后t 值计算新值:值计算新值:r=(步步4),),q(步步6),),u(步步7)和)和c(步步4)。)。u(11)计算装载计算装载t tXtp/偏最小二乘方法偏最小二乘方法46(12)计算残差:计算残差:p tXEqrYF以残差以残差E与与F作为作为X和和Y,若迭代次数若迭代次数 a,回到步回到步(3),否则到,否则到步步(13)。(13)y值预测值预测(i)将新矢量将新矢量x进行预处理(减去模型建立中的平均值,并进行预处理(减去模型建立中的平均值,并被标准偏差除),置预测的被标准偏差除),置预测的y0为为0。(ii)对于每一因子对于每一因子(a)计算:计算:xt 2210tctccruatpXEarqYY