2/5日更新：
2月4日匯總的“疑似病例”人數(shù)有點(diǎn)奇怪，有的地方顯示累計(jì)疑似病例23260，對比昨天只新增了46例。另一些地方顯示對比昨日新增3971例，那么就肯定不止23260例。
如果說累計(jì)疑似人數(shù)是對的，那么防控就起了很好的效果-->假定所有的疑似都轉(zhuǎn)成確診，人數(shù)也控制在4萬多了。如果新增疑似人數(shù)是對的，那么現(xiàn)在還在持續(xù)增加。
加入了2/4日的數(shù)據(jù)，預(yù)測最終確診人數(shù)在3.8萬人左右。

2/4日更新：
有幾天沒來了。起因在于我看了黃岡市唐主任一問三不知的視頻以后，意識到可能很多地方對于疫情人數(shù)的匯報(bào)不準(zhǔn)確。因此在整體意義上進(jìn)行預(yù)測意義就不那么大了。（借口）
經(jīng)過幾天線性增長以后，我把模型換成了Logistic模型，也就是累計(jì)病例會(huì)成為一個(gè)S形曲線。
這個(gè)曲線有一個(gè)上限，也就是最終累計(jì)確診人數(shù)，目前我的模型預(yù)測最終大約會(huì)有31240人確診。大約12天后趨近這個(gè)數(shù)值。

Figure_1.png

1/29日更新：
今天一起來就聽到好消息，1/28日新增確診人數(shù)為1459，少于1/27日的新增人數(shù)17xx人。不過新增疑似人數(shù)3248例仍在上升。（注意這里我討論的是，“每日增加人數(shù)”的變化，不是累計(jì)人數(shù)的變化。）
這樣就有兩種可能：

1. 新增確診人數(shù)下降，但新增疑似人數(shù)仍在上升，有一種可能性是確診人數(shù)是受到試劑盒制約而下降，而更多的感染者在“疑似人數(shù)”中無法被驗(yàn)證。這樣的話，接下來可能新增人數(shù)還會(huì)增加；
2. 假如沒有受到制約，1/29日真的是增量拐點(diǎn)，那么應(yīng)該有一個(gè)短期的拉鋸?？梢苑赐破骄鶟摲诩s為5天。
   要不要把指數(shù)模型換成logistic模型不好說，我決定觀察一天再定。數(shù)據(jù)本身的干擾因素太多，預(yù)測也不一定準(zhǔn)確，只能大概參考一下。

--------------------------------------差點(diǎn)把分割線同志忘記了---------------------------------------------------
1/28日更新：
真香~
事實(shí)證明，當(dāng)模型基本管用的時(shí)候，就不應(yīng)該修改。我前天預(yù)測的1/27日確診人數(shù)為4531，今天報(bào)道的確診人數(shù)為4515例，誤差只有0.3%，昨天修改之后，反而誤差更大。
我認(rèn)為1/26日的確診人數(shù)有問題，因?yàn)榉N種原因，少報(bào)了。
預(yù)測1/28日確診人數(shù)：6849

--------------------------------------我是一條分割線---------------------------------------------------
1/27日更新：
1/26日確診人數(shù)2744人，與預(yù)測2982人誤差8.6%，覺得主要有兩個(gè)原因：

1. 檢測試劑盒不夠 --> 確診人數(shù)受到限制
2. 床位不夠 --> 能留下來被檢測的人數(shù)受到限制
   修正1/27確診人數(shù)預(yù)測：3721

其他
現(xiàn)在比較顧慮的是1/23日武漢封城之前離開的大量人口，其中可能有攜帶者。按照7-10天潛伏期來算，確診人數(shù)可能在1/30-2/3日到達(dá)頂峰，之后預(yù)防措施見效，應(yīng)該就會(huì)下降了。
--------------------------------------我是一條分割線---------------------------------------------------
原文，寫于1/26日
最近新型冠狀病毒的爆發(fā)肯定牽動(dòng)了無數(shù)人的心。筆者在這里構(gòu)建了一個(gè)數(shù)學(xué)模型，希望能夠預(yù)測感染人數(shù)趨勢以及給大家提供一些防控指導(dǎo)。
太長不看：

1. 傳染率大約為1.5，也就是每個(gè)人會(huì)傳播1.5個(gè)人左右
2. 未來兩天的預(yù)測結(jié)果：
   1月26日確診人數(shù)：2982
   1月27日確診人數(shù)：4531
   由于樣本數(shù)很?。ㄖ挥?0天），同時(shí)早期數(shù)據(jù)可能存在較大偏差或誤報(bào)，模型有可能存在overfitting，還需要不斷修正。
3. 關(guān)于干預(yù)防控，目前的病例應(yīng)該都是23日武漢封城之前的，那時(shí)候的干預(yù)應(yīng)該還是線性的（影響k*x項(xiàng)）。1/30號以后的病例，應(yīng)該就是影響指數(shù)函數(shù)的系數(shù)q了，到時(shí)候的控制效果會(huì)好很多。

一、數(shù)據(jù)

首先搜集數(shù)據(jù)，來源是各大網(wǎng)站公開的國內(nèi)各省市確診人數(shù)和疑似病例人數(shù)（不包括外國）。

1月15日前沒有找到相關(guān)數(shù)據(jù)。筆者搜集到的數(shù)據(jù)見表1。

image.png

Date：日期
True_num：確診病例
Pos_num：疑似病例
Total：true_num + sus_num

二、模型

傳染病一般按照指數(shù)傳染，但是因?yàn)榘l(fā)現(xiàn)后大家有救治和預(yù)防，因此需要修正。在預(yù)防中，有一些本會(huì)感染的人沒有被感染，模型中假定這部分人數(shù)與時(shí)間成正比。

模型為：

image.png

其中：
· 第x天的感染人數(shù)Nx
· 疾病傳播率為q，q>1表示疾病在擴(kuò)張
· 第x天本應(yīng)感染而沒有感染的人數(shù)為k*x

三、確定病發(fā)起點(diǎn)（Zero-Day）

由于網(wǎng)上搜集到的數(shù)從1月16日開始，為了使指數(shù)模型更精確，需要確定病發(fā)起點(diǎn)，即Zero-Day.

在網(wǎng)上查到了以下資料：
公開消息顯示，疫情初期（12月8日），由武漢宣傳通稿爆出正式收治第一例感染病人，當(dāng)時(shí)網(wǎng)絡(luò)已爆出有多個(gè)疑似感染個(gè)案，正在入院檢查及治療，期間有感染者處于危重狀態(tài)。二十天后的12月28日，武漢市衛(wèi)建委勉為其難通報(bào)疫情信息，稱已有27例收治記錄，初步認(rèn)定為病毒性肺炎，并由專家強(qiáng)調(diào)“未發(fā)現(xiàn)有人傳人現(xiàn)象”。

由于12月8日當(dāng)時(shí)的檢測手段有限，感染人數(shù)不大可靠，因此本人采用了12月28日的數(shù)據(jù)，在此模型下，1月16號的為第20個(gè)數(shù)據(jù)節(jié)點(diǎn)。

因此現(xiàn)在的數(shù)據(jù)見表2，其中date_id欄為發(fā)現(xiàn)疾病后第date_id天。

image.png

四、程序

因?yàn)楝F(xiàn)在有些疑似病例因?yàn)槿狈υ噭o法確診，因此在模型中，我既預(yù)測了確診病例"Confirmed Num of Patients"，也計(jì)算了全部病例"Total Num of Patients"（確診+疑似）
得出的對于全部病例的擬合曲線如下：

Figure_1.png

import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np
from scipy.optimize import curve_fit
import math

def exp_with_control(x_list, a, q, k):

    result_list = []
    for x in x_list:
        exp_part = a * (1 - q ** x) / (1 - q)
        i = 1
        while i < x:
            sum_q = 0
            sum_q = sum_q + (i * q ** (x - 1 - i))
            i += 1
        poly_part = k * sum_q * x
        result_list.append(exp_part - poly_part)

    return result_list

# import
data = pd.read_csv('data.csv')
date_id = [i for i in range(20, 30)]
data['date_id'] = date_id

# curve fit
popt_ec, pcov_ec = curve_fit(exp_with_control, xdata=date_id, ydata=data['total'], p0=[1.5, 1.5, 1.5])
y_hat_ec = exp_with_control(date_id, popt_ec[0],popt_ec[1],popt_ec[2])

# plot results

plt.plot(date_id, data['total'], label="Total Num of Patients")
plt.plot(date_id, y_hat_ec,'r--',label="Expected Total")
print("Coefficients", popt_ec)
RMSE_ec = math.sqrt(np.mean((data['total'] - y_hat_ec)**2))
print("RMSE_total", RMSE_ec)
plt.xlabel("Days")
plt.xlabel("Num of Patients")
plt.legend()
plt.show()

對于確診病例的擬合曲線：

Figure_2.png

這部分程序如下：

popt_ec_true, pcov_ec_true = curve_fit(exp_with_control, xdata=date_id, ydata=data['true_num'], p0=[1.5, 1.5, 1.5])
y_hat_ec_true = exp_with_control(date_id, popt_ec_true[0],popt_ec_true[1],popt_ec_true[2])
plt.plot(date_id, y_hat_ec_true,'r--', label="Confirmed Num of Patients")
print("Coefficients", popt_ec_true)

RMSE_ec_true = math.sqrt(np.mean((data['true_num'] - y_hat_ec_true)**2))
print("RMSE Confirmed", RMSE_ec_true)

plt.plot(date_id, data['true_num'], label="Expected Confirmed")
plt.legend()
plt.show()

五、討論

1. 可以看出，模型對于確診病例的預(yù)測誤差很小，遠(yuǎn)好于對于全部病例的預(yù)測。對于全部病例的誤差，可能是由于疑似病例中很多不是新型肺炎的人誤報(bào)導(dǎo)致。
2. 傳染率大約為1.5，也就是每個(gè)人會(huì)傳播1.5個(gè)人左右
3. 未來兩天的預(yù)測結(jié)果：
   1月26日確診人數(shù)：2982
   1月27日確診人數(shù)：4531
4. 關(guān)于防控干預(yù)，目前的病例應(yīng)該都是23日武漢封城之前的，所以那時(shí)候的干預(yù)應(yīng)該還是線性的（影響k*x項(xiàng)）。1/30號以后的病例，應(yīng)該就是影響指數(shù)函數(shù)的系數(shù)q了，那時(shí)候的控制效果會(huì)好很多。
5. 由于樣本數(shù)很?。ㄖ挥?0天），同時(shí)早期數(shù)據(jù)可能存在較大偏差或誤報(bào)，模型有可能存在overfitting，還需要不斷修正。