新冠疫情预测建模记录

持续到现在的新冠疫情，已经改变了千万人的命运，其中也包括我

面对这么残暴的天灾，我做了一个预测模型，希望能贡献自己一点微薄的力量，让天灾早点过去

有首歌唱得好：

让我悲也好，让我悔也好，让疫情过去没烦恼；让我苦也好，让我累也好，让我看见大家的笑……

好了，废话结束，进入正题：

在这之前，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，每天自动拟合

这样，就能不断的逼近最客观的结果！

以上就是我分享的思路，有不正确的地方还请大家多多批评！