持续到现在的新冠疫情,已经改变了千万人的命运,其中也包括我
面对这么残暴的天灾,我做了一个预测模型,希望能贡献自己一点微薄的力量,让天灾早点过去
有首歌唱得好:
让我悲也好,让我悔也好,让疫情过去没烦恼;让我苦也好,让我累也好,让我看见大家的笑……
好了,废话结束,进入正题:
在这之前,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,每天自动拟合
这样,就能不断的逼近最客观的结果!
以上就是我分享的思路,有不正确的地方还请大家多多批评!