一种基于蒙特卡罗及LASSO的双集成偏最小二乘建模方法与流程

本发明属于分析化学技术领域，具体涉及一种基于蒙特卡罗及LASSO的双集成偏最小二乘建模方法。

背景技术：

光谱分析技术以其简便、快速、绿色以及无损等优点在农业、食品、医药、环境等领域得到了广泛应用。然而，由于光谱吸收峰重叠现象严重，信号吸收较弱，背景干扰严重等问题，因此，在进行复杂样品定性、定量分析时需借助化学计量学手段。而建立一个稳定性好、预测精度高的模型一直是复杂样品定量分析的关键。

传统建模方法采用单一模型在光谱和待测目标值之间建立定量分析模型，预测效果往往不尽如人意。因此，发展了集成建模技术，该技术通过将多个模型的结果进行融合得到最终预测结果。目前主要有基于样品及基于变量方法方向的集成建模方法，前者如boosting PLS（X.G. Shao, X.H. Bian, W.S. Cai, An improved boosting partial least squares method for near-infrared spectroscopic quantitative analysis, Anal. Chim. Acta,2010, 666, 32-37）、bagging ELM（卞希慧，李淑娟，谭小耀，王江江，王治国，刘维国，陈宗蓬，王晨，Bagging极限学习机集成建模方法，2015，中国专利，ZL201510466504），后者如变量聚类-多元线性回归（C. Tan, X. Qin, M. Li, Subspace regression ensemble method based on variable clustering for near-infrared spectroscopic calibration, Analytical Letters,2009, 42,1693-1710）。上述方法只在样品或者变量方向进行单一的集成，如果在这两个方向同时集成建模（G.R.Du, W.S.Cai, X.G.Shao,A variable differential consensus method for improving thequantitative near-infrared spectroscopic analysis,Science China,2012, 55,1946-1952），则有望进一步提高模型的预测精度。

蒙特卡罗（Monte Carlo，MC）方法，也称统计模拟方法、随机抽样技术，是基于“随机数”、以概率统计理论为基础的一种非常重要且被广泛应用的数值计算方法。LASSO变量选择方法是在回归系数的绝对值之和小于一个常数的约束条件下，使残差平方和最小化，从而较严格地使系数本该为0的化为0，相应的变量被删除，进而达到变量选择的目的。通过LASSO方法进行变量选择建立偏最小二乘模型（卞希慧，颜鼎荷，李淑娟，谭小耀，李翔，一种基于LASSO的近红外光谱变量选择方法，2015，中国专利，ZL201510581659），既克服了传统方法在选择变量上的不足，又可以处理高维问题。本发明引进蒙特卡罗及LASSO变量选择，在样品及变量方向进行双集成建模，来进一步提高模型的预测精度。

技术实现要素：

本发明的目的是针对上述存在的问题，提供一种基于蒙特卡罗（MC）及LASSO的样品和变量方向进行双集成的偏最小二乘建模方法（简记为MC-LASSO-PLS），并用于复杂样品定量分析。该方法既保留了两种方法的优势，又克服了单一方法建模预测精确度不高的缺点。

本发明提供的基于蒙特卡罗（MC）及LASSO的样品和变量方向进行双集成的偏最小二乘建模方法，首先采用蒙特卡罗技术选取一定数目的样品作为样品子集，然后采用LASSO技术从该样品子集中选择部分变量，作为样品变量子子集，重复（即迭代）多次，建立多个子模型，将这些模型的预测结果直接作算术平均，即得到最终预测结果。

具体步骤为：

（1）收集一定数目的复杂样品，测量其光谱，确定并测定目标分析物含量，采用一定的分组方式将数据集划分为训练集和预测集；这里，所述复杂样品是指组分多、已知信息少的复杂混合物，如石油、牛奶、中药、农产品、血液等样品；所述分组方式如KS分组、随机分组等；

（2）确定PLS模型的因子数LV：

计算不同因子数下的交叉验证均方根误差（RMSECV），最小的RMSECV对应的因子数即为最佳因子数LV。RMSECV计算公式如下：

其中，m为总样品数，为第i个样品的预测值，为第i个样品的真实值；

（3）确定迭代次数T：

迭代次数（即模型个数）从1变化到500，间隔为1，依次计算不同迭代次数下模型对训练集预测的预测均方根误差（RMSEP），当RMSEP相对稳定时对应的迭代次数为最佳迭代次数T；

（4）确定训练子集样品百分数P：

分别选取训练子集样品总数的5-100%，每隔5%取样一次，依次计算不同训练子集样品百分数下模型对预测集预测的RMSEP，RMSEP最小值对应的训练子集样品百分数为最佳训练子集样品百分数P；

（5）建立模型并预测未知样品

采用确定好的因子数LV、迭代次数T、训练子集样品百分数P，通过MC方法从训练集随机选取最佳训练子集样品百分数P对应的样品数目作为样品子集，然后对样品子集进行LASSO变量选择，得到样品变量子子集，对该子子集建立PLS子模型；

对上述过程重复T次建立T个子模型，并对这T个子模型的预测结果直接平均，得到最终预测结果。

本发明中，所述的复杂样品，可以是石油、烟草、中药、农产品、食品、生物或环境领域的样品等。

本发明方法能有效地改善模型的预测能力，提高模型的预测精度，在预测精度和稳定性方面具有明显优势。本发明适用于石油、烟草、食品、中药等复杂样品定量分析。

附图说明

图1是MC-LASSO-PLS流程图。

图2是玉米近红外数据的训练集光谱图。

图3是玉米近红外数据预测均方根误差值随着迭代次数的变化图。

图4是玉米近红外数据预测均方根误差值随着训练子集样品百分数的变化图。

图5是玉米近红外数据采用不同建模方法预测的预测值与真实值的关系图，其中(a) PLS，(b) LASSO-PLS，(c) MC-LASSO-PLS。

图6是汽油紫外数据的训练集光谱图。

图7是汽油紫外数据预测均方根误差值随着迭代次数的变化图。

图8是汽油紫外数据预测均方根误差值随着训练子集样品百分数变化图。

图9是汽油紫外数据采用不同建模方法预测的预测值与真实值的关系图，其中(a) PLS，(b) LASSO-PLS，(c) MC-LASSO-PLS。

图10是四元调和油近红外的训练集光谱图。

图11是四元调和油近红外数据预测均方根误差值随着迭代次数的变化图。

图12是四元调和油近红外数据预测均方根误差值随着训练集样品百分数的变化图。

图13是四元调和油近红外光谱数据采用不同建模方法预测的预测值与真实值的关系图，其中(a) PLS，(b) LASSO-PLS，(c) MC-LASSO-PLS。

图14是烟草近红外数据的训练集光谱图。

图15是烟草近红外数据预测均方根误差值随着迭代次数的变化图。

图16是烟草近红外数据预测均方根误差值随着训练集样品百分数的变化图。

图17是烟草近红外数据采用不同建模方法预测的预测值与真实值的关系图，其中(a) PLS，(b) LASSO-PLS，(c) MC-LASSO-PLS。

具体实施方式

为更好地理解本发明，下面结合实施例对本发明做进一步地详细说明，但是本发明要求保护的范围并不局限于实施例所表示的范围。

实施例1：

本实施例是应用于近红外光谱数据分析，对玉米样品中油含量进行测定。具体的步骤如下：

1）收集80个玉米样品，采用三个不同的近红外光谱仪（M5，MP5，MP6）测量玉米的近红外光谱数据，将油含量作为目标值。近红外光谱的波长范围为2498～1100nm（4003～9091cm^-1），采样间隔为2nm，共700个波长数据点。该数据从http: // software .eigenvector. com /Data/ Corn / index. html下载。采用KS分组方法，53个样品用作训练集，剩余27个品用作预测集。该数据的训练集近红外光谱如图2所示。

2）确定PLS模型的因子数LV

计算不同因子数下的交叉验证均方根误差（RMSECV），最小的RMSECV对应的因子数即为最佳因子数。其中RMSECV计算公式如下：

其中m为总样品数，为第i个样品的预测值，为第i个样品的真实值。该实施例中最佳因子数LV为6。

3）确定迭代次数T

迭代次数（即模型个数）从1变化到500，间隔为1，依次计算不同迭代次数下模型对训练集预测的预测均方根误差（RMSEP）。图3显示了RMSEP值随着迭代次数的变化，当RMSEP相对稳定时对应的迭代次数500即为最佳迭代次数T。

4）确定训练子集样品百分数P

分别选取训练子集样品总数的5-100%，每隔5%取样一次，依次计算不同训练子集样品百分数下模型对训练集预测的RMSEP。图4显示了RMSEP值随着迭代次数的变化，RMSEP最小值对应的训练子集样品百分数55%即为最佳训练子集样品百分数P。

5）建立模型并预测未知样品

采用确定好的因子数6、迭代次数500、训练子集样品百分数55%，通过MC方法从训练集随机选取最佳训练子集样品百分数55%对应的30个样品作为样品子集，然后对样品子集进行LASSO变量选择，得到样品变量子子集，对该子子集建立PLS子模型。对上述过程重复500次建立500个子模型。对这500个子模型的预测结果直接平均得到最终预测结果。

图5 (a)、(b)及(c)分别显示了PLS、LASSO-PLS及MC-LASSO-PLS三种不同建模方法对预测集进行预测的预测值与真实值的关系图。结果表明，MC-LASSO明显的提高了PLS模型的预测能力，达到较好的预测效果。

实施例2：

本实施例是应用于紫外光谱数据分析，对汽油样品中的单环芳烃含量进行测定。具体的步骤如下：

1）收集115个轻汽油和柴油燃料样品，紫外光谱波长范围为200～400nm，采样间隔约为0.35nm，共572个波长数据点。单环芳烃含量通过HP model G1205A超临界流体色谱仪器（Hewlett-Packard，Palo Alto，Calif）测定。该数据从http://myweb.dal.ca/pdwentze/downloads.html下载。训练集及预测集的划分按照网上的说明，前70个样品用作训练集，后44个样品用作预测集。该数据的训练集紫外光谱如图6所示。

2）确定PLS模型的因子数LV

计算不同因子数下的交叉验证均方根误差（RMSECV），最小的RMSECV对应的因子数即为最佳因子数。其中RMSECV计算公式如下：

其中m为总样品数，为第i个样品的预测值，为第i个样品的真实值。该实施例中最佳因子数LV为7。

3）确定迭代次数T

迭代次数（即模型个数）从1变化到500，间隔为1，依次计算不同迭代次数下模型对训练集预测的预测均方根误差（RMSEP），图7显示了RMSEP值随着迭代次数的变化，当RMSEP相对稳定时对应的迭代次数500即为最佳迭代次数T。

4）确定训练子集样品百分数P

分别选取训练子集样品总数的5-100%，每隔5%取样一次，依次计算不同训练子集样品百分数下模型对训练集预测的RMSEP，图8显示了RMSEP值随着训练子集样品百分数变化，RMSEP最小值对应的训练子集样品百分数40%即为最佳训练子集样品百分数P。

5）建立模型并预测未知样品

采用确定好的因子数7、迭代次数500、训练子集样品百分数40%，通过MC方法从训练集随机选取最佳训练子集样品百分数40%对应的28个样品作为样品子集，然后对样品子集进行LASSO变量选择，得到样品变量子子集，对该子子集建立PLS子模型。对上述过程重复500次来建立500个子模型。对这500个子模型的预测结果直接平均得到最终预测结果。

图9 (a)、(b)及(c)分别显示了PLS、LASSO-PLS及MC-LASSO-PLS三种不同建模方法对预测集进行预测的预测值与真实值的关系图。结果表明，MC-LASSO明显提高了PLS模型的预测能力，达到较好的预测效果。

实施例3：

本实施例是应用于近红外光谱数据分析，对四元调和油样品中的香油品含量进行测定。具体的步骤如下：

1）收集51个含有香油、玉米油、大豆油和稻米油的四元调和油样品。使用Vertex70多波段红外/近红外光谱仪（德国Bruker公司）进行近红外光谱数据测量，波数范围为4000～12000cm^-1，采样间隔为1.93cm^-1，共4148个数据点。将香油含量作为目标值。采用KS分组方法，34个样品用作训练集，剩余17个样品用作预测集。该数据的训练集近红外光谱如图10所示。

2）确定PLS模型的因子数LV

计算不同因子数下的交叉验证均方根误差（RMSECV），最小的RMSECV对应的因子数即为最佳因子数。其中RMSECV计算公式如下：

其中m为总样品数，为第i个样品的预测值，为第i个样品的真实值。该实施例中最佳因子数LV为8。

3）确定迭代次数T

迭代次数（即模型个数）从1变化到500，间隔为1，依次计算不同迭代次数下模型对训练集预测的预测均方根误差（RMSEP），图11显示了RMSEP值随着迭代次数的变化，当RMSEP相对稳定时对应的迭代次数500即为最佳迭代次数T。

4）确定训练子集样品百分数P

分别选取训练子集样品总数的5-100%，每隔5%取样一次，依次计算不同训练子集样品百分数下模型对训练集预测的RMSEP，图12显示了RMSEP值随着训练集样品百分数的变化，RMSEP最小值对应的训练子集样品百分数65%即为最佳训练子集样品百分数P。

5）建立模型并预测未知样品

采用确定好的因子数8、迭代次数500、训练子集样品百分数65%，通过MC方法从训练集随机选取最佳训练子集样品百分数65%对应的23个样品作为样品子集，然后对样品子集进行LASSO变量选择，得到样品变量子子集，对该子子集建立PLS子模型。对上述过程重复500次建立500个子模型。对这500个子模型的预测结果直接平均得到最终预测结果。

图13 (a)、(b)及(c)分别显示了PLS、LASSO-PLS及MC-LASSO-PLS三种不同建模方法对预测集进行预测的预测值与真实值的关系图。结果表明，MC-LASSO明显提高了PLS模型的预测能力，达到较好的预测效果。

实施例4：

本实施例是应用于近红外光谱数据分析，对烟草样品中的尼古丁含量进行测定。具体的步骤如下：

1）收集373个烟草叶片样品。采用Bruker Vector 22/N的傅里叶近红外光谱仪（德国Bruker公司）进行测定波数范围为4000～9000cm^-1，采样间隔为3.9cm^-1，共1296个数据点。烟草样品中尼古丁含量采用AAⅢ型连续流动分析仪（德国Bran Luebbe公司）按照标准方法测定。采用KS分组方法，249个样品用作训练集，剩余124个样品用作预测集。该数据的训练集近红外光谱如图14所示。

2）确定PLS模型的因子数LV

计算不同因子数下的交叉验证均方根误差（RMSECV），最小的RMSECV对应的因子数即为最佳因子数。其中RMSECV计算公式如下：

其中m为总样品数，为第i个样品的预测值，为第i个样品的真实值。该实施例中最佳因子数LV为8。

3）确定迭代次数T

迭代次数（即模型个数）从1变化到500，间隔为1，依次计算不同迭代次数下模型对训练集预测的预测均方根误差（RMSEP），图15显示了RMSEP值随着迭代次数的变化，当RMSEP相对稳定时对应的迭代次数500即为最佳迭代次数T。

4）确定训练子集样品百分数P

分别选取训练子集样品总数的5-100%，每隔5%取样一次，依次计算不同训练子集样品百分数下模型对训练集预测的RMSEP，图16显示了RMSEP值随着训练集样品百分数的变化，RMSEP最小值对应的训练子集样品百分数45%即为最佳训练子集样品百分数P。

5）建立模型并预测未知样品

采用确定好的因子数8、迭代次数500、训练子集样品百分数45%，通过MC方法从训练集随机选取最佳训练子集样品百分数45%对应的23个样品作为样品子集，然后对样品子集进行LASSO变量选择，得到样品变量子子集，对该子子集建立PLS子模型。对上述过程重复500次建立500个子模型。对这500个子模型的预测结果直接平均得到最终预测结果。

图17中(a)、(b)及(c)分别显示了PLS、LASSO-PLS及MC-LASSO-PLS三种不同建模方法对预测集进行预测的预测值与真实值的关系图。结果表明，MC-LASSO提高了PLS模型的预测能力，达到较好的预测效果。