新冠疫情预测建模记录
持续到现在的新冠疫情，已经改变了千万人的命运，其中也包括我
面对这么残暴的天灾，我做了一个预测模型，希望能贡献自己一点微薄的力量，让天灾早点过去
有首歌唱得好：
让我悲也好，让我悔也好，让疫情过去没烦恼；让我苦也好，让我累也好，让我看见大家的笑……
好了，废话结束，进入正题：
在这之前，github上已经有一个预测项目了：
https://github.com/YiranJing/Coronavirus-Epidemic-2019-nCov
于是我将他的源码down下来仔细看了看
该项目作者实现了自己的SIR和SEIR模型
（两种经典的传染病模型：https://baike.baidu.com/item/%E4%BC%A0%E6%9F%93%E7%97%85%E6%A8%A1%E5%9E%8B/5130035?fr=aladdin）
而模型的关键部分是定义了一个分布函数：
def Binomial_log_like(n: int, p: float, k: int):
loglik = lgamma(n+k+1) - lgamma(k+1) - lgamma(n+1) + k*np.log(p) + n*np.log(1-p)上述公式定义了对疫情传染人数趋势的一个假设公式
顾名思义，他将疫情的传染人数曲线看作二元对数分布
使用SIR和SEIR构建模型是否恰当呢？
这两种传染病动力学模型，其实都未考虑隔离措施的因素
都是建立在自然传播和自然康复的基础上，所以如果要做更客观的预测，需要考虑人类隔离活动的影响
而疫情数据实际的分布是否满足这种分布呢？我们取出历史数据看看：
国内每天新增病例趋势
从上图可以发现，新增病例曲线总体上更有点类似于伽马分布
而曲线中残缺的部分，猜测则是由于政府相关部门采取了强有力的隔离措施，导致了病毒传播的中断
所以，从宏观上看，新增病例曲线就是两种力量不断对抗的结果：
新增感染数量 = 病毒自然传播力 - 人类阻断摩擦力如果非要用一个物理现象做类比，就好像给往下滚动的雪球踩刹车
假设人类阻断力在23号封城以后是一个常量（已达到人类能实现的最高级别阻断）
那么其实我们就可以用伽马曲线近似的逼近23号之后的曲线
所以，这里我自定义了一个控制函数：
def control_curve(x, i0, missing, latency=7, infection_num=2):
"""
传染曲线，假设过了潜伏期即被隔离
今天传染人数 = (昨天实际传染人数 - 昨天隔离人数) * 平均传染人数
昨天隔离人数 = 七天前实际传染人数
今天传染人数 = (昨天实际传染人数 - 七天前实际传染人数) * 平均传染人数
:param x: 传染时间
:param i0: 初始传染人数
:param missing: 遗漏率
:param latency: 潜伏期（确切的说，这里代表平均发病时间，假设发病时间内每天都具有传染性）
:param infection_num: 平均传染人数
"""
if x >= latency:
isolated = i0 * (infection_num ** (x - latency))
else:
isolated = 0
day_num = i0 * missing * (infection_num ** x) - isolated
if day_num < 0:
day_num = 0
return day_num该函数中存在4个未知参数，我们可以利用scipy的优化器自动求解最优值：
from scipy.optimize import curve_fit
def control_curve_list(x, i0, missing, latency, infection_num):
y = list(map(control_curve, x, [i0] * len(x), [missing] * len(x), [latency] * len(x), [infection_num] * len(x)))
return y
popt,pcov = curve_fit(f=control_curve_list, xdata=x_data[:-1], ydata=y_data[:-1])通过 curve_fit 曲线拟合函数，我们就可以快速的得到一组近似最优解
拟合出来后，看看效果：
最优参数得到的模型预测120天结果
看起来很光滑呀，实际数据恐怕会比这个粗糙很多
那这个预测值，与实际是否匹配呢？
看起来，方差还是比较大的，有没有办法得到一个更精确的预测呢？
我们可以大胆假设，人类的隔离力度，在每个阶段、每个区域可能都不尽相同
每个时期的参数值可能会略有不同
所以接下来就用贝叶斯优化器对最优参数继续做更深入的探索：
这里先引入JS散度，方便贝叶斯优化器对结果进行比较
from bayes_opt import BayesianOptimization
def JS_divergence(p,q):
M=(p+q)/2
return 0.5*scipy.stats.entropy(p, M)+0.5*scipy.stats.entropy(q, M)再构建一个目标函数，返回可比较的分值
def result_function(nu, close_contacts, k, I0, E0, T, death_rate):
N = 14 * (10 ** 8)
econ = len(t)
R0 = I0 * (death_rate) * 2
try:
Est = Estimate_parameter(nu=nu, k=close_contacts, t=t, I=I)
eso = Estimate_Wuhan_Outbreak(Est, k=k, N=N, E0=E0, I0=I0, R0=R0, T=T, econ=econ)
result = eso._run_SIER('Forecat', 'China', 'Days after 2019-12-1', death_rate=death_rate, show_Plot=False)
score = JS_divergence(I,result['Infected'])
except Exception as e:
score = 999
print(e)
try:
score = float(score)
except:
score = 999
if score > 999:
score = 999
return -score然后，根据前面的最优参数，手工设置一些搜寻区间：
rf_bo = BayesianOptimization(
result_function,
{'nu': (1/20, 1/10),
'close_contacts': (5, 15),
'k': (1, 3),
'I0': (1, 10000),
'E0': (0,10000),
'T':(3,14),
'death_rate':(0.2,0.3)
}
)最后，把我们的优化器跑起来，让它自己找找
rf_bo.maximize(init_points=5,
n_iter=1000,)最终，我们就得到了一个相对靠谱的结果：
使用这个参数进行预测，与实际结果的差距就不大了
整套逻辑写好后，也可以把流程写成一个pipline，每天自动拟合
这样，就能不断的逼近最客观的结果！
以上就是我分享的思路，有不正确的地方还请大家多多批评！
