新冠疫情预测建模记录

简介: 持续到现在的新冠疫情,已经改变了千万人的命运,其中也包括我 面对这么残暴的天灾,我做了一个预测模型,希望能贡献自己一点微薄的力量,让天灾早点过去 有首歌唱得好: 让我悲也好,让我悔也好,让疫情过去没烦恼;让我苦也好,让我累也好,让我看见大家的笑…… 好了,废话结束,进入正题: 在这之前,github上已经有一个预测项目了: https://github.

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

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

有首歌唱得好:

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

好了,废话结束,进入正题:

在这之前,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构建模型是否恰当呢?

这两种传染病动力学模型,其实都未考虑隔离措施的因素

都是建立在自然传播和自然康复的基础上,所以如果要做更客观的预测,需要考虑人类隔离活动的影响

而疫情数据实际的分布是否满足这种分布呢?我们取出历史数据看看:

image
国内每天新增病例趋势

从上图可以发现,新增病例曲线总体上更有点类似于伽马分布

而曲线中残缺的部分,猜测则是由于政府相关部门采取了强有力的隔离措施,导致了病毒传播的中断

所以,从宏观上看,新增病例曲线就是两种力量不断对抗的结果:

新增感染数量 = 病毒自然传播力 - 人类阻断摩擦力

如果非要用一个物理现象做类比,就好像给往下滚动的雪球踩刹车

假设人类阻断力在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 曲线拟合函数,我们就可以快速的得到一组近似最优解

拟合出来后,看看效果:
image
最优参数得到的模型预测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,)

最终,我们就得到了一个相对靠谱的结果:
image

使用这个参数进行预测,与实际结果的差距就不大了

整套逻辑写好后,也可以把流程写成一个pipline,每天自动拟合

这样,就能不断的逼近最客观的结果!

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

相关文章
|
6月前
|
机器学习/深度学习 数据可视化 安全
Python随机森林、线性回归对COVID-19疫情、汇率数据预测死亡率、病例数、失业率影响可视化(下)
Python随机森林、线性回归对COVID-19疫情、汇率数据预测死亡率、病例数、失业率影响可视化
|
6月前
|
机器学习/深度学习 数据可视化 Python
Python随机森林、线性回归对COVID-19疫情、汇率数据预测死亡率、病例数、失业率影响可视化(上)
Python随机森林、线性回归对COVID-19疫情、汇率数据预测死亡率、病例数、失业率影响可视化
|
3月前
|
机器学习/深度学习 安全 算法
【2023年第十一届泰迪杯数据挖掘挑战赛】A题:新冠疫情防控数据的分析 32页和40页论文及实现代码
本文总结了2023年第十一届泰迪杯数据挖掘挑战赛A题的新冠疫情防控数据分析,提供了32页和40页的论文以及实现代码,涉及密接者追踪、疫苗接种影响分析、重点场所管控以及疫情趋势研判等多个方面,运用了机器学习算法和SEIR传染病模型等方法。
53 0
【2023年第十一届泰迪杯数据挖掘挑战赛】A题:新冠疫情防控数据的分析 32页和40页论文及实现代码
|
5月前
|
数据采集 机器学习/深度学习 数据可视化
使用决策树对金融贷款数据进行分析
使用决策树对金融贷款数据进行分析
|
6月前
R语言生存分析: 时变竞争风险模型分析淋巴瘤患者
R语言生存分析: 时变竞争风险模型分析淋巴瘤患者
R语言生存分析: 时变竞争风险模型分析淋巴瘤患者
|
6月前
|
机器学习/深度学习 安全 算法
Matlab基于SEIRD模型,NSIR预测模型,AHP层次分析法新冠肺炎预测与评估分析
Matlab基于SEIRD模型,NSIR预测模型,AHP层次分析法新冠肺炎预测与评估分析
|
6月前
|
机器学习/深度学习 数据挖掘
R语言逻辑回归模型的移动通信客户流失预测与分析
R语言逻辑回归模型的移动通信客户流失预测与分析
|
6月前
|
机器学习/深度学习 数据可视化 数据挖掘
R语言用CPV模型的房地产信贷信用风险的度量和预测
R语言用CPV模型的房地产信贷信用风险的度量和预测
|
6月前
|
算法 安全 数据可视化
用航空公司复杂网络对疫情进行建模
用航空公司复杂网络对疫情进行建模
|
6月前
|
安全 算法 数据可视化
在R语言中使用航空公司复杂网络对疫情进行建模
在R语言中使用航空公司复杂网络对疫情进行建模