python模拟预测一下新型冠状病毒肺炎的数据

发布时间:2020-02-05 22:07 来源:互联网 当前栏目:网站服务器

大家还好吗?

背景就不用多说了吧?本来我是初四上班的,现在延长到2月10日了。这是我工作以来时间最长的一个假期了。可惜哪也去不了。待在家里,没啥事,就用python模拟预测一下新冠病毒肺炎的数据吧。要声明的是本文纯属个人自娱自乐,不代表真实情况。

采用SIR模型,S代表易感者,I表示感染者,R表示恢复者。染病人群为传染源,通过一定几率把传染病传给易感人群,ta自己也有一定的几率被治愈并免疫,或死亡。易感人群一旦感染即成为新的传染源。

模型假设:

①不考虑人口出生、死亡、流动等情况,即人口数量保持常数。

②一个病人一旦与易感者接触就必然具有一定的传染力。假设 t 时刻单位时间内,一个病人能传染的易感者数目与此环境内易感者总数s(t)成正比,比例系数为β,从而在t时刻单位时间内被所有病人传染的人数为βs(t)i(t)。

③ t 时刻,单位时间内从染病者中移出的人数与病人数量成正比,比例系数为γ,单位时间内移出者的数量为γi(t)。
模型为



其中,β为感染系数,代表易感人群与传染源接触被感染的概率。γ为隔离(恢复)系数,我们对其倒数1/γ更感兴趣,代表了平均感染时间(average infectious period)。S(0)为初始易感人数,I(0)为初始感染人数。

按照[1]里面的代码模型的感染人数是这样的


现在的问题就是利用现有的数据找到新冠肺炎的β值,γ值等数据了。先把数据拔下来吧。从[3]上扒数据,由于数据不多,就手工完成吧。保存到csv文件里。

然后把数据作图


还有一个指标是再生数R0=β/γ,大于1时人群中大部分才被感染[4]。世卫组织1月23日的估计是R0在1.4到2.5之间[5],最新的根据前425例发病数据的估计值为2.2[6]。

文章[7]中的按一般病毒性肺炎恢复期25天计算得到的γ值为0.04。

关于β值和初始易感人群,[7]的作者采用的方法是先估计一个区间,然后用最小二乘法找到最佳参数,β≈3.57*10^-5。S[0]的范围为5000-30000人。[7]文章里有matlab代码,我用python改写一下,由于对最小二乘法法的实现比较陌生,尝试了半天,最后我决定用最笨的办法——穷举法。就是用两个嵌套循环将范围内所有β值和S0值都试一遍,计算每次尝试结果与实际数据之间差值的平方和,平方和最小的一组β值和S0值用来做预测。代码如下:

γ值设定为0.04,即一般病程25天

用最小二乘法估计β值和初始易感人数

gamma = 0.04
S0 = [i for i in range(20000, 40000, 1000)]
beta = [f for f in np.arange(1e-7, 1e-4, 1e-7)]
# 定义偏差函数
def error(res):
 err = (data["感染者"] - res)**2
 errsum = sum(err)
 return errsum

# 穷举法,找出与实际数据差的平方和最小的S0和beta值
minSum = 1e10
minS0 = 0.0
minBeta = 0.0
bestRes = None

for S in S0:
 for b in beta:
  # 模型的差分方程
  def diff_eqs_2(INP, t):
   Y = np.zeros((3))
   V = INP
   Y[0] = -b * V[0] * V[1]
   Y[1] = b * V[0] * V[1] - gamma * V[1]
   Y[2] = gamma * V[1]
   return Y

  # 数值解模型方程
  INPUT = [S, I0, 0.0]
  RES = spi.odeint(diff_eqs_2, INPUT, t_range)
  errsum = error(RES[:21, 1])
  if errsum < minSum:
   minSum = errsum
   minS0 = S
   minBeta = b
   bestRes = RES
   print("S0=%d beta=%f minErr=%f" % (S, b, errsum))
print("S0 = %d β = %f" % (minS0, minBeta))
        
 
  • 1、
  • 2、
  • 3、
  • 4、
  • 5、
  • 6、
  • 7、
  • 8、
  • 9、
  • 10、
  • 11、
  • 12、
  • 13、
  • 14、
  • 15、
  • 16、
  • 17、
  • 18、
  • 19、
  • 20、
  • 21、
  • 22、
  • 23、
  • 24、
  • 25、
  • 1、
  • 2、
  • 3、
  • 4、
  • 5、
  • 6、
  • 7、
  • 8、
  • 9、
  • 10、
  • 11、
  • 12、
  • 13、
  • 14、
  • 15、
  • 16、
  • 17、
  • 18、
  • 19、
  • 20、
  • 21、
  • 22、
  • 23、
  • 24、
  • 25、