基于Python軟件的蒙特卡洛法不確定度評(píng)定

王舵

(遼寧省計(jì)量科學(xué)研究院，沈陽(yáng) 110004)

蒙特卡洛法(MCM)作為《測(cè)量不確定度表示指南》（GUM）的重要補(bǔ)充方法，可有效地解決測(cè)量模型非線性等諸多不確定度評(píng)定問(wèn)題［1］。近年來(lái)，國(guó)內(nèi)外對(duì)MCM評(píng)定不確定度做了許多研究，在物理學(xué)、生物學(xué)、工程科學(xué)等領(lǐng)域MCM得到了廣泛的應(yīng)用［2］。研究結(jié)果表明，與GUM相比，MCM沒(méi)有明顯的不同，具有很好的評(píng)定效果［3–6］。但MCM計(jì)算簡(jiǎn)單，不用考慮各分量之間的相關(guān)性，不需為簡(jiǎn)化評(píng)定過(guò)程而忽略某些影響因素，評(píng)定過(guò)程更加方便易行［7］。目前，相關(guān)的研究主要實(shí)施在MATLAB計(jì)算平臺(tái)上［8–11］。MATLAB作為一款商業(yè)軟件在各個(gè)方面都十分完善，但其使用費(fèi)用并不能被人們廣泛接受。隨著人工智能的飛速發(fā)展，Python軟件作為一款開(kāi)源軟件在科學(xué)計(jì)算領(lǐng)域占據(jù)了主導(dǎo)地位。與商業(yè)軟件MATLAB相比，Python具有免費(fèi)、跨平臺(tái)、函數(shù)庫(kù)完善等諸多優(yōu)點(diǎn)，非常適用于科學(xué)計(jì)算［12］。毛細(xì)管電泳是一種效率極高的分離分析技術(shù)，在很多領(lǐng)域都有廣泛的應(yīng)用［13］。多數(shù)的研究都集中在分析毛細(xì)管電泳分離過(guò)程中各因素對(duì)實(shí)驗(yàn)結(jié)果不確定度的影響［14］，但對(duì)毛細(xì)管電泳法檢出限的不確定度評(píng)定卻鮮有報(bào)道。與直接以標(biāo)準(zhǔn)物質(zhì)為研究對(duì)象的檢出限不確定度評(píng)定不同，毛細(xì)管電泳儀檢定規(guī)程［15］中的計(jì)量方法需要采用多級(jí)稀釋的方式配制溶液，這為檢出限不確定度的評(píng)定引入了更多的影響因素。筆者以毛細(xì)管電泳儀檢出限的不確定評(píng)定為例，以GUM評(píng)定過(guò)程為對(duì)照，詳細(xì)討論了Python在MCM上的具體應(yīng)用與分析。

1 GUM與MCM介紹

GUM是采用“不確定度傳播率”計(jì)算合成被測(cè)量估計(jì)值的測(cè)量不確定度評(píng)定方法。通過(guò)建立測(cè)量模型，分析不確定度來(lái)源，計(jì)算傳播系數(shù)，評(píng)定輸入量的標(biāo)準(zhǔn)不確定度，計(jì)算合成標(biāo)準(zhǔn)不確定度，確定擴(kuò)展不確定度等步驟完成不確定度的評(píng)定。

MCM為Monte Carlo Method 的簡(jiǎn)寫，是采用“概率分布傳播”的通用方法來(lái)評(píng)估被測(cè)量的測(cè)量不確定度。MCM要求建立測(cè)量模型，確定輸入量的概率密度函數(shù)（PDF），并通過(guò)隨機(jī)抽取每個(gè)輸入量的值，重復(fù)計(jì)算多次（通常為1×105到1×106次），分析計(jì)算結(jié)果，確定被測(cè)量的估計(jì)值及其包含區(qū)間。

2 實(shí)驗(yàn)部分

2.1 主要儀器與試劑

毛細(xì)管電泳儀：K1050型，北京凱奧科技發(fā)展有限公司。

維生素B6純度標(biāo)準(zhǔn)物質(zhì)：純度為99.9%，擴(kuò)展不確定度為0.4%(k=2)，編號(hào)為GBW(E) 100183，中國(guó)計(jì)量科學(xué)研究院。

2.2 標(biāo)準(zhǔn)溶液配制

對(duì)于2×10–6g／mL的維生素B6溶液的配制，JJG 964–2001 《毛細(xì)管電泳儀》檢定規(guī)程中給出了具體方法。

采用兩級(jí)稀釋的方式配制維生素B6溶液：首先，稱取10.0 mg溶于100 mL容量瓶中，配制1.0×10–4g／mL維生素B6溶液溶液。然后，取2 mL 1.0×10–4g／mL維生素B6溶液溶液稀釋到100 mL容量瓶中，配制成2.0×10–6g／mL維生素B6溶液。按式(1)計(jì)算維生素B6溶液的質(zhì)量濃度。

式中：ρ——維生素B6標(biāo)準(zhǔn)溶液的質(zhì)量濃度，g／mL；

m——維生素B6的質(zhì)量，mg；

P——維生素B6的純度；

V——溶液定容體積，mL；

1 000——單位換算系數(shù)。

2.3 毛細(xì)管電泳儀檢出限的測(cè)量

依據(jù)毛細(xì)管電泳儀檢定規(guī)程，采用質(zhì)量濃度為2×10–6g／mL的維生素B6溶液，連續(xù)進(jìn)樣3次測(cè)定峰高，按式(2)計(jì)算檢出限。

式中：D——檢出限，g／mL；

HN——噪聲峰高；

ρ——維生素B6標(biāo)準(zhǔn)溶液的質(zhì)量濃度；

H——樣品峰高。

采用《化學(xué)分析中不確定度的評(píng)估指南》［16］中的方法，按式(3)以相對(duì)不確定度計(jì)算傳播率：

3 GUM與MCM在不確定度評(píng)定中的應(yīng)用比較

以毛細(xì)管電泳儀檢出限的測(cè)量不確定度為例，對(duì)GUM與MCM在不確定度評(píng)定中的應(yīng)用進(jìn)行比較，結(jié)果列于表1。

表1 GUM與MCM不確定度評(píng)定應(yīng)用比較

續(xù)表1

MCM法不確定度評(píng)定的完整程序代碼如下。

以上代碼可直接輸入Python編譯器中運(yùn)行。注意代碼中的縮進(jìn)不可省略。

4 計(jì)算結(jié)果的分析

MCM的計(jì)算結(jié)果與GUM結(jié)果相同，表明MCM的評(píng)估方法具有很好的效果。采用Python進(jìn)行MCM計(jì)算分析，與傳統(tǒng)的GUM相比，代碼實(shí)現(xiàn)更加簡(jiǎn)單，且無(wú)關(guān)數(shù)學(xué)模型的復(fù)雜性，不需求導(dǎo)不確定度傳播系數(shù)，可直接通過(guò)計(jì)算獲得結(jié)果，計(jì)算速度快，計(jì)算精度高。

對(duì)計(jì)算結(jié)果做進(jìn)一步分析。圖1為MCM計(jì)算結(jié)果的概率分布。從圖1可以看出，毛細(xì)管電泳儀檢出限的計(jì)算結(jié)果具有正態(tài)分布的形態(tài)，符合線性模型下GUM分析的理論基礎(chǔ)，是對(duì)GUM分析結(jié)果的進(jìn)一步驗(yàn)證。由于計(jì)算結(jié)果具有對(duì)稱的分布形態(tài)，可以依據(jù)需求選取概率范圍，確定正確的k值用于計(jì)算擴(kuò)展不確定度。

圖1 檢出限MCM計(jì)算結(jié)果的PDF和CDF

圖2為計(jì)算結(jié)果與計(jì)算次數(shù)的關(guān)系展示。由圖2可以看出，隨著計(jì)算次數(shù)的增加，計(jì)算結(jié)果的平均值（檢出限）與標(biāo)準(zhǔn)偏差（標(biāo)準(zhǔn)不確定度）越接近理論計(jì)算值。當(dāng)計(jì)算次數(shù)大于1×104次試驗(yàn)時(shí)，計(jì)算結(jié)果足夠穩(wěn)定。更大數(shù)量的試驗(yàn)?zāi)軌蛱峁└交腜DF，得到更準(zhǔn)確的結(jié)果。

5 基于實(shí)際采樣值的MCM

以上是基于理想條件下的分析，直接將峰高的測(cè)量值理論化為正態(tài)分布，但實(shí)際的測(cè)量為有限次的測(cè)量，不可能獲取到正態(tài)分布內(nèi)的全部值，因此，需要分析在實(shí)際采樣條件下計(jì)算結(jié)果的分布情況。

在Python中提供了random.choices()的方法，即在已知采樣序列中隨機(jī)選取數(shù)值的方法。將11次測(cè)量的樣本作為研究對(duì)象，每次從中隨機(jī)取出3個(gè)數(shù)值作為一次試驗(yàn)的可能取值用于計(jì)算。將相關(guān)代碼更改如下：

圖2 檢出限MCM計(jì)算結(jié)果與計(jì)算次數(shù)的關(guān)系

經(jīng)過(guò)M=1×106次循環(huán)計(jì)算，得出計(jì)算結(jié)果：

檢出限（計(jì)算結(jié)果平均值）：1.13×10–8g／mL。

檢出限的標(biāo)準(zhǔn)不確定度（計(jì)算結(jié)果標(biāo)準(zhǔn)偏差）：1.83×10–9g／mL。

基于實(shí)際采樣計(jì)算結(jié)果與理論計(jì)算結(jié)果相同。

圖3為基于采樣值模擬計(jì)算后的PDF和CDF?？梢钥闯?，基于實(shí)際采樣進(jìn)行計(jì)算的結(jié)果，盡管得到的檢出限和標(biāo)準(zhǔn)不確定度與原計(jì)算相同，但計(jì)算結(jié)果的分布卻非常不同，不再完全符合正態(tài)分布。由于計(jì)算結(jié)果不再有對(duì)稱分布，則不能直接確定k值計(jì)算擴(kuò)展不確定度。擴(kuò)展不確定度應(yīng)采用區(qū)間的表示方式，結(jié)合累計(jì)概率密度分布函數(shù)(CDF)，通過(guò)百分位點(diǎn)獲取擴(kuò)展不確定度的上限和下限。

圖3 基于采樣值的MCM計(jì)算結(jié)果的PDF和CDF

實(shí)際上對(duì)從11次測(cè)量值隨機(jī)選取3個(gè)測(cè)量值這種方式做進(jìn)一步分析表明，經(jīng)過(guò)1×106次循環(huán)計(jì)算可以得到與正態(tài)分布相同的平均值和標(biāo)準(zhǔn)偏差。平均值與標(biāo)準(zhǔn)偏差相同，而計(jì)算結(jié)果的分布有巨大差異的原因是數(shù)據(jù)的有限性，即有限的數(shù)據(jù)未能正確地反應(yīng)正態(tài)分布特性。多數(shù)儀器的重復(fù)測(cè)量在理論上應(yīng)符合正態(tài)分布，但如果不符合正態(tài)分布，應(yīng)采用實(shí)際數(shù)據(jù)進(jìn)行計(jì)算，這樣才能獲得正確的不確定度評(píng)定結(jié)果。

6 結(jié)語(yǔ)

從計(jì)算結(jié)果可以看出，MCM能夠獲得與GUM相同的結(jié)果。但MCM不需計(jì)算傳播系數(shù)，不需精簡(jiǎn)模型，能夠獲得包含實(shí)際采樣信息在內(nèi)的更為全面準(zhǔn)確的計(jì)算結(jié)果。該方法基本涵蓋了Python軟件在不確定度計(jì)算時(shí)所需的全部代碼，方便在其它研究中參考使用。