李春玉,安博文,王浩東,牛麗娜

（1.河北經(jīng)貿(mào)大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)學(xué)院,河北 石家莊 050061；2.華僑大學(xué) 經(jīng)濟(jì)與金融學(xué)院,福建 泉州 362021；3.新疆理工學(xué)院 理學(xué)院,新疆 阿克蘇 843000）

新型冠狀病毒肺炎現(xiàn)存確診人數(shù)是評(píng)估疫情發(fā)展形勢(shì)的重要指標(biāo),對(duì)于我國(guó)疫情防控工作開展具有重要意義.關(guān)于新冠病毒肺炎確診人數(shù)的研究多集中在傳播和預(yù)測(cè)方面.嚴(yán)閱等[1]率先基于時(shí)滯動(dòng)力學(xué)系統(tǒng)建立新冠病毒傳染病動(dòng)力學(xué)模型,通過(guò)反演模型參數(shù)模擬了湖北省疫情的發(fā)展趨勢(shì).白寧等[2]根據(jù)湖北省疫情發(fā)展相關(guān)數(shù)據(jù)建立動(dòng)力學(xué)模型,并通過(guò)MCMC 算法估計(jì)得到相關(guān)參數(shù),對(duì)湖北省未來(lái)疫情的發(fā)展趨勢(shì)作出評(píng)估預(yù)測(cè).林秋實(shí)等[3]結(jié)合湖北省確診病例的協(xié)變量信息,采用雙重穩(wěn)健模型對(duì)缺失數(shù)據(jù)進(jìn)行填補(bǔ),提高了動(dòng)力學(xué)模型的估計(jì)預(yù)測(cè)精度.黃森忠等[4]基于SEIR 模型厘清了新冠疫情的流行病學(xué)基本參數(shù),評(píng)估了疾控策略的效率.王霞等[5]構(gòu)建了離武漢較近城市的新冠疫情傳播網(wǎng)絡(luò)模型,分析了武漢及周邊地區(qū)復(fù)工復(fù)產(chǎn)的可能時(shí)間點(diǎn).郭尊光等[6]基于反應(yīng)擴(kuò)散方程采用最小二乘法估計(jì)得到新冠病毒傳染率,發(fā)現(xiàn)累計(jì)確診人數(shù)會(huì)隨傳染率和擴(kuò)散速率的增大而逐漸增加.丁志偉等[7]基于Back-projection模型推算得到我國(guó)部分地區(qū)新冠病毒肺炎的感染曲線,并給出了每日新增感染人數(shù)的期望值和置信區(qū)間.

在考察新冠疫情的傳播途徑與速度之余,充分了解當(dāng)前現(xiàn)存確診人數(shù)的分布狀況也同等重要.Ali等[8]首次將Poisson 分布引入傳染病人數(shù)研究中,白永昕等[9]構(gòu)造慢性病發(fā)病率和風(fēng)險(xiǎn)差時(shí)將傳染病確診人數(shù)假定服從Poisson 分布,王維賢等[10]也基于Poisson 分布構(gòu)造了慢性病相關(guān)差的置信區(qū)間.Poisson 分布是Negative Binomial 分布的極限分布, 意思就是在樣本量巨大的情況下,NegativeBinomial 分布與Poisson 分布是漸進(jìn)等價(jià)的[11].當(dāng)考察我國(guó)34 個(gè)地區(qū)新冠病毒肺炎確診人數(shù)的分布情況時(shí),樣本量為34,樣本量較少,此時(shí)假定現(xiàn)存確診人數(shù)服從Poisson 分布不免出現(xiàn)偏差,因此需要采用更為合適的Negative Binomial 分布.田茂再等[12]基于Negative Binomial 分布構(gòu)造了發(fā)病率的置信區(qū)間,并且模擬效果良好.之前學(xué)者研究新冠疫情的感染率時(shí)發(fā)現(xiàn),確診病例的感染率不盡相同,意味著Negative Binomial 分布中的感染率參數(shù)p 不再是一個(gè)恒定不變的常數(shù),而可能是一個(gè)服從某種統(tǒng)計(jì)分布的隨機(jī)變量.借鑒Lee 等[13]和Wang 等[14]的研究思路, 假定感染率p 服從Beta 分布, 那么現(xiàn)存確診人數(shù)便服從Beta-Negative Binomial 分布,本文就基于Beta-Negative Binomial 分布展開研究.

與此同時(shí),考慮到我國(guó)疫情防治工作效果明顯,部分地區(qū)現(xiàn)存確診人數(shù)持續(xù)為0 或持續(xù)為1,而且這樣的地區(qū)占比較高,這種情況與Lambert 等[15]提出的0 膨脹分布十分相近,也與Melkersson 等[16]提出的0-1 膨脹分布非常相近.從目前我國(guó)34 個(gè)地區(qū)現(xiàn)存確診人數(shù)的樣本數(shù)據(jù)來(lái)看,現(xiàn)存確診人數(shù)可近似看作一種膨脹分布.可將現(xiàn)存確診人數(shù)假定服從0 膨脹Beta-Negative Binomial 分布,也可將現(xiàn)存確診人數(shù)假定服從0-1 膨脹Beta-Negative Binomial 分布,對(duì)于分布中的參數(shù)可以借助EM算法進(jìn)行求解.Zhang 等[17]采用極大似然估計(jì)下的Fisher 信息矩陣和EM算法估計(jì)了0-1 膨脹分布的參數(shù);Tang 等[18]借助隱變量寫出了0-1 膨脹的似然函數(shù),并基于EM算法給出了參數(shù)的極大似然估計(jì)值,肖翔[19]又在此基礎(chǔ)上提供了0-1 膨脹參數(shù)估計(jì)的EM算法和M-H 抽樣算法.

綜合之前學(xué)者的研究成果,本文將基于Beta-Negative Binomial 分布研究我國(guó)34 個(gè)地區(qū)新冠病毒肺炎現(xiàn)存確診人數(shù)的分布情況.本文其余部分的布局如下：第1 部分依次構(gòu)建基于Beta-Negative Binomial的新冠肺炎感染人數(shù)分布、基于0 膨脹Beta-Negative Binomial 的新冠肺炎感染人數(shù)分布和基于0-1 膨脹Beta-Negative Binomial 的新冠肺炎感染人數(shù)分布;第2 部分對(duì)構(gòu)建的3 種分布進(jìn)行參數(shù)估計(jì),并給出擬合分布相對(duì)最優(yōu)的檢驗(yàn)方法;第3 部分對(duì)理論結(jié)果進(jìn)行數(shù)值模擬,分析估計(jì)算法的收斂性與精度;第4 部分對(duì)我國(guó)2021年1月1日～6日的新冠肺炎感染人數(shù)分布情況進(jìn)行擬合;第5 部分給出研究結(jié)論.

1 分布構(gòu)建

本部分針對(duì)新冠病毒肺炎現(xiàn)存確診人數(shù)構(gòu)建以下3 種分布,分別為：基于Beta-Negative Binomial 的新冠肺炎感染人數(shù)分布（簡(jiǎn)記為BNB 分布）、基于0 膨脹Beta-Negative Binomial 的新冠肺炎感染人數(shù)分布（簡(jiǎn)記為ZIBNB 分布）以及基于0-1 膨脹Beta-Negative Binomial 的新冠肺炎感染人數(shù)分布（簡(jiǎn)記為ZOIBNB 分布）.

用X 表示單個(gè)人的新冠病毒肺炎感染概率,當(dāng)現(xiàn)存確診人數(shù)Y 服從Negative Binomial 分布時(shí),則有YNB～NB (r ,x ),其概率質(zhì)量函數(shù)為

隨著現(xiàn)存確診人數(shù)變動(dòng)的復(fù)雜化,單一的Negative Binomial 分布并不能全面揭示全國(guó)34 個(gè)地區(qū)現(xiàn)存確診人數(shù)的分布狀況.若假設(shè)單個(gè)人肺炎病毒感染概率X 服從Beta 分布,即X～Beta(α,β),則對(duì)應(yīng)的概率密度函數(shù)為

此時(shí),現(xiàn)存確診人數(shù)Y 便服從具有分層性質(zhì)的離散分布,即Beta-Negative Binomial 分布,記為～BNB (r ,α,β),根據(jù)式（1）和式（2）可以寫出Y 與X 的聯(lián)合質(zhì)量函數(shù),為

通過(guò)對(duì)式（3）中的X 從0 到1 積分便可得到現(xiàn)存確診人數(shù)YBNB的概率質(zhì)量函數(shù),為

為刻畫較多地區(qū)現(xiàn)存確診人數(shù)為0 的情況,參考Tang 等[18]的研究方法,引入隱變量B～B(1,p ),令YZIBNB=YBNB(1?B),則現(xiàn)存確診人數(shù)Y 便服從0 膨脹Beta-Negative Binomial 分布,記為YZIBNB～ZIBNB (r,α,β,p ),其中p 表示0 膨脹參數(shù),概率質(zhì)量函數(shù)為

式（5）的具體證明過(guò)程可仿照文獻(xiàn)[20]中定理1 的證明.

為進(jìn)一步刻畫部分地區(qū)現(xiàn)存確診人數(shù)為1 的情形, 因此在0 膨脹Beta-Negative Binomial 分布的基礎(chǔ)上再引入隱變量C～B(1,q ),令YZOIBNB=YZIBNB+B(1,?C)則現(xiàn)存確診人數(shù)便服從0-1 膨脹Beta-Negative Binomial 分布,記為YZOIBNB～ZOIBNB (r,α,β, p ,q ),其中q 表示1 膨脹參數(shù),概率質(zhì)量函數(shù)為

式（6）的具體證明過(guò)程可仿照文獻(xiàn)[20]中定理2 的證明.下面將對(duì)本節(jié)建立的3 種分布進(jìn)行參數(shù)估計(jì).

2 參數(shù)估計(jì)

2.1 點(diǎn)估計(jì)

假設(shè)有n個(gè)地區(qū),現(xiàn)存確診人數(shù)的數(shù)據(jù)結(jié)構(gòu)為y=(y1,???,yn),決定0膨脹的隱變量B 數(shù)據(jù)結(jié)構(gòu)為b=(b1,???,bn),決定1膨脹的隱變量C數(shù)據(jù)結(jié)構(gòu)為c=(c1,???,cn),這里采用極大化似然法依次對(duì)上述3種分布的參數(shù)進(jìn)行估計(jì).

基于Beta-Negative Binomial 的新冠肺炎感染人數(shù)分布的對(duì)數(shù)似然函數(shù)為

由極大似然估計(jì)原理可知,式（7）最大化下的估計(jì)值α和β即為方程組（8）的解.

式（8）中：ψ(?)表示對(duì)數(shù)伽馬函數(shù)的導(dǎo)數(shù),即.經(jīng)過(guò)后續(xù)的數(shù)值模擬發(fā)現(xiàn),與同類計(jì)算非線性方程組的迭代算法比較（例如Newton-Raphso 迭代）,不動(dòng)點(diǎn)迭代更適用于求解式（8）,因此這里采用不動(dòng)點(diǎn)迭代計(jì)算式（8）的數(shù)值解.對(duì)式（8）進(jìn)行等價(jià)變形,有

的譜半徑.由不動(dòng)點(diǎn)定理可知,當(dāng)譜半徑ρ(φ′ ) <1時(shí),式（8）具有收斂的數(shù)值解,由此可以得到λ的取值.

最終計(jì)算出的數(shù)值解表示為

式中t 表示迭代次數(shù).

基于0 膨脹Beta-Negative Binomial 的新冠肺炎感染人數(shù)分布的對(duì)數(shù)似然函數(shù)為

式（9）中含有隱變量,因此采用EM算法求解似然的最大化.E 步計(jì)算隱變量B 的期望,為

式中隱變量B 數(shù)學(xué)期望的具體計(jì)算過(guò)程詳見文獻(xiàn)[20]中命題1 的證明.M 步計(jì)算似然函數(shù)最大化,式（9）最大化下的估計(jì)值α、β和 p分別為方程組（10）和方程（11）的解.

由式（11）解得

式（10）的求解方法與式（8）同理,最終的數(shù)值結(jié)果表示為

式中t 表示迭代次數(shù).

基于0-1 膨脹Beta-Negative Binomial 的新冠肺炎感染人數(shù)分布的對(duì)數(shù)似然函數(shù)為

式（12）中含有隱變量,因此采用EM 算法求解似然的最大化.E 步計(jì)算隱變量B 和隱變量C 的期望,為

式中隱變量B與隱變量C 數(shù)學(xué)期望的具體計(jì)算過(guò)程詳見文獻(xiàn)[20]中命題2 的證明.M 步計(jì)算似然函數(shù)最大化,式（12）最大化下的估計(jì)值α、β、p和 q分別為方程組（13）和方程（14）的解.

由式（14）解得

式（13）的求解方法與式（8）同理,最終的數(shù)值結(jié)果表示為

式中t 表示迭代次數(shù).

2.2 置信區(qū)間構(gòu)造

2.1 中計(jì)算出了分布參數(shù)的點(diǎn)估計(jì)值,由此可以分別寫出3 種分布不含隱變量的對(duì)數(shù)似然函數(shù).用n0表示現(xiàn)存確診人數(shù)為0 的地區(qū)個(gè)數(shù),用 n1表示現(xiàn)存確診人數(shù)為1 的地區(qū)個(gè)數(shù),為了方便利用Matlab 編程計(jì)算,分別記

基于Beta-Negative Binomial 的新冠肺炎感染人數(shù)分布的對(duì)數(shù)似然函數(shù)為

基于0 膨脹Beta-Negative Binomial 的新冠肺炎感染人數(shù)分布的對(duì)數(shù)似然函數(shù)為

基于0-1 膨脹Beta-Negative Binomial 的新冠肺炎感染人數(shù)分布的對(duì)數(shù)似然函數(shù)為

用l*表示各個(gè)分布的對(duì)數(shù)似然函數(shù),依次取為lBNB、lZIBNB和lZOIBNB,使用Matlab 軟件分別計(jì)算3 種分布對(duì)數(shù)似然函數(shù)對(duì)參數(shù)的二階偏導(dǎo),進(jìn)而求得3 種分布對(duì)應(yīng)的Fisher 信息矩陣.考慮如下形式的Fisher信息矩陣,將第（假設(shè)所有參數(shù)個(gè)數(shù)為 k,則i= 1,???,k）個(gè)參數(shù) oi對(duì)應(yīng)的Fisher 信息矩陣記為,其中,

由此,可以計(jì)算出回歸系數(shù) oi的方差為

從而,回歸系數(shù) oi的(1?α)置信區(qū)間為

式（15）中：zα2表示標(biāo)準(zhǔn)正態(tài)分布的分位點(diǎn).

2.3 分布比較

本文構(gòu)建了3 種新冠肺炎感染人數(shù)的分布,在實(shí)際情況中僅需要選擇一個(gè)相對(duì)最優(yōu)分布解釋現(xiàn)實(shí)問題.基于此,這里給出一種相對(duì)最優(yōu)分布的比較辦法.分布的比較原則為：當(dāng)兩個(gè)分布所含有的信息量差別不大時(shí),分布參數(shù)越少則分布越優(yōu).分布的信息量可以通過(guò)似然函數(shù)進(jìn)行反映,那么分布的比較就轉(zhuǎn)化為似然函數(shù)差別的比較,因此構(gòu)造似然比檢驗(yàn)（LR 檢驗(yàn)）統(tǒng)計(jì)量進(jìn)行假設(shè)檢驗(yàn)[11].

ZIBNB 分布與BNB 分布比較.原假設(shè)為兩個(gè)分布的信息量差別不大,即選擇BNB 分布;當(dāng)拒絕原假設(shè)時(shí),說(shuō)明ZIBNB 分布所含信息量要遠(yuǎn)高于BNB 分布,即ZIBNB 分布.構(gòu)造的LR 檢驗(yàn)統(tǒng)計(jì)量為

其中,kZIBNB表示ZIBNB 分布的參數(shù)個(gè)數(shù),kBNB表示BNB 分布的參數(shù)個(gè)數(shù).

ZOIBNB 分布與BNB 分布比較.原假設(shè)為兩個(gè)分布的信息量差別不大,即選擇BNB 分布;當(dāng)拒絕原假設(shè)時(shí),說(shuō)明ZOIBNB 分布所含信息量要遠(yuǎn)高于BNB 分布,即ZOIBNB 分布.構(gòu)造的LR 檢驗(yàn)統(tǒng)計(jì)量為

其中,kZOIBNB表示ZOIBNB 分布的參數(shù)個(gè)數(shù).

ZOIBNB 分布與ZIBNB 分布比較.原假設(shè)為兩個(gè)分布的信息量差別不大,即選擇ZIBNB 分布;當(dāng)拒絕原假設(shè)時(shí),說(shuō)明ZOIBNB 分布所含信息量要遠(yuǎn)高于ZIBNB 分布,即ZOIBNB 分布.構(gòu)造的LR 檢驗(yàn)統(tǒng)計(jì)量為

為了方便編程計(jì)算,分別記

則有

3 數(shù)值模擬

下面采用Matlab R2016a 編程進(jìn)行數(shù)值模擬.

3.1 Beta-Negative Binomial 分布

產(chǎn)生40 個(gè)隨機(jī)樣本.假設(shè)地區(qū)數(shù)為n=40;生成單個(gè)人肺炎病毒感染概率X～Beta(4,9),再生成Negative Binomial 分布隨機(jī)數(shù)Y～NB(8 ,x );最后生成現(xiàn)存確診人數(shù)YBNB～BNB(8,4,9).將現(xiàn)存確診人數(shù)YBNB和r=8作為已知條件進(jìn)行參數(shù)估計(jì),估計(jì)結(jié)果見表1.

表1 BNB 分布參數(shù)的數(shù)值模擬結(jié)果Tab.1 Numerical simulation results of BNB distribution parameters

表1數(shù)值模擬結(jié)果顯示,當(dāng)觀測(cè)樣本量為40 時(shí),估計(jì)值的相對(duì)誤差大小均控制在1%左右,且真實(shí)值都落入95%置信區(qū)間,說(shuō)明估計(jì)精度較高.圖1展示了參數(shù)估計(jì)值50 次迭代過(guò)程的軌跡圖,可以發(fā)現(xiàn),參數(shù)估計(jì)值均在迭代10 次之內(nèi)收斂,說(shuō)明估計(jì)結(jié)果的收斂性較強(qiáng).

圖1 BNB 分布參數(shù)估計(jì)的迭代軌跡圖Fig.1 Iterative trajectory diagram for parameter estimation of BNB distribution

3.2 0 膨脹Beta-Negative Binomial 分布

產(chǎn)生40 個(gè)隨機(jī)樣本.在3.1 的基礎(chǔ)上,生成決定0膨脹部分的隱變量B～B(1,0.3),進(jìn)而生成現(xiàn)存確診人數(shù)YZIBNB～ ZIBNB (8,4,9,0.3).將現(xiàn)存確診人數(shù)YZIBNB和 r=8作為已知條件進(jìn)行參數(shù)估計(jì),估計(jì)結(jié)果見表2.

表2 ZIBNB 分布參數(shù)的數(shù)值模擬結(jié)果Tab.2 Numerical simulation results of ZIBNB distribution parameters

表2數(shù)值模擬結(jié)果顯示,當(dāng)觀測(cè)樣本量為40 時(shí),估計(jì)值的相對(duì)誤差大小均控制在5%以內(nèi),且真實(shí)值都落入95%置信區(qū)間,說(shuō)明估計(jì)精度較高.圖2展示了參數(shù)估計(jì)值50 次迭代過(guò)程以及隱變量估計(jì)值49 次迭代過(guò)程的軌跡圖,可以發(fā)現(xiàn),分布參數(shù)和隱變量的估計(jì)值均在迭代10 次之內(nèi)收斂,說(shuō)明估計(jì)結(jié)果的收斂性較強(qiáng).

圖2 ZIBNB 分布參數(shù)和隱變量估計(jì)的迭代軌跡圖Fig.2 Iterative trajectory diagram for parameters and hidden variable estimatesof ZIBNB distribution

3.3 0-1 膨脹Beta-Negative Binomial 分布

產(chǎn)生40 個(gè)隨機(jī)樣本.在3.2 的基礎(chǔ)上,生成決定1膨脹部分的隱變量C～B(1,0.4),進(jìn)而生成現(xiàn)存確診人數(shù)YZOIBNB～ZOIBNB (8,4,9,0.3,0.4).將現(xiàn)存確診人數(shù)YZOIBNB和r=8作為已知條件進(jìn)行參數(shù)估計(jì),估計(jì)結(jié)果見表3.

表3 ZOIBNB 分布參數(shù)的數(shù)值模擬結(jié)果Tab.3 Numerical simulation results of ZOIBNB distribution parameters

表3數(shù)值模擬結(jié)果顯示,當(dāng)觀測(cè)樣本量為40 時(shí),估計(jì)值的相對(duì)誤差大小均控制在5%以內(nèi),且真實(shí)值都落入95%置信區(qū)間,說(shuō)明估計(jì)精度較高.圖3展示了參數(shù)估計(jì)值50 次迭代過(guò)程以及隱變量估計(jì)值49 次迭代過(guò)程的軌跡圖,可以發(fā)現(xiàn),分布參數(shù)和隱變量的估計(jì)值均在迭代20 次之內(nèi)收斂,說(shuō)明估計(jì)結(jié)果的收斂性較強(qiáng).

圖3 ZOIBNB 分布參數(shù)和隱變量估計(jì)的迭代軌跡圖Fig.3 Iterative trajectory diagram for parameters and hidden variable estimatesof ZOIBNB distribution

4 新冠病毒肺炎感染人數(shù)的分布擬合

這部分對(duì)我國(guó)34 個(gè)地區(qū)2021年1月1日～6日新冠病毒肺炎現(xiàn)存確診人數(shù)的分布狀況進(jìn)行擬合,樣本數(shù)據(jù)來(lái)源于國(guó)家及各省、市、自治區(qū)衛(wèi)健委,數(shù)據(jù)截止日期為2021年1月7日0 時(shí).

4.1 樣本數(shù)據(jù)描述

本文選取的數(shù)據(jù)指標(biāo)為新冠病毒肺炎現(xiàn)存確診人數(shù)（Y,單位：人）,作出2021年1月1日～6日每一天的現(xiàn)存確診人數(shù)分布狀況,如圖4所示.從1月1日～6日的數(shù)據(jù)比較來(lái)看,香港的現(xiàn)存確診人數(shù)最多,湖北、新疆、江西、重慶、甘肅、海南、吉林、貴州、寧夏、澳門、青海和西藏的現(xiàn)存確診人數(shù)持續(xù)為0,湖南、安徽和廣西的現(xiàn)存確診人數(shù)持續(xù)為1.其中,1月1日有13 個(gè)地區(qū)現(xiàn)存確診人數(shù)為0,有3 個(gè)地區(qū)現(xiàn)存確診人數(shù)為1,二者占比為47.06%;1月2日、3日、5日和6日有12 個(gè)地區(qū)現(xiàn)存確診人數(shù)為0,有4 個(gè)地區(qū)現(xiàn)存確診人數(shù)為1,二者占比為47.06%;1月4日有12 個(gè)地區(qū)現(xiàn)存確診人數(shù)為0,有3 個(gè)地區(qū)現(xiàn)存確診人數(shù)為1,二者占比為44.12%.我國(guó)34 個(gè)地區(qū)現(xiàn)存確診人數(shù)中0 值和1 值占比較高,這也間接證明了構(gòu)建膨脹分布描述現(xiàn)存確診人數(shù)分布狀況的必要性.

圖4 2021年1月1日～6日現(xiàn)存確診人數(shù)條形圖Fig.4 Bar chart of the number of existing confirmed cases from January 1 to 6,2021

4.2 估計(jì)結(jié)果分析

假定新冠病毒肺炎現(xiàn)存確診人數(shù)依次服從Beta-Negative Binomial 分布、0 膨脹Beta-Negative Binomial 分布和0-1 膨脹Beta-Negative Binomial 分布,根據(jù)前文推導(dǎo)的理論公式分別對(duì)3 個(gè)分布的參數(shù)進(jìn)行估計(jì),以下結(jié)果均采用Matlab R2016a 編程得到.

采用Matlab 計(jì)算數(shù)值解時(shí)發(fā)現(xiàn),若假定現(xiàn)存確診人數(shù)服從Beta-Negative Binomial 分布,Matlab 編程結(jié)果一直報(bào)錯(cuò),這就意味著無(wú)法得到數(shù)值結(jié)果,說(shuō)明不能找到Beta-Negative Binomial 分布對(duì)應(yīng)的相關(guān)參數(shù),即現(xiàn)存確診人數(shù)不會(huì)服從Beta-Negative Binomial 分布,也就是說(shuō)無(wú)法用Beta-Negative Binomial 分布描述現(xiàn)存確診人數(shù)的分布狀況.出現(xiàn)這一計(jì)算結(jié)果其實(shí)并不意外,從5.1 的數(shù)據(jù)結(jié)構(gòu)分析中就可以了解到,每一日現(xiàn)存確診人數(shù)中0 值和1 值的占比均超過(guò)40%,甚至其中有5 天時(shí)間現(xiàn)存確診人數(shù)接近50%,因此需要用膨脹的Beta-Negative Binomial 分布進(jìn)行擬合.當(dāng)然這也是本文所作的主要工作,將Beta-Negative Binomial 分布推廣到了0 膨脹Beta-Negative Binomial 分布和0-1 膨脹Beta-Negative Binomial 分布.

表4給出的是0 膨脹Beta-Negative Binomial 分布和0-1 膨脹Beta-Negative Binomial 分布的參數(shù)估計(jì)結(jié)果,并給出了這兩種分布的比較檢驗(yàn)結(jié)果.由LR 檢驗(yàn)結(jié)果可知,每一日的LR 檢驗(yàn)均在5%顯著水平下拒絕原假設(shè),即0-1 膨脹Beta-Negative Binomial 分布更適合描述我國(guó)新冠病毒肺炎現(xiàn)存確診人數(shù)的分布狀況.1月1日現(xiàn)存確診人數(shù)服從ZOIBNB（8,1.22,2.54,0.46,0.84）,1月2日現(xiàn)存確診人數(shù)服從ZOIBNB（8,1.22,2.60,0.46,0.77）,1月3日現(xiàn)存確診人數(shù)服從ZOIBNB（8,1.38,3.33,0.46,0.76）,1月4日現(xiàn)存確診人數(shù)服從ZOIBNB（8,1.27,2.82,0.43,0.82）,1月5日現(xiàn)存確診人數(shù)服從ZOIBNB（8,1.42,3.79,0.47,0.76）,1月6日現(xiàn)存確診人數(shù)服從ZOIBNB（8,1.29,3.34,0.46,0.76）.

表4 分布參數(shù)的估計(jì)結(jié)果Tab.4 Estimation results of distribution parameters

通過(guò)比較1月1日～6日的分布參數(shù)可以發(fā)現(xiàn),決定肺炎感染概率的參數(shù)α和β變化最大,連續(xù)6 天α變動(dòng)的標(biāo)準(zhǔn)差為8.32%、β變動(dòng)的標(biāo)準(zhǔn)差為49.46%;決定1 膨脹的參數(shù) q變化次之,連續(xù)6 天變動(dòng)的標(biāo)準(zhǔn)差為3.56%;決定0 膨脹的參數(shù) p變化最小,連續(xù)6 天變動(dòng)的標(biāo)準(zhǔn)差為1.37%.由此可見,新冠病毒肺炎感染概率每日波動(dòng)較大,但由于我國(guó)疫情防治工作一直有條不紊的開展,我國(guó)近一半地區(qū)現(xiàn)存確診人數(shù)持續(xù)為0 或持續(xù)為1,且波動(dòng)情況較小.

5 小結(jié)

本文將Beta-Negative Binomial 分布推廣到0 膨脹Beta-Negative Binomial 分布和0-1 膨脹Beta-Negative Binomial 分布,分別構(gòu)建新冠病毒肺炎現(xiàn)存確診人數(shù)BNB 分布、ZIBNB 分布以及ZOIBNB 分布;采用EM算法給出分布參數(shù)的點(diǎn)估計(jì)值,此部分工作的重點(diǎn)在于M步選用了不動(dòng)點(diǎn)迭代算法,之前學(xué)者多采用Newton-Raphso 迭代算法[20-21],但經(jīng)過(guò)數(shù)值模擬發(fā)現(xiàn),Newton-Raphso 迭代無(wú)法計(jì)算出最終結(jié)果,因此這里選用不動(dòng)點(diǎn)迭代,使迭代算法更適用于本文所提到的3 個(gè)分布;基于Fisher 信息矩陣計(jì)算估計(jì)參數(shù)的方差,進(jìn)而給出參數(shù)的區(qū)間估計(jì)值;基于LR 檢驗(yàn)給出BNB 分布、ZIBNB 分布和ZOIBNB 分布中相對(duì)最優(yōu)分布的遴選方法.通過(guò)數(shù)值模擬發(fā)現(xiàn),本文給出的估計(jì)算法收斂性很強(qiáng)、精度較高.新冠病毒肺炎感染人數(shù)分布的擬合結(jié)果顯示,ZOIBNB 分布更適用于描述我國(guó)34 個(gè)地區(qū)現(xiàn)存確診人數(shù)的分布情況,通過(guò)比較2021年1月1日～6日現(xiàn)存確診人數(shù)的分布參數(shù)發(fā)現(xiàn),單個(gè)人的每日病毒感染率波動(dòng)較大,我國(guó)有近一半地區(qū)現(xiàn)存確診人數(shù)為0 或?yàn)?.