基于R 語言的蒙特卡洛法在不確定度評定中的應用

黃歡，黃宇

（華測檢測認證集團股份有限公司，標準物質(zhì)研究中心，深圳 518000）

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

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

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

目前，MCM 評定不確定度主要是在MATLAB（矩陣實驗室）軟件平臺上實現(xiàn)的［10-12］，但由于其付費使用門檻而未能被廣泛使用。R 語言作為一種開源的統(tǒng)計編程語言，經(jīng)過20 多年的發(fā)展，已成為數(shù)據(jù)科學領(lǐng)域最受歡迎的工具之一［13-14］。相比于MATLAB，R 語言具有開源、硬件條件要求低、跨平臺、程序庫包豐富、可視化功能完備等優(yōu)勢，完全適用于實驗人員對于不確定度評定的數(shù)據(jù)計算［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 鹽酸溶液標定涉及的相關(guān)參數(shù)

天平線性分量：±15 mg；KHP 純度（P）：(100±0.05)%；滴定管最大允許誤差：±0.03 mL；相對分子質(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 溶液標定的計算模型

（1）準備0.1 mol/L 數(shù)量級的待標定NaOH 溶液和HCl 溶液。

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

（3）將滴定標準物KHP 溶解于約50 mL 的去

2 結(jié)果與討論

2.1 GUM 法測量不確定度評定

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

表1 GUM 法輸入量及其不確定度來源分析

進一步計算合成相對標準不確定度：

2.2 MCM 測量不確定度評定

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

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

表2 MCM 輸入量及其不確定度來源分析

使用MCM 獲得輸出量估計值及其不確定度，關(guān)鍵是要確定輸入量的參數(shù)信息，包括概率分布類型和相關(guān)參數(shù)。如重復性服從均值為1，標準差為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 語言實現(xiàn)蒙特卡洛評定是模擬對輸入量的PDF 離散抽樣，獲得輸出量的PDF 的離散抽樣值，進而可以計算輸出量的最佳估計值、標準不確定度和包含區(qū)間。

R 代碼實現(xiàn)過程如下：

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

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

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

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

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

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

3 結(jié)論

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

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