基于R 语言的蒙特卡洛法在不确定度评定中的应用

黄欢，黄宇

（华测检测认证集团股份有限公司，标准物质研究中心，深圳 518000）

测量不确定度的评定常因其较高的数学和统计学要求，让实验室工作人员望而却步。然而，GB/T 27025—2019 和RB/T 214—2017 等标准规定了实验室需要评估和报告测量不确定度，以证明其操作能力符合要求以及测试结果有效［1-2］。

目前，测量不确定度的评定使用较多的仍是JJF 1059.1—2012 《测量不确定度评定与表示》中规定的方法（《测量不确定度表示指南》GUM 法）［3］。GUM 法主要根据不确定度传播定律要求，通过建立计算模型，分析不确定度来源，确定传播系数，从而计算输入量贡献的不确定度分量，合成输出量标准不确定度及扩展不确定度等过程来评定测量不确定度［4］。

近年来，随着计算机和计算软件的迭代升级，蒙特卡洛法（MCM）评定不确定度在物理、化学、材料等领域得到了广泛的应用，其评定结果与GUM法基本一致［5-8］。MCM 作为GUM 法的重要补充，无需计算传播系数或近似处理计算模型，可得到包含实际采样信息在内的，且更为全面准确的计算结果。根据JJF1059.2—2012 要求，MCM 采用“概率分布传播”的通用方法来评估被测量的测量不确定度。它通过对输入量的概率密度函数（PDF）离散抽样，由测量模型传播输入量的分布，计算获得输出量PDF 的离散抽样值，进而由输出量的离散分布数值直接获取输出量的最佳估计值、标准不确定度和包含区间［9］。

目前，MCM 评定不确定度主要是在MATLAB（矩阵实验室）软件平台上实现的［10-12］，但由于其付费使用门槛而未能被广泛使用。R 语言作为一种开源的统计编程语言，经过20 多年的发展，已成为数据科学领域最受欢迎的工具之一［13-14］。相比于MATLAB，R 语言具有开源、硬件条件要求低、跨平台、程序库包丰富、可视化功能完备等优势，完全适用于实验人员对于不确定度评定的数据计算［15］。截至目前，尚无采用基于R 语言的MCM 评定不确定度的研究报道。笔者采用盐酸溶液浓度标定实验为例［4］，以GUM 法评定作为参照，研究了基于R 语言的蒙特卡洛法（MCM）在测量不确定度评定中的实际应用。

1 盐酸溶液浓度的标定实验

1.1 测试流程

HCl 溶液标定的流程如图1 所示。离子水中，以NaOH 溶液进行滴定。滴定装置自动控制NaOH 的加入量，同时绘出pH 曲线。通过记录的pH 曲线形状确定终点。滴定消耗的氢氧化钠溶液体积为18.64 mL。

图1 HCl 溶液标定流程

（4）用移液管移取15 mL HCl 溶液，用去离子水稀释至约50 mL 于滴定瓶中，以标定后的NaOH溶液进行滴定。用同一台自动滴定装置，绘制pH标准工作曲线并判定滴定终点。滴定消耗的氢氧化钠体积为14.89 mL。

1.2 盐酸溶液标定涉及的相关参数

天平线性分量：±15 mg；KHP 纯度（P）：(100±0.05)%；滴定管最大允许误差：±0.03 mL；相对分子质量：M(C) = 12.010 7(8)，M(H) = 1.0079 4(7)，M(O) = 15.999 4(3)，M(K) = 39.098 3(1)；环境温度变化范围：±4 ℃；水的膨胀系数：（2.1×10-4） ℃-1。

1.3 HCl 溶液标定的计算模型

（1）准备0.1 mol/L 数量级的待标定NaOH 溶液和HCl 溶液。

（2）干燥滴定所需的一级标准物邻苯二甲酸氢钾（KHP），以确保其纯度符合供应商提供的证书上所标数值。称取大约0.388 g 干燥的标准物KHP 用于标定19 mL NaOH 溶液。

（3）将滴定标准物KHP 溶解于约50 mL 的去

2 结果与讨论

2.1 GUM 法测量不确定度评定

GUM 法输入量及其不确定度来源分析如表1所示。表1 中列出了各输入量数值、不确定度来源、标准不确定度u(x)和相对标准不确定度u(x)/x，将表1 中各数值代入到公式(1)，计算得：

表1 GUM 法输入量及其不确定度来源分析

进一步计算合成相对标准不确定度：

2.2 MCM 测量不确定度评定

R 语言是一种用于统计计算和画图的开源语言和环境，它可用于Windows（视窗操作系统）、MacOS（苹果操作系统）和各种UNIX（多任务多用户操作系统）平台。反删除和数据恢复软件系列RStudio为R 语言开发提供了一个易操作的、集成的开发环境，可实现代码编写、变量值监控、控制台和图像输出。试验中提供的所有代码均是在RStudio 软件中开发的。

在R 语言中，可以非常容易地通过特定函数程序包来实现特定PDF 的随机数样本抽取。如使用runif 函数产生均匀分布随机数、rnorm 函数产生正态分布随机数、rtriangle 函数产生三角分布随机数、rexp 函数产生指数分布随机数、RT 函数产生t分布随机数等，这些均为实现蒙特卡洛模拟提供了便利。MCM 各输入量数值及其不确定度来源分析如表2所示，分别给出了不同输入量的数值、不确定度来源以及输入量的PDF 离散抽样。

表2 MCM 输入量及其不确定度来源分析

使用MCM 获得输出量估计值及其不确定度，关键是要确定输入量的参数信息，包括概率分布类型和相关参数。如重复性服从均值为1，标准差为0.001 的正态分布；KHP 质量服从中心为0.388 8 g，半宽为0.15 mg 的均匀分布；滴定HCl 所用NaOH溶液的体积存在两种概率分布：三角分布（中心为14.89 mL，半宽为0.03 mL）和均匀分布（中心为14.89 mL，半宽为0.013 mL）。R 语言实现蒙特卡洛评定是模拟对输入量的PDF 离散抽样，获得输出量的PDF 的离散抽样值，进而可以计算输出量的最佳估计值、标准不确定度和包含区间。

R 代码实现过程如下：

这个过程主要通过查找排序后的输出量中是否存在宽度不大于95%对称包含区间的任一95%包含区间来实现的。结果表明，输出量95%对称包含区间就是最短包含区间，输出量的概率分布是对称的。图2 为输出量概率密度函数。

图2 输出量概率函数与理论正态分布函数

从图2 中可以看出，输出量概率密度函数与理论正态分布函数几近重合。将输出量95%包含区间的半宽作为扩展不确定度U与输出量标准差（标准不确定度u）进行比较，包含因子k=U/u=1.95，与正态分布函数包含因子k（95%）=1.96 非常接近，代码如#4：

结果表明，输出量的不确定度受到服从正态分布的输入量（重复性）的影响最大。采用GUM 法评估的不确定度分量也是重复性不确定度分量最大，这与MCM 的结果是一致的。

上述不确定度的蒙特卡洛评定在许多编程语言上均能够实现，如C、C++、python 以及MATLAB等［16-17］。但试验中采用R 语言编程来实现蒙特卡洛评定，不仅因其功能的丰富性和开源性，还在于其在向量和矩阵运算中的高效性和简洁性，这在一定程度上缩短了大量复杂模型和模拟数据的计算时间，同时还规避了一些逻辑语句的使用，对于初学者十分友好。因此，试验中也使用向量、矩阵和数据框运算方式的R 代码来实现蒙特卡洛评定，代码如下：

在代码#5 中，不再使用for 循环语句，而是直接离散抽样或离散抽样后通过矩阵或数据框行加和得到输入量离散抽样数值，将这些数值向量代入计算模型计算，进而获得输出量PDF 的离散抽样值，最后计算输出量的最佳估计值、标准不确定度和包含区间。从计算结果来看，采用向量、矩阵和数据框运算方式编写的代码#5 与常规编程方式编写的代码#1 相同。#5 代码编写方式充分利用了R 语言的优势，在一定程度上降低了编码难度的同时，提升了计算的效率。

3 结论

以盐酸溶液浓度标定试验为例，研究了基于R语言的蒙特卡洛法（MCM）在测量不确定度评定中的实际应用。MCM 的计算过程简洁易行，无需考虑各分量间的相关性以及简化影响因素，善于非线性模型问题的计算；而R 语言作为一种开源的统计编程语言，计算高效且编程友好，两者的实践结合在本研究的范例中得到了有效的验证。结果表明，基于R 语言的MCM 与传统的GUM 法对于盐酸溶液浓度标定实验的测量不确定度的评定具有一致性。同时，MCM 算得的输出量PDF 与理论正态分布函数基本重合，且输出量PDF 半宽与标准不确定度的比值k=1.95，表明服从正态分布的重复性不确定度分量对计算结果影响最大。

采用R 语言擅长的向量和矩阵运算方式编程可以提升计算效率、降低编程难度。随着R 语言的的快速发展和广泛使用，采用基于R 语言的MCM来评定不确定度的研究必将会越来越多，以期本研究使用的范例以及R 代码能为类似研究提供一些参考。