本文通信作者李志宏，為中國原子能科學研究院核物理研究所研究員，主要研究方向：核天體物理。
本站（www.paper.edu.cn）歡迎更多疫情相關學術投稿，將第一時間審核發布，免費分享，一起攜手，共克時疫！
本文為預印本論文，尚未經同行評議
新型冠狀病毒傳播的高斯分布模型
李星君1，李志宏2
1 北京化工大學化學工程學院，北京 100029
2 中國原子能科學研究院核物理研究所，北京 102413
摘要：本文利用三種數學分布模型擬合了我國多個地區新型冠狀病毒肺炎的確診病例數據，結果發現高斯分布模型給出的疫情發展曲線與實際數據的吻合度較高。高斯分布模型不但可以確定新增病例增長曲線的拐點時間、病毒的感染期和最大累計確診病例數，而且能夠很好地預測疫情的走勢，為國家的抗疫決策提供理論支撐。把高斯分布模型和二項式模型相結合起來，我們還導出了這次新型冠狀病毒肺炎的基本繁殖率(R0) 在1.68 至3.56之間，結果對該病毒的傳播機制研究有重要意義。
關鍵詞：流行病與衛生統計學，新型冠狀病毒肺炎，疾病傳播，基本繁殖率
基金項目： 中國原子能科學研究院院長基金項目，穩定性支持項目：WDJC-2019-13
0 引言
對疫情發展狀況的科學推斷不但是傳染病預測、預警和抗擊疫情決策關鍵依據，也是減小恐慌心理，消滅謠言的重要利器。因此，每逢疫情到來的時候，都會有大量疫情預測文章出現。2019年末從湖北武漢開始的新型冠狀病毒肺炎是一種雖然致死率低，但傳播非常迅速的惡性傳播疾病。截至2020年3月2日24時，全國共確診新型冠狀病毒肺炎80175例，疑似715例，死亡2915例。其傳播速度超過了2003年的SARS。對該疫情進行模型分析，有助於我們了解該新型冠狀病毒的傳播特點和發展規律，對此類病毒的防治和研究都具有重要意義。
人們對流行病傳播的數學模型研究迄今已經有一百多年的歷史了。1911年Ross[1, 2]首先建立了流行病傳播的數學模型，奠定了傳染病研究的動力學基礎。經過一百多年的發展，傳染病模型已經研究得相當細緻。在當前流行的模型中，主要有SI[3]、SIR[4]、SIRS[5]、SEIR[6] 等。隨著模型考慮的因素越多，模型越複雜，擬合參數也就越多，導致模型預測的誤差很大。比如SEIR 模型，其待定參數高達7個以上，要想得到比較好的模型參數至少需要十餘天數據累計。如果傳染病的周期較短，用這種較多參數的模型做預測就顯得很力不從心，甚至會出現在沒有獲得理想擬合參數的情況下，疫情已趨於結束。因此，對於傳染病的預測，最好能找到一種擬合參數少、預測效果可靠的數學模型。
在作者2020年2月7日完成的「新型冠狀病毒的傳播模型研究」[7]論文中，我們使用感染人數的增長速率為二項式分布的模型預測了全國和北京地區的疫情發展趨勢。從當前的進展情況來看，我們的預測結果總體表現良好：北京地區到2月20日時的新增確診病例小於5例/日，累計最大確診病例數為380例，與實際情況基本符合；全國疫情發展的趨勢也基本符合預測，但是受湖北地區2月12日起執行新確診標準的影響，作者預測的最大病例數有較大的出入。按照新的診斷標準，全國的最大確診病例數比以前的標準約增大1.6倍。因此，2月7日的預測結果換算成新檢測標準後應該是最大累計確診病例為7.8萬人。需要指出的是，我們原來的估算是建立在孤立模型的基礎上的。如果孤立系統因擴散或確診病例的流入事件被打破，就需要對原有的模型進行修正。修正的方法之一是把來自孤立系統之外的引入病例去掉，比如北京地區2020年3月8日的累計確診病例數為428例，其中就包括了13例國外輸入和23例外地來京病例。除去這些外部引入病例，結果和我們2020年2月7日預測的最大累計確診病例380 例[7]基本一樣。這樣做當然是可以的，不過對照起來不夠直觀，也會引起不必要的誤解。為了避免這種情況，可以考慮進行一下模型改進，使其包含外部引入因素。
在這篇文章中，我們嘗試把新增病例的二項式分布模型改為高斯分布模型，並用其解釋本次國內疫情較重的幾個地區的確診病例數據，評估利用高斯分布模型進行疫情預測的可行性。
1 高斯分布模型
對國內各地的疫情數據分析發現，當有外部輸入確診病例的情況出現時，會導致新增病例數的統計分布出現拖尾現象。主要表現在快到平臺期的新增病例數目高於預測結果(如圖1)，導致模型預測的累計確診病例和新增病例都存在一定程度的偏低現象。為了解決這個問題，需要進行模型改進。
圖1: 北京地區每日新增確診病例的二項式分布與最佳擬合分布的比較
首先，我們來挖掘一下二項式分布的內涵，把文獻[7]中的公式(2)做進一步的變換，令r=kM，則有：
此方程正是通常所說的Logistic統計模型，其回歸運算常用於流行病學的研究。通過logistic回歸分析，人們可以大致了解導致病患的危險因素。Logistic模型中的因變量為二項式分布，其增長率曲線為拋物線。方程(1)的解為：
該解就是我們平時所說的S型曲線，又叫皮爾曲線、理察曲線或者增長曲線等[8]。它描繪了傳染病的發生和發展過程：初始階段緩慢發展，接著是快速的增長階段，最後是平穩的發展時期一直到飽和狀態。可以看出，在式(1)中M為最大累計確診病例數，它等於r/k。也就是說，當達到飽和狀態時，孤立體系的所有人最終都將被感染的。由於接觸時間，患病時間的隨機性，導致新增確診病例隨時間的變化曲線接近高斯分布。因此有：
其中，M為最大累計確診病例數；μ為發病時間的數學期望值，它代表新增病例數的處於拐點的時間；σ為高斯分布的標準差，它決定了高斯分布的幅度，也和疫情的持續時間相關。我們猜測，2σ比較接近病毒的有效傳染期(潛伏期)，這可以通過各地的疫情數據加以確認。
當然，在新增病例分布偏離對稱的情況下，也可以考慮向長時間方向拖尾的指數分布模型：
上述三種分布模型都能擬合2003年非典的疫情的數據(二項式分布模型見引文[7]中的圖4，高斯分布和指數分布模型見本文圖2)，但它們的預測能力、穩定度，以及能否正確地反映疫情走向，需要進行檢驗。
圖2: 高斯分布和指數分布擬合2003年北京市非典累計和新增確診病例數據的結果。
我們使用三種模型分別擬合了北京、上海、廣東、浙江、河南、安徽和江西等地區的確診病例數據，並比較了使用不同數據點的擬合優度。結果列於表1 中。其中M1 代表多項式模型，M2 代表高斯分布模型，M3 代表指數分布模型，表中的日期代表選用該日期之前的數據進行擬合，擬合優度的表達式為：
其中，n為擬合使用的樣本數目，Ni為第i天的累計確診病例，Nmi為模型計算的第i天的累計確診病例，σi為第i天累計確診病例的統計誤差。
圖3: 高斯分布擬合北京市2020年新型冠狀病毒肺炎累計和新增確診病例的結果
從表1中的擬合優度數據可以看出：指數分布模型除在浙江省的數據擬合中表現稍好外，其它地區都差於二項式分布模型與高斯分布模型。這說明指數模型不能反映疫情發展趨勢的細節，因此可以棄之不用。二項式分布模型在北京、上海的確診病例擬合中表現較好，而高斯分布模型在其它地區的表現佔優。圖3為用高斯分布擬合北京地區累計和新增確診病例的結果，從中可以看出，高斯分布可以重現新增病例分布的尾部。擬合得到北京地區的最大累計確診病例數為429例，也與北京市近期公布的累計確診病例數符合很好。從結果看，高斯分布的模型計算結果已包含了外部引入病例情況，可以與官方公布的數據直接進行比對。
表1: 三種模型擬合結果的比較
其中的數值代表模型預測曲線相對於總數據集的擬合優度，M1指多項式模型，M2為高斯分布模型，M3代表指數分布模型。日期表示僅擬合該日期之前的病例數據。
下面來分析一下高斯分布模型的預測效果。我們假定使用前期疫情數據獲得的擬合優度χ2值比使用全部數據的擬合結果大十倍的預測效果是可以接受的。從表1 中可以看出，高斯分布模型在2月6日左右就可以預測大部分地區確診病例數據。模型計算與公布數據在全時間段內都符合很好，這表明模型參數相對穩定，全國的疫情防控工作總體向好，模型具有很好的預測效果。湖北地區的情況比較複雜，其病號太多，醫療資源不夠充分，感染者的確診時間可能延遲於感染時間。另外，檢測標準的改變與核減病例等不確定性也會影響模型的預測效果。我們因而使用改變標準前後的不同數據給出兩種模型計算結果，如圖4所示。新舊兩種標準之間的最大累計病例數目之比為1.68，新增確診病例數的拐點相差3天，分別對應於2 月6日和2月9日。如果不出現鬆懈和意外感染事件發生，湖北地區的累計確診病例平臺也會在3月中旬到來，最大的累計確診病例數為6.6萬人。本估計比王霞等人[9]的3月2日要晚一點。考慮到病毒的潛伏期，等到最後一位病人的所有密切接觸者結束14 天的隔離，疫情才算真正結束。
圖4: 湖北地區累計確診病例數據與模型計算結果的比較
表2列出了使用高斯分布模型擬合出的北京、廣東等地區新型冠狀肺炎確診病例的模型參數。其中，M代表最大累計確診病例，σ代表高斯分布的標準差，μ代表發病時間的數學期望值。模型給出的最大累計確診病例數很好地符合了當前各地區公布的累計確診病例的數值，各地方差參量的平均值為6.92，該數值的2倍非常吻合新型冠狀病毒14 天的感染期。
表2: 利用高斯分布模型擬合得到的國內幾個地區的高斯分布模型參數
2 結論
本文分析了三組參數較少的新型冠狀病毒肺炎疫情的預測模型，結果發現：高斯分布模型可以很好地符合新增確診病例和累計確診病例的發展曲線，能夠確定各地病毒的感染期，新增病例曲線的拐點時間，並可以預測最大累計確診病例數。基於這些特點，可以根據疫情前期的數據來預測疫情的發展趨勢。
二項式分布模型和高斯分布模型結合起來可以更好地解釋疫情傳播特徵。二項式分布模型可以給出新增病例的拐點時間，最大累計確診病例數以及病毒感染的加速因子r。高斯分布可以給出最大累計確診病例數，新增確診病例的拐點時間和病毒的傳染周期。感染加速因子與傳染期的乘積即為傳染病的基本繁殖率R0，它是人們研究傳染病傳播的最重要的指標之一。本工作得到的這次新型冠狀病毒肺炎的基本繁殖率的取值範圍在1.68至3.56 之間。山東省的數值最小為1.68，湖北省的最大為3.56。本結果與趙等人[10]等人利用時間序列模型研究的結果(R0的範圍為2.24-3.58)符合很好。
致謝
感謝中國原子能科學研究院核物理研究所的領導對本工作的鼓勵和大力支持；感謝中國原子能科學研究院院長基金的資助。
參考文獻（References）
[1] Ross R. The prevention of malaria[M]. London: John Murray; 1911, 651-686.
[2] Ross R. Some quantitative studies in epidemiology[J] Nature. 1911, 87:466-467.
[3] Bailey N T J, The Mathematical Theory of Infectious Diseases[M], 2nd ed. Macmillan, New York, 1975
[4] Luz P M, Struchiner C J, Galvani A P, et al. Modeling Transmission Dynamics and Control of Vector-Borne Neglected Tropical Diseases[J] . PLoS Neglected Tropical Diseases, 2010, 4(10):e761.
[5] Busenberg S and Cooke K, Vertically Tmnsmitted Diseases: Models and Dynamics[M], Springer, Berlin, (1993).
[6] Audrey M. Dor_elien, Ballesteros S, Grenfell B T. Impact of Birth Seasonality on Dynamics of Acute Immunizing Infections in Sub-Saharan Africa[J]. PLOS ONE, 2013, 8.
[7] 李星君，李志宏. 新型冠狀病毒的傳播模型研究[EB/OL]. 北京：中國科技論文在線[2020-02-12]. http://www.paper.edu.cn/releasepaper/content/202002-42.
[8] Kucharavy D, Guio D G. Application of S-shaped curves[J]. Procedia Engineering, 2011,
9:559-572.
[9] 王霞，唐三一，陳勇等，新型冠狀病毒肺炎疫情下武漢及周邊地區何時復工？[J]. 數據驅動的網絡模型分析.DOI://10.1360/SSM-2020-0037, 2020.
[10] Zhao S, Lin Q, Ran J et al. reliminary estimation of the basic reproduction number of novel coronavirus (2019-nCoV) in China, from 2019 to 2020: A data-driven analysis in the early phase of the outbreak[J]. Int J Infect Dis. 2020, doi: 10.1016/j.ijid.2020.01.050.