🍠Python流感传播感染康复图模型计算和算法

关键词

Python | 流感 | 传播 | 疑似 | 感染 | 康复 | 图模型 | 算法 | 数学 | 网络 | 规模分布 | 随机 | 动态 | 预期 | 节点 | 同质 | 异质 | 预期 | 双峰 | 幂律 | 建模 | 方程 | 概率 | 估计 | 泊松 | 贝叶斯 | 动画 | 非马尔可夫 | 分层

Python流感传播感染康复图模型计算和算法 | Notion亚图跨际 on Notion

🏈指点迷津 | Brief

🎯要点

🎯流感传播网络：🖊截断幂律网络、埃尔多斯-雷尼网络、同构网络 | 🖊绘制最终流感规模分布 | 🖊集总系统和随机模拟流感动态特征 | 🖊绘制模型的预期流行率。

🎯受感染节点数量时间依赖性同质模型：🖊埃尔多斯-雷尼随机图 | 🖊精确系统预期受感染节点数量| 🖊常规随机网络受感染节点预期数量 | 🖊随机模拟受感染节点数量 | 🖊随机模拟埃尔多斯-雷尼受感染节点预期数量 | 🖊随机模拟双峰受感染节点预期数量 | 🖊随机模拟规则、双峰、埃尔多斯-雷尼和幂律预期受感染节点数量。

🎯受感染节点数量时间依赖性异质模型：🖊随机模拟、紧凑成对模型和齐次成对模型 | 🖊埃尔多斯-雷尼随机网络 | 🖊幂律配置模型 | 🖊双峰随机网络。

🎯渗透疾病建模 | 🎯疑似感染康复模型分层 | 🎯非马尔可夫流感 | 🎯动画感染力传播。

🎯 流感模型：Python流感常微分方程房室数学模型 | 🎯图模型：Python图嵌入信息潜在表征算法、Python非线性图嵌入和降维技术拉普拉斯特征图算法。

🍇Python传染病传播概率估计

模型用常微分方程给出了不同人群的动态，假设整个人口是一个常数N，忽略流行期间的出生率和死亡率。即：

$$\left\{\begin{array}{l} \frac{d S}{d t}=-\frac{\beta I S}{N}, \\ \frac{d I}{d t}=\frac{\beta I S}{N}-\gamma I, \\ \frac{d R}{d t}=\gamma I, \end{array}\right.$$

β和γ这两个参数分别代表感染率和清除率。 它们的比率是基本再生数，其值决定了传染病的渐近行为。 一旦给出这两个参数，模型所描述的动态就可以确定。

依据常微分方程给出，传染病传播模型在放入网络之前应进行离散化。 在这里，我们将使用经典的四阶龙格-库塔积分，这是一种采用逐次逼近求解形式为 dy/dx=f(x,y) 的微分方程的数值方法。

$$\begin{aligned} & K_1=h f\left(x_n, y_n\right) \\ & K_2=h f\left(x_n+\frac{h}{2}, y_n+\frac{k_1}{2}\right) \\ & K_3=h f\left(x_n+\frac{h}{2}, y_n+\frac{k_2}{2}\right) \\ & K_4=h f\left(x_n+h, y_n+k_3\right) \\ & y_{n+1}=y_n+k_1 / 6+k_2 / 3+k_3 / 3+k_4 / 6+O\left(h^5\right) \end{aligned}$$

• 将 beta 和 gamma 的先验设置为 [0,1] 内的均匀分布；

• 给定模型的初始状态（用 y_0​ 表示），使用龙格-库塔方法获取未来实例的状态序列。

• 为每个实例添加逐项泊松噪声。

• 进行贝叶斯推理以获得 beta 和 gamma 的后验。

在 100 个数据点上进行训练后，参数收敛到以 0.6 和 0.3 为中心的正态分布。

计算并绘制图形

 import numpy as np
 ​
 def rk_4(y_0,t,f,args=()):
   n = len(t)
   y = np.zeros((n, len(y0)))
   y[0] = y0
   for i in range(n - 1):
       h = t[i+1] - t[i]
       k1 = f(y[i], t[i], *args)
       k2 = f(y[i] + k1 * h / 2., t[i] + h / 2., *args)
       k3 = f(y[i] + k2 * h / 2., t[i] + h / 2., *args)
       k4 = f(y[i] + k3 * h, t[i] + h, *args)
       y[i+1] = y[i] + (h / 6.) * (k1 + 2*k2 + 2*k3 + k4)
   return y, k2
 ​
 ​
 N=1e4
 y0=np.array([N-100, 100, 0])
 beta=0.6 
 gamma=0.3  
 T_max= 100 
 dt=1
 T_num=int(T_max/dt)
 t= np.linspace(start=0,stop=T_max,num=int(T_num))
 ​
 ​
 def f_SIR(y,t,N,beta, gamma):
   f_SIR=np.zeros(3)
   f_SIR[0]= - beta/N*y[0]*y[1]
   f_SIR[1]= beta/N*y[0]*y[1]-gamma*y[1]
   f_SIR[2]= gamma*y[1]
   return f_SIR
 ​
 y=rk_4(y0, t, f_SIR,args=(N, beta, gamma))
 import matplotlib.pyplot as plt
 for i in range(3):
   plt.plot(t, y[:,i])
 plt.title("SIR model")

读取数据：

 import pandas as pd
 cases=pd.DataFrame(
     {
         "reported_infection": infections,
      "recovered": recovered,
      "population":N,
 ​
     }
 )
 dates=pd.date_range(start="2022-01-01", periods=len(cases), freq='D')
 cases["date"]=dates
 cases.to_csv("daily_cases.csv")

现在考虑使用贝叶斯模型来预测。

 y_real=y.astype(int)
 ​
 np.random.seed(0)
 y_observed=[]
 t_observed=[]
 for i in range(len(y_real)):
   if i % (len(y_real)/100)==0:
     y_observed.append(y_real[i])
     t_observed.append(t[i])
 y_observed=np.vstack(y_observed)
 y_observed=np.random.poisson(y_observed)

 import matplotlib.pyplot as plt
 plt.scatter(t_observed, y_observed[:,0],color=(0,1,0),label="S",s=1)
 plt.scatter(t_observed, y_observed[:,1],color=(1,0,0),label="I",s=1)
 plt.scatter(t_observed, y_observed[:,2],color=(0,0,1),label="R",s=1)
 ​
 plt.legend()