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

关键词

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

🏈指点迷津 | 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()

PreviousPython射频电磁肿瘤热疗数学模型和电磁爆炸性变化统计推理模型 NextPython流感常微分方程房室数学模型

Last updated 6 months ago

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