壽紀(jì)綱，霍志鵬，王玉，何毅

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

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

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

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

1 材料與儀器

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

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

2 方法與結(jié)果

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

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

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

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

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

從抽樣結(jié)果的產(chǎn)生的過(guò)程分析，上述抽樣服從二項(xiàng)分布規(guī)律，基于二項(xiàng)分布的理論概率對(duì)此進(jìn)行驗(yàn)證[7]。將混合飲片抽樣模型的抽樣次數(shù)設(shè)為100 000，得到抽樣模擬的結(jié)果后使用GraphPad Prism 軟件分別對(duì)模擬抽樣結(jié)果和基于二項(xiàng)分布理論概率預(yù)測(cè)的結(jié)果繪制頻數(shù)分布曲線如圖3b所示?？梢钥吹絻蓷l頻數(shù)分布趨勢(shì)線基本一致，這說(shuō)明對(duì)于兩批飲片混合得到的抽樣總體，當(dāng)抽樣總體量與抽樣量相差較大時(shí)，抽樣過(guò)程符合二項(xiàng)分布的抽樣模型。

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

2.2 基于隨機(jī)模擬探討“四分過(guò)程”對(duì)抽樣結(jié)果的影響

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

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

圖4 四分過(guò)程抽樣模型運(yùn)行流程圖Fig.4 Operation process of quartering process sampling model

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

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

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

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

2.3 基于隨機(jī)模擬探討重復(fù)抽樣過(guò)程對(duì)抽樣結(jié)果的影響

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

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

圖6 重復(fù)過(guò)程抽樣模型運(yùn)行流程圖Fig.6 Operation process of repeated sampling process sampling model

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

圖7 對(duì)包含5 760 個(gè)飲片的抽樣大樣本循環(huán)抽樣1 次、2次、3 次、4 次、5 次的抽樣結(jié)果對(duì)比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 對(duì)包含5 760 個(gè)飲片的抽樣大樣本循環(huán)抽樣不同次數(shù)時(shí)抽樣結(jié)果落于0.5±5%區(qū)間的累計(jì)頻數(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

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

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

表2 每輪抽樣不同次數(shù)且不同循環(huán)抽樣數(shù)時(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%

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

圖9 重復(fù)過(guò)程模擬抽樣與直接隨機(jī)模擬抽樣的結(jié)果對(duì)比Fig.9 Comparison between results of repeated process simulation sampling and direct random simulation sampling

3 結(jié)論

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

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

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

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