基于隨機模擬探討《中國藥典》四分法取樣的注意事項

壽紀綱，霍志鵬，王玉，何毅

（1.天士力醫(yī)藥集團股份有限公司研究院，創(chuàng)新中藥關鍵技術國家重點實驗室，天津 300410；2.中國藥科大學 中藥學院，江蘇 南京 211198；3.天津中醫(yī)藥大學 中藥學院，天津 301617；4.天津大學 藥物科學與技術學院，天津 300072）

中藥材和飲片由于自然或人為原因不可避免地在外觀和成分含量方面存在差異，未必能保證批內(nèi)或批間質(zhì)量一致性，這給其抽樣檢測帶來了較大困難[1-2]。而無論是在中藥現(xiàn)代化研究中還是在藥企的實際生產(chǎn)中，采用科學的抽樣方法獲得能夠代表物料總體含量的樣本，是決定后續(xù)過程的結果是否具有科學性和可靠性的前提條件。2020年11月4日，國家藥監(jiān)局藥審中心發(fā)布了《中藥均一化研究技術指導原則（試行）》，孫昱等[3]指出了獲得均化物料必然包含的五個步驟，其中待均化物料的取樣檢測與均化效果確認兩個步驟均涉及混合飲片的抽樣評價，但目前尚無發(fā)現(xiàn)有明確提出可用于中藥均一化的混合飲片抽樣評價方法。

按照2020年版《中國藥典》四部通則0211 藥材與飲片取樣法規(guī)定，供檢驗用藥材或飲片的取樣方法為四分法，其具體操作為：當從藥材包件中抽取的樣品總量超過檢驗用量的數(shù)倍時，可按四分法再取樣，即將所有的樣品攤成正方形，依對角線畫“×”，使分成四等份，取用對角兩份；再如上操作，反復數(shù)次，直至最后剩余量能滿足供檢驗樣品用量；最終抽取的供檢驗用樣品量，一般不得少于檢驗所需用量的三倍即1/3 供實驗室分析用，另1/3 供復核用，其余1/3 留樣保存[4]?？梢钥吹?，四分法是一個逐步減小抽樣樣本量的過程，然而對于取樣樣本量大小與樣本組成的變異卻少見報道。

蒙特卡洛方法，也稱統(tǒng)計模擬方法，于二十世紀四十年代被“曼哈頓計劃”的成員John von Neumann 等人首先提出，是一種通過生成合適的隨機數(shù)和觀察一些服從特定性質(zhì)或?qū)傩缘臄?shù)據(jù)來解決問題的方法，可通過計算機進行統(tǒng)計抽樣試驗來提供近似解[5-6]。因此本文依據(jù)蒙特卡洛方法對四分法的取樣過程進行隨機模擬，首先按照四分取樣法的操作流程建立隨機抽樣模型，然后進行多輪模擬抽樣并以抽樣結果的頻數(shù)分布代表抽樣模型的概率基礎，最后通過對抽樣結果的頻數(shù)分布情況進行對比與分析，探討取樣量與抽樣結果變異大小的關系。

1 材料與儀器

Windows 10 系統(tǒng)計算機（處理器：Core i7-855U四核1.8 GHz，內(nèi)存：8G）；Python3.8；GraphPad Prism7。

白芍飲片（天津天士力現(xiàn)代中藥資源有限公司，2007002）；焦糖色（上海愛普食品工業(yè)有限公司）。

2 方法與結果

2.1 混合飲片抽樣模型的建立

取同一批白芍飲片1 800 g，其中900 g 用焦糖色染成黑色，代表一批白芍飲片，剩余900 g 代表另一批白芍飲片，將兩批飲片充分混合后作為抽樣總體，混合效果如圖1所示。抽樣樣本的大小依照《中國藥典》一部白芍飲片項下的【鑒別】、【檢查】、【含量測定】、【浸出物】設置為30 g。基于抽樣總體的平均片重（0.25 g/片）可以將實際抽樣過程轉化為從7 200 片隨機混合的飲片總體中抽120 片的過程。飲片取樣模型采用編程軟件Python 3.8 建立，建模流程如圖2所示。

圖1 染色白芍飲片與未染色白芍飲片混合效果圖Fig.1 Mixed effect of dyed and non-dyed pieces of Radix Paeoniae Alba

圖2 混合飲片抽樣模型運行流程圖Fig.2 Operation process of sampling model of mixed decoction pieces

將實際抽樣與模擬抽樣各100 次的結果的頻數(shù)分布情況進行對比，如圖3a所示。從圖中可以看到兩條曲線的形狀近似且均在未染色飲片占比為0.5 時累計頻數(shù)達到最大，這說明將混合飲片抽象化為均勻的抽樣單位后使用計算機進行模擬抽樣可以對實際抽樣的情況進行初步預測。

從抽樣結果的產(chǎn)生的過程分析，上述抽樣服從二項分布規(guī)律，基于二項分布的理論概率對此進行驗證[7]。將混合飲片抽樣模型的抽樣次數(shù)設為100 000，得到抽樣模擬的結果后使用GraphPad Prism 軟件分別對模擬抽樣結果和基于二項分布理論概率預測的結果繪制頻數(shù)分布曲線如圖3b所示。可以看到兩條頻數(shù)分布趨勢線基本一致，這說明對于兩批飲片混合得到的抽樣總體，當抽樣總體量與抽樣量相差較大時，抽樣過程符合二項分布的抽樣模型。

圖3 實際抽樣與模擬抽樣各100 次的結果的頻數(shù)分布情況進行對比Fig.3 Comparison of frequency distribution between actual sampling and simulated sampling

2.2 基于隨機模擬探討“四分過程”對抽樣結果的影響

基于隨機模擬探討《中國藥典》0211 藥材與飲片抽樣法中的“將所有的樣品攤成正方形，依對角線畫“×”，使分成四等份，取用對角兩份”過程對抽樣結果的影響。使用Python 3.8 編程語言分別建立四分過程抽樣模型與直接隨機抽樣模型。

建立四分過程抽樣模型思路如下：使用random庫中的隨機數(shù)函數(shù)生成一個長度為900 的數(shù)組并使每一個數(shù)代表染色飲片與未染色飲片的概率各為1/2，表示從兩種飲片1∶1 混合的無限大的抽樣總體中獲得一個包含900 個飲片的抽樣大樣本；同時將獲得的抽樣大樣本排列為30×30 的正方形平面并使用matplotlib 庫繪圖表示；然后從“正方形”的對角線處將平面分為四個區(qū)域且使每個區(qū)域內(nèi)的飲片數(shù)相等；最后將“正方形”上下兩個對角區(qū)域的飲片合并作為模擬四分過程抽樣的結果，建模流程如圖4所示。

圖4 四分過程抽樣模型運行流程圖Fig.4 Operation process of quartering process sampling model

直接隨機抽樣模型按照相同方式獲得抽樣大樣本，但將四分過程取消，而是直接隨機從抽樣大樣本中取出半數(shù)飲片作為抽樣結果。

將兩個抽樣模型的抽樣次數(shù)均設為100 000，得到抽樣模擬的結果后使用GraphPad Prism 軟件分別對它們繪制頻數(shù)分布曲線，如圖5a所示；運行四分過程抽樣模型后可同時得到四分過程繪圖的結果如圖5b所示。黑色圓點表示染色飲片，白色圓點代表未染色飲片；按照四分過程的步驟，從對角線處將“正方形”分為四個區(qū)域并保證每個區(qū)域中的飲片數(shù)相等，如圖5c所示。對照圖5c 中的白色區(qū)域，從圖5b 中取出對應位置的飲片作為四分過程的抽樣結 果。

從圖5a 中可以看到兩條頻數(shù)分布趨勢線基本一致，這說明當抽樣總體是隨機分布時，使用包含四分過程的取樣方法進行抽樣與直接隨機抽樣產(chǎn)生的結果是一致的。

圖5 抽樣模型軟件模擬圖Fig.5 Sampling model software simulation diagram

2.3 基于隨機模擬探討重復抽樣過程對抽樣結果的影響

基于隨機模擬探討《中國藥典》0211 藥材與飲片抽樣法中的“如上操作（四分過程），反復數(shù)次，直至最后剩余量能滿足供檢驗用樣品量”過程對抽樣結果的影響。

重復過程抽樣模型采用編程軟件Python 3.8 建立，使前一次抽樣得到的樣本作為下一次抽樣的抽樣總體并繼續(xù)按照相同的抽樣過程進行下一次抽樣，直至獲得最終樣本，建模流程如圖6所示。

圖6 重復過程抽樣模型運行流程圖Fig.6 Operation process of repeated sampling process sampling model

重復抽樣過程是一個逐步減小樣本量的過程，為探究逐步減少抽樣樣本的過程對抽樣結果的影響，對重復抽樣四次過程中初始樣本到最終樣本的均值變化情況進行分析。將重復過程抽樣模型的初始大樣本量設置為5 760，抽樣次數(shù)設為100 000，循化抽樣數(shù)依次設為1、2、3、4、5，即對包含5 760 個抽樣單位的抽樣大樣本進行反復抽樣分別獲得最終樣本量為2 880、1 440、720、360、180 的抽樣結果并繪制頻數(shù)分布曲線，如圖7所示；然后分別對上述抽樣結果落于期望值±5%區(qū)間內(nèi)的累計頻數(shù)進行統(tǒng)計，如表1所示。可以看到隨著抽樣樣本量減小，模擬抽樣結果的頻數(shù)分布區(qū)間逐漸增大；而抽樣結果落入期望值左右固定區(qū)間內(nèi)的累計頻數(shù)逐漸減小。

圖7 對包含5 760 個飲片的抽樣大樣本循環(huán)抽樣1 次、2次、3 次、4 次、5 次的抽樣結果對比Fig.7 The sampling results of one,two,three,four and five times of cyclic sampling for large sample of 5 760 decoction pieces

表1 對包含5 760 個飲片的抽樣大樣本循環(huán)抽樣不同次數(shù)時抽樣結果落于0.5±5%區(qū)間的累計頻數(shù)及其占比Tab.1 The cumulative frequency and proportion of the sampling results of 5 760 decoction pieces in the 0.5±5% with different times of cyclic sampling

對循環(huán)抽樣數(shù)不同時抽樣結果的變異進行考察，分別以2 次、3 次、4 次、5 次、6 次、10 次抽樣作為一輪并計算每輪結果的相對標準偏差（RSD），通過計算機模擬得到10 000 輪抽樣模擬的RSD 并對其小于5%的輪數(shù)占比進行統(tǒng)計，如表2所示?？梢钥吹皆诿枯喅闃哟螖?shù)不同時，10 000 抽樣模擬的RSD 小于5%的比例都會隨抽樣樣本量減小而減小，以5 次抽樣作為一輪為例，對循環(huán)抽樣數(shù)不同時10 000 輪抽樣模擬的RSD 的頻數(shù)分布情況進行作圖，如圖8a～e所示；當循環(huán)抽樣數(shù)不同時，隨每輪抽樣次數(shù)的增加，10 000 輪抽樣模擬的RSD小于5%的比例的變化規(guī)律不一致，分別對循環(huán)抽樣數(shù)為3、4、5 時，抽樣模擬的RSD 在每輪抽樣次數(shù)不同時的分布情況進行作圖，如圖8f～h所示。此外根據(jù)二項分布的數(shù)學模型可以推導得到抽樣量一定時理論RSD 的計算公式為RSD=σ（x）/μ（x）=[7]，計算得循環(huán)抽樣數(shù)為3、4、5時理論RSD 分別為3.726%、5.270%、7.453%，并將它們標注到圖8（f～h）中，可以看到在循環(huán)抽樣數(shù)一定時，隨每輪抽樣次數(shù)的增加，RSD 分布的擬合曲線近似為鐘形且分布范圍逐漸減小，其對稱軸逐漸趨近于理論RSD。

圖8 抽樣模擬的RSD 分布情況Fig.8 RSD distribution of sampling simulation

表2 每輪抽樣不同次數(shù)且不同循環(huán)抽樣數(shù)時10 000 輪抽樣模擬的RSD 小于5%的占比Tab.2 When the number of cyclic sampling and sampling number of each round is different，the proportion of RSD simulated by 10 000 rounds of sampling that is less than 5%

最后對重復抽樣方式與直接隨機抽樣方式進行對比。將2.2 中直接隨機抽樣模型的抽樣大樣本數(shù)與抽樣樣本數(shù)分別設置為5 760 與180，使之能夠通過一次抽樣從抽樣大樣本中直接獲得最終樣本；然后將重復抽樣模型與直接隨機抽樣模型的抽樣次數(shù)均設為100 000，得到它們抽樣模擬的結果后使用GraphPad Prism 軟件繪制頻數(shù)分布曲線，如圖9所示。從圖中可以看到兩條頻數(shù)分布趨勢線基本一致。這說明當抽樣總體是隨機分布且最終樣本量一定時，使用重復抽樣與直接隨機抽樣產(chǎn)生的結果是一致的。

圖9 重復過程模擬抽樣與直接隨機模擬抽樣的結果對比Fig.9 Comparison between results of repeated process simulation sampling and direct random simulation sampling

3 結論

本文按照《中國藥典》通則中的飲片取樣方法四分法的操作步驟使用計算機語言Python 建立了混合飲片的抽樣模型，然后基于蒙特卡羅方法得到了四分法模擬抽樣結果的概率分布情況，并對其進行了分析，旨在探討四分法的抽樣方式與抽樣量對抽樣結果代表性以及變異性的影響。相對于繁重的人工取樣統(tǒng)計，計算機模擬混合飲片取樣可在短時間內(nèi)獲得大量模擬數(shù)據(jù)，有助于快速了解抽樣的概率分布特征[8]。

由模擬結果分析發(fā)現(xiàn)，當抽樣總體是隨機分布且最終樣本量一定時，計算機模擬四分法抽樣與模擬直接隨機抽樣產(chǎn)生的結果是一致的，且四分法抽樣結果的代表性同樣會由于樣本量的減小而降低[9]。實際上，隨機分布狀態(tài)是一種完全混合狀態(tài)[10]，當物料未達到完全混合的狀態(tài)時，由于四分取樣法包含對中間樣本重新混合的過程，可能對減少抽樣變異有一定改善效果；在已經(jīng)充分混合的總體中，為使抽樣樣本具有充分的代表性，需要根據(jù)抽樣的具體情況評估合理的抽樣量大小。

考慮到變異性也是反映抽樣代表性的指標之一[11]，針對不同抽樣量時抽樣結果的RSD 分析，發(fā)現(xiàn)當抽樣樣本量減小時，樣本的RSD 增加，每輪抽樣模擬RSD 小于5%的概率會減?。划斆枯喅闃哟螖?shù)增加時，RSD 的分布范圍逐漸變小且其對稱軸逐漸趨近于理論RSD，這提示設計適當?shù)某闃哟螖?shù)才能獲得接近理論分布的RSD。以循環(huán)抽樣數(shù)為5 時為例，兩次抽樣作為一輪時，模擬抽樣的RSD 落在理論RSD 周圍的概率遠小于10 次抽樣作為一輪時，且兩次抽樣結果的均值實質(zhì)上有更大的變異，此時以2次抽樣的RSD 評價抽樣代表性顯然是不合理的。

由于中藥原料來源的復雜性，很多情況下無法得知混合飲片的具體情況，很難綜合判斷抽樣方法本身對抽樣結果產(chǎn)生的影響，從影響混合的主要因素入手對抽樣過程進行簡化或為一種可行的解決方法[12]。本研究發(fā)現(xiàn)實際抽樣與模擬抽樣的結果頻數(shù)分布曲線相似度較高，這表明隨機因素為影響物料混合的主要因素，因此本文從隨機混合的角度入手，對抽樣方式與抽樣量對抽樣結果代表性以及變異性的影響進行了初步推測，有助于為中藥均一化發(fā)展中更合理的混合飲片抽樣評價方法提供參考。