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

🎯要点

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

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

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

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

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

🍇Python传染病传播概率估计

模型用常微分方程给出了不同人群的动态，假设整个人口是一个常数N，忽略流行期间的出生率和死亡率。即：
$$\{\begin{array}{l}\frac{dS}{dt}=-\frac{\beta IS}{N},\\ \frac{dI}{dt}=\frac{\beta IS}{N}-\gamma I,\\ \frac{dR}{dt}=\gamma I,\end{array}$$
β和γ这两个参数分别代表感染率和清除率。 它们的比率是基本再生数，其值决定了传染病的渐近行为。 一旦给出这两个参数，模型所描述的动态就可以确定。

依据常微分方程给出，传染病传播模型在放入网络之前应进行离散化。 在这里，我们将使用经典的四阶龙格-库塔积分，这是一种采用逐次逼近求解形式为 $dy/dx=f(x,y)$ 的微分方程的数值方法。
$$\begin{array}{rl}& K_1=hf\left(x_n,y_n\right)\\ & K_2=hf\left(x_n+\frac{h}{2},y_n+\frac{k_1}{2}\right)\\ & K_3=hf\left(x_n+\frac{h}{2},y_n+\frac{k_2}{2}\right)\\ & K_4=hf\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{array}$$

• 将 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()

🔗参阅一：计算思维

🔗参阅二：亚图跨际