黃歡，黃宇

（華測(cè)檢測(cè)認(rèn)證集團(tuán)股份有限公司，標(biāo)準(zhǔn)物質(zhì)研究中心，深圳 518000）

測(cè)量不確定度的評(píng)定常因其較高的數(shù)學(xué)和統(tǒng)計(jì)學(xué)要求，讓實(shí)驗(yàn)室工作人員望而卻步。然而，GB/T 27025—2019 和RB/T 214—2017 等標(biāo)準(zhǔn)規(guī)定了實(shí)驗(yàn)室需要評(píng)估和報(bào)告測(cè)量不確定度，以證明其操作能力符合要求以及測(cè)試結(jié)果有效［1-2］。

目前，測(cè)量不確定度的評(píng)定使用較多的仍是JJF 1059.1—2012 《測(cè)量不確定度評(píng)定與表示》中規(guī)定的方法（《測(cè)量不確定度表示指南》GUM 法）［3］。GUM 法主要根據(jù)不確定度傳播定律要求，通過(guò)建立計(jì)算模型，分析不確定度來(lái)源，確定傳播系數(shù)，從而計(jì)算輸入量貢獻(xiàn)的不確定度分量，合成輸出量標(biāo)準(zhǔn)不確定度及擴(kuò)展不確定度等過(guò)程來(lái)評(píng)定測(cè)量不確定度［4］。

近年來(lái)，隨著計(jì)算機(jī)和計(jì)算軟件的迭代升級(jí)，蒙特卡洛法（MCM）評(píng)定不確定度在物理、化學(xué)、材料等領(lǐng)域得到了廣泛的應(yīng)用，其評(píng)定結(jié)果與GUM法基本一致［5-8］。MCM 作為GUM 法的重要補(bǔ)充，無(wú)需計(jì)算傳播系數(shù)或近似處理計(jì)算模型，可得到包含實(shí)際采樣信息在內(nèi)的，且更為全面準(zhǔn)確的計(jì)算結(jié)果。根據(jù)JJF1059.2—2012 要求，MCM 采用“概率分布傳播”的通用方法來(lái)評(píng)估被測(cè)量的測(cè)量不確定度。它通過(guò)對(duì)輸入量的概率密度函數(shù)（PDF）離散抽樣，由測(cè)量模型傳播輸入量的分布，計(jì)算獲得輸出量PDF 的離散抽樣值，進(jìn)而由輸出量的離散分布數(shù)值直接獲取輸出量的最佳估計(jì)值、標(biāo)準(zhǔn)不確定度和包含區(qū)間［9］。

目前，MCM 評(píng)定不確定度主要是在MATLAB（矩陣實(shí)驗(yàn)室）軟件平臺(tái)上實(shí)現(xiàn)的［10-12］，但由于其付費(fèi)使用門檻而未能被廣泛使用。R 語(yǔ)言作為一種開(kāi)源的統(tǒng)計(jì)編程語(yǔ)言，經(jīng)過(guò)20 多年的發(fā)展，已成為數(shù)據(jù)科學(xué)領(lǐng)域最受歡迎的工具之一［13-14］。相比于MATLAB，R 語(yǔ)言具有開(kāi)源、硬件條件要求低、跨平臺(tái)、程序庫(kù)包豐富、可視化功能完備等優(yōu)勢(shì)，完全適用于實(shí)驗(yàn)人員對(duì)于不確定度評(píng)定的數(shù)據(jù)計(jì)算［15］。截至目前，尚無(wú)采用基于R 語(yǔ)言的MCM 評(píng)定不確定度的研究報(bào)道。筆者采用鹽酸溶液濃度標(biāo)定實(shí)驗(yàn)為例［4］，以GUM 法評(píng)定作為參照，研究了基于R 語(yǔ)言的蒙特卡洛法（MCM）在測(cè)量不確定度評(píng)定中的實(shí)際應(yīng)用。

1 鹽酸溶液濃度的標(biāo)定實(shí)驗(yàn)

1.1 測(cè)試流程

HCl 溶液標(biāo)定的流程如圖1 所示。離子水中，以NaOH 溶液進(jìn)行滴定。滴定裝置自動(dòng)控制NaOH 的加入量，同時(shí)繪出pH 曲線。通過(guò)記錄的pH 曲線形狀確定終點(diǎn)。滴定消耗的氫氧化鈉溶液體積為18.64 mL。

圖1 HCl 溶液標(biāo)定流程

（4）用移液管移取15 mL HCl 溶液，用去離子水稀釋至約50 mL 于滴定瓶中，以標(biāo)定后的NaOH溶液進(jìn)行滴定。用同一臺(tái)自動(dòng)滴定裝置，繪制pH標(biāo)準(zhǔn)工作曲線并判定滴定終點(diǎn)。滴定消耗的氫氧化鈉體積為14.89 mL。

1.2 鹽酸溶液標(biāo)定涉及的相關(guān)參數(shù)

天平線性分量：±15 mg；KHP 純度（P）：(100±0.05)%；滴定管最大允許誤差：±0.03 mL；相對(duì)分子質(zhì)量：M(C) = 12.010 7(8)，M(H) = 1.0079 4(7)，M(O) = 15.999 4(3)，M(K) = 39.098 3(1)；環(huán)境溫度變化范圍：±4 ℃；水的膨脹系數(shù)：（2.1×10-4） ℃-1。

1.3 HCl 溶液標(biāo)定的計(jì)算模型

（1）準(zhǔn)備0.1 mol/L 數(shù)量級(jí)的待標(biāo)定NaOH 溶液和HCl 溶液。

（2）干燥滴定所需的一級(jí)標(biāo)準(zhǔn)物鄰苯二甲酸氫鉀（KHP），以確保其純度符合供應(yīng)商提供的證書(shū)上所標(biāo)數(shù)值。稱取大約0.388 g 干燥的標(biāo)準(zhǔn)物KHP 用于標(biāo)定19 mL NaOH 溶液。

（3）將滴定標(biāo)準(zhǔn)物KHP 溶解于約50 mL 的去

2 結(jié)果與討論

2.1 GUM 法測(cè)量不確定度評(píng)定

GUM 法輸入量及其不確定度來(lái)源分析如表1所示。表1 中列出了各輸入量數(shù)值、不確定度來(lái)源、標(biāo)準(zhǔn)不確定度u(x)和相對(duì)標(biāo)準(zhǔn)不確定度u(x)/x，將表1 中各數(shù)值代入到公式(1)，計(jì)算得：

表1 GUM 法輸入量及其不確定度來(lái)源分析

進(jìn)一步計(jì)算合成相對(duì)標(biāo)準(zhǔn)不確定度：

2.2 MCM 測(cè)量不確定度評(píng)定

R 語(yǔ)言是一種用于統(tǒng)計(jì)計(jì)算和畫(huà)圖的開(kāi)源語(yǔ)言和環(huán)境，它可用于Windows（視窗操作系統(tǒng)）、MacOS（蘋(píng)果操作系統(tǒng)）和各種UNIX（多任務(wù)多用戶操作系統(tǒng)）平臺(tái)。反刪除和數(shù)據(jù)恢復(fù)軟件系列RStudio為R 語(yǔ)言開(kāi)發(fā)提供了一個(gè)易操作的、集成的開(kāi)發(fā)環(huán)境，可實(shí)現(xiàn)代碼編寫(xiě)、變量值監(jiān)控、控制臺(tái)和圖像輸出。試驗(yàn)中提供的所有代碼均是在RStudio 軟件中開(kāi)發(fā)的。

在R 語(yǔ)言中，可以非常容易地通過(guò)特定函數(shù)程序包來(lái)實(shí)現(xiàn)特定PDF 的隨機(jī)數(shù)樣本抽取。如使用runif 函數(shù)產(chǎn)生均勻分布隨機(jī)數(shù)、rnorm 函數(shù)產(chǎn)生正態(tài)分布隨機(jī)數(shù)、rtriangle 函數(shù)產(chǎn)生三角分布隨機(jī)數(shù)、rexp 函數(shù)產(chǎn)生指數(shù)分布隨機(jī)數(shù)、RT 函數(shù)產(chǎn)生t分布隨機(jī)數(shù)等，這些均為實(shí)現(xiàn)蒙特卡洛模擬提供了便利。MCM 各輸入量數(shù)值及其不確定度來(lái)源分析如表2所示，分別給出了不同輸入量的數(shù)值、不確定度來(lái)源以及輸入量的PDF 離散抽樣。

表2 MCM 輸入量及其不確定度來(lái)源分析

使用MCM 獲得輸出量估計(jì)值及其不確定度，關(guān)鍵是要確定輸入量的參數(shù)信息，包括概率分布類型和相關(guān)參數(shù)。如重復(fù)性服從均值為1，標(biāo)準(zhǔn)差為0.001 的正態(tài)分布；KHP 質(zhì)量服從中心為0.388 8 g，半寬為0.15 mg 的均勻分布；滴定HCl 所用NaOH溶液的體積存在兩種概率分布：三角分布（中心為14.89 mL，半寬為0.03 mL）和均勻分布（中心為14.89 mL，半寬為0.013 mL）。R 語(yǔ)言實(shí)現(xiàn)蒙特卡洛評(píng)定是模擬對(duì)輸入量的PDF 離散抽樣，獲得輸出量的PDF 的離散抽樣值，進(jìn)而可以計(jì)算輸出量的最佳估計(jì)值、標(biāo)準(zhǔn)不確定度和包含區(qū)間。

R 代碼實(shí)現(xiàn)過(guò)程如下：

這個(gè)過(guò)程主要通過(guò)查找排序后的輸出量中是否存在寬度不大于95%對(duì)稱包含區(qū)間的任一95%包含區(qū)間來(lái)實(shí)現(xiàn)的。結(jié)果表明，輸出量95%對(duì)稱包含區(qū)間就是最短包含區(qū)間，輸出量的概率分布是對(duì)稱的。圖2 為輸出量概率密度函數(shù)。

圖2 輸出量概率函數(shù)與理論正態(tài)分布函數(shù)

從圖2 中可以看出，輸出量概率密度函數(shù)與理論正態(tài)分布函數(shù)幾近重合。將輸出量95%包含區(qū)間的半寬作為擴(kuò)展不確定度U與輸出量標(biāo)準(zhǔn)差（標(biāo)準(zhǔn)不確定度u）進(jìn)行比較，包含因子k=U/u=1.95，與正態(tài)分布函數(shù)包含因子k（95%）=1.96 非常接近，代碼如#4：

結(jié)果表明，輸出量的不確定度受到服從正態(tài)分布的輸入量（重復(fù)性）的影響最大。采用GUM 法評(píng)估的不確定度分量也是重復(fù)性不確定度分量最大，這與MCM 的結(jié)果是一致的。

上述不確定度的蒙特卡洛評(píng)定在許多編程語(yǔ)言上均能夠?qū)崿F(xiàn)，如C、C++、python 以及MATLAB等［16-17］。但試驗(yàn)中采用R 語(yǔ)言編程來(lái)實(shí)現(xiàn)蒙特卡洛評(píng)定，不僅因其功能的豐富性和開(kāi)源性，還在于其在向量和矩陣運(yùn)算中的高效性和簡(jiǎn)潔性，這在一定程度上縮短了大量復(fù)雜模型和模擬數(shù)據(jù)的計(jì)算時(shí)間，同時(shí)還規(guī)避了一些邏輯語(yǔ)句的使用，對(duì)于初學(xué)者十分友好。因此，試驗(yàn)中也使用向量、矩陣和數(shù)據(jù)框運(yùn)算方式的R 代碼來(lái)實(shí)現(xiàn)蒙特卡洛評(píng)定，代碼如下：

在代碼#5 中，不再使用for 循環(huán)語(yǔ)句，而是直接離散抽樣或離散抽樣后通過(guò)矩陣或數(shù)據(jù)框行加和得到輸入量離散抽樣數(shù)值，將這些數(shù)值向量代入計(jì)算模型計(jì)算，進(jìn)而獲得輸出量PDF 的離散抽樣值，最后計(jì)算輸出量的最佳估計(jì)值、標(biāo)準(zhǔn)不確定度和包含區(qū)間。從計(jì)算結(jié)果來(lái)看，采用向量、矩陣和數(shù)據(jù)框運(yùn)算方式編寫(xiě)的代碼#5 與常規(guī)編程方式編寫(xiě)的代碼#1 相同。#5 代碼編寫(xiě)方式充分利用了R 語(yǔ)言的優(yōu)勢(shì)，在一定程度上降低了編碼難度的同時(shí)，提升了計(jì)算的效率。

3 結(jié)論

以鹽酸溶液濃度標(biāo)定試驗(yàn)為例，研究了基于R語(yǔ)言的蒙特卡洛法（MCM）在測(cè)量不確定度評(píng)定中的實(shí)際應(yīng)用。MCM 的計(jì)算過(guò)程簡(jiǎn)潔易行，無(wú)需考慮各分量間的相關(guān)性以及簡(jiǎn)化影響因素，善于非線性模型問(wèn)題的計(jì)算；而R 語(yǔ)言作為一種開(kāi)源的統(tǒng)計(jì)編程語(yǔ)言，計(jì)算高效且編程友好，兩者的實(shí)踐結(jié)合在本研究的范例中得到了有效的驗(yàn)證。結(jié)果表明，基于R 語(yǔ)言的MCM 與傳統(tǒng)的GUM 法對(duì)于鹽酸溶液濃度標(biāo)定實(shí)驗(yàn)的測(cè)量不確定度的評(píng)定具有一致性。同時(shí)，MCM 算得的輸出量PDF 與理論正態(tài)分布函數(shù)基本重合，且輸出量PDF 半寬與標(biāo)準(zhǔn)不確定度的比值k=1.95，表明服從正態(tài)分布的重復(fù)性不確定度分量對(duì)計(jì)算結(jié)果影響最大。

采用R 語(yǔ)言擅長(zhǎng)的向量和矩陣運(yùn)算方式編程可以提升計(jì)算效率、降低編程難度。隨著R 語(yǔ)言的的快速發(fā)展和廣泛使用，采用基于R 語(yǔ)言的MCM來(lái)評(píng)定不確定度的研究必將會(huì)越來(lái)越多，以期本研究使用的范例以及R 代碼能為類似研究提供一些參考。