sir模型matlab案例_直播案例 | 人群接触网络中的 SIR 疫情模拟
![](/assets/blank.gif)
查看本案例完整的数据、代码和报告请登录数据酷客(http://cookdata.cn)案例板块。 快速获取案例链接、直播课件:数据科学人工智能公众号内发送“疫情案例”。
![](/assets/blank.gif)
https://www.zhihu.com/video/1239240729668898816
获取案例链接、直播课件、数据集在本公众号内发送“疫情数据”。
如何用网络来表示人之间的接触关系?在接触网络中,如何通过 SIR 模型模拟疫情的发展趋势?
本案例将介绍SIR模型,图和网络的基本知识。然后使用 networkx 工具,在生成的随机网络和真实的网络数据上,实现网络中的 SIR 模型进行疫情模拟。
1 SIR 模型介绍
SIR 模型用于计算封闭人群中随着时间推移感染传染性疾病的人数。最早提出来解释在瘟疫(1665-1666年伦敦,1906年孟买)和霍乱(1865年伦敦)等流行病中观察到的感染病人数量的迅速上升和下降。
![](/assets/blank.gif)
伦敦大瘟疫是1665年~1666年间,爆发在英国伦敦的大规模传染病,超过8万人死于瘟疫,相当于当时城市人口的1/5。
它假定人口规模是固定的(即无出生和自然死亡等)。传染源的潜伏期是瞬时的,传染持续时间与疾病的长度相同。人群是没有结构的,人之间的接触是完全随机的。恢复之后即获得免疫力,不会继续染病。
1.1 SIR 模型
SIR 是传染病中的最基础最核心的模型,可以使用下面的图形来表示。
![](/assets/blank.gif)
如图所示,在 SIR 模型中,人群被划分为以下三类:
- 易感者(Susceptible):未感染人群,与感染者接触后可能感染疾病,记为
或
;
- 感染者(Infectives):感染人群,感染者可能将疾病传染给易感者,也会以一定概率恢复健康,记为
或
;
- 恢复者(Recovered):感染后恢复健康的人,在可能致死的疾病中,也可能包括死亡人群,记为
或
。
SIR 模型包括两个重要的参数:感染率
在某个时间点,多少人会不幸染病呢?SIR 模型假设人之间的接触是随机的,人群中感染者数量是
SIR 模型中,感染者恢复的概率是
SIR 假设人群总数不变,即
可以看到,SIR 模型可以表示成一个常微分方程组:
1.2 使用 Scipy 求解 SIR 模型
上述常微分方程组要直接得到解析解(即把
SIR
来表示 SIR 模型。
def SIR(y,t,beta,gamma):S,I,R = ydSdt = - S*(I / (S + I + R))*betadIdt = beta*S*I/(S + I + R) - gamma*IdRdt = gamma*Ireturn [dSdt,dIdt,dRdt]
接下来我们就开始运用 scipy.integrate.odeint()
函数,获得微分方程组的解。
seed = 123 #本案例中将会使用的随机数种子
days = 100 #设置模拟的天数
beta = 0.30 #感染率
gamma = 0.10 #恢复率
N = 150 #人群大小
I0 = 1 #初始感染人数
R0 = 0 #初始恢复人数
S0 = N - I0 - R0# 设置初始值
y0 = [S0, I0, R0]
from scipy.integrate import odeint
# 求解
solution = odeint(SIR, y0, range(0,days), args = (beta, gamma))
将模拟的结果使用matplotlib
工具绘制出来,这里我们直接使用DataFrame对象封装的plot
方法。
import matplotlib.pyplot as plt
import pandas as pd
#plt.rcParams['figure.dpi'] = 100
solution_df = pd.DataFrame(solution, columns = ["S","I","R"])
#设置不同人群的显示颜色,易感者为橘色,感染者为红色,恢复者为绿色
color_dict = {"S":"orange","I":"red","R":"green"}
solution_df.plot(figsize=(9,6),color=[color_dict.get(x) for x in solution_df.columns])
![](/assets/blank.gif)
![](/assets/blank.gif)
为了更好的演示效果,我们使用pyecharts
将图动态地演示出来。
from pyecharts.charts import Line, Grid
import pyecharts
import pyecharts.options as opts
SIR_line = Line().add_xaxis(xaxis_data = solution_df.index)
for col in solution_df.columns:SIR_line.add_yaxis(series_name = col, y_axis = solution_df[col].values, # 添加上证指数数据symbol_size = 3,color = color_dict[col],label_opts = opts.LabelOpts(is_show = False),linestyle_opts = opts.LineStyleOpts(width = 1.5), # 设置线宽is_smooth = True)
# 在notebook中进行渲染
SIR_line.render_notebook()
![](/assets/blank.gif)
2 网络中的SIR模型
网络可以表示成图
在 SIR 模型中,假设人之间是随机接触的。如果人之间的接触关系不是随机的,而是形成了一个接触网络。那么在这个网络中,每个人接触到感染者的概率不再相等,而与他在网络中的位置相关。
![](/assets/blank.gif)
如上图的例子,人群中一共有8个人,其中
与传统 SIR 模型类似,有两个重要的参数:感染率
- 如果当前节点是恢复者,则下一步,节点状态依然是恢复者。
![](/assets/blank.gif)
- 如果当前节点是感染者,则下一步,
的概率转化为恢复者。
的概率依然是感染者。
![](/assets/blank.gif)
- 如果当前节点是易感者,则需要计算其邻居节点中感染者的数量,假设其有
个邻居为感染者。则该节点下一步转化为感染者的概率为
,否则继续保持易感者状态。
![](/assets/blank.gif)
对于本次新冠病毒肺炎疫情,如果要使用以上网络的方法对疫情进行模拟,我们还需要构建一个人之间的接触网络。从新闻的疫情报道以及本系列案例中第一次用爬虫抓取的数据,是完全不足以构建接触网络的。需要更多的数据,而这在当前阶段是十分困难的。
本案例中我们采用两种办法简单地构建一个网络结构:使用随机图生成算法生成一个无标度网络;使用一个真实的小型人群接触网络数据集。
3 生成无标度网络进行 SIR 疫情模拟
3.1 无标度网络
统计物理学家把服从幂律分布的现象称为无标度现象,即系统中个体的尺度相差悬殊,缺乏一个优选的标度。于是,满足幂律分布的网络也被称为无标度网络(scale-free network)。
无标度网络中,节点的度
其中
匈牙利科学家 Barabási 和 Albert 提出了一个 BA 模型(Barabási–Albert model)来生成无标度网络。 BA 模型整体流程如下:
- 首先,初始化
个相互连通的点;
- 然后,逐个地增加点,每增加一个点,这个新增的点与已有的点相连的概率正比于已有的点的度数。形式化地来说,对于点
而言,在新增的点与
之间增加一条边的概率
为
其中
3.2 使用 Networkx 生成无标度网络
Python 中的 Networkx 包提供了方便的随机网络生成函数。其中 barabasi_albert_graph
函数实现了 BA 模型生成无标度网络。主要的参数有网络节点数
draw_networkx
函数,它的常用参数包括网络对象、是否显示节点标签(with_labels)、网络的布局(pos)等。
#生成100个节点的BA无标度网络
import warnings
warnings.filterwarnings("ignore")import networkx as nx #导入networkx包,命名为nx
random_network = nx.barabasi_albert_graph(100,2) # 生成无标度网络,节点数和每个节点边数分别为100和2
#网络可视化
nx.draw_networkx(random_network,with_labels = True,pos = nx.spring_layout(random_network,seed = 1))
![](/assets/blank.gif)
3.3 SIR 疫情模拟与分析
节点状态模拟更新函数 updateNodeState
我们使用一个简单的函数来实现一个节点的状态的更新。首先,如果一个节点是恢复者,那么下一步还是恢复者,其节点状态保持不变。 如果一个节点是感染者,那么其恢复的概率是
![](/assets/blank.gif)
如一个节点是易感者,先要去其邻居节点中看看一共有多少个邻居是感染者,有
![](/assets/blank.gif)
updateNodeState
函数实现如下所示,其中输入的参数为网络
import random
# 根据 SIR 模型,更新节点的状态
def updateNodeState(G,node, beta, gamma):if G.nodes[node]["state"] == "I": #感染者p = random.random() # 生成一个0到1的随机数if p < gamma: # gamma的概率恢复G.nodes[node]["state"] = "R" #将节点状态设置成“R”elif G.nodes[node]["state"] == "S": #易感者p = random.random() # 生成一个0到1的随机数k = 0 # 计算邻居中的感染者数量for neibor in G.adj[node]: # 查看所有邻居状态,遍历邻居用 G.adj[node]if G.nodes[neibor]["state"] == "I": #如果这个邻居是感染者,则k加1k = k + 1if p < 1 - (1 - beta)**k: # 易感者被感染G.nodes[node]["state"] = "I"
网络中所有节点状态模拟更新函数 updateNetworkState
有了单个节点状态更新的函数,为了便于后续使用,我们再实现一个整个网络状态进行更新的函数。
def updateNetworkState(G, beta, gamma):for node in G: #遍历图中节点,每一个节点状态进行更新updateNodeState(G,node, beta, gamma)
为了动态地追踪易感者、感染者和恢复者数量,我们实现一个函数 countSIR
统计三类人群的数量。
# 计算三类人群的数量
def countSIR(G):S = 0;I = 0for node in G:if G.nodes[node]["state"] == "S":S = S + 1elif G.nodes[node]["state"] == "I":I = I + 1return S,I, len(G.nodes) - S - I
为了更清楚地演示在网络中疾病的传播过程,我们分别将图中的节点使用不同的颜色进行展示。易感者为橘色,感染者为红色,恢复者为绿色。
def get_node_color(G): #返回每一个节点的颜色组成的列表color_list = []for node in G:#使用我们前面创建的状态到颜色的映射字典 color_dict color_list.append(color_dict[G.nodes[node]["state"]])return color_list
在正式开始模拟之前,我们需要进行一些初始化工作,包括节点状态的初始化和 SIR 模型的参数
ba = nx.barabasi_albert_graph(N,2,seed=seed)
#初始化节点 state 属性
for node in ba:ba.nodes[node]["state"] = "S"
#随机选取一个节点为初始感染者
ba.nodes[55]["state"] = "I"
在图中开始 SIR 模型的模拟。设置模拟天数,开始执行模拟过程。
# 模拟天数为days,更新节点状态
import matplotlib.pyplot as plt
#fig,ax = plt.subplots(111)
%matplotlib inline
import time
SIR_list = []
for t in range(0,days):updateNetworkState(ba,beta,gamma) #对网络状态进行模拟更新SIR_list.append(list(countSIR(ba))) #计算更新后三种节点的数量
对模拟的结果进行可视化,查看易感者、感染者和恢复者人数的变化趋势。
df = pd.DataFrame(SIR_list,columns=["S","I","R"])
df.plot(figsize=(9,6),color=[color_dict.get(x) for x in df.columns])
![](/assets/blank.gif)
![](/assets/blank.gif)
我们再次画出在相同的
solution_df.plot(figsize=(9,6),color=[color_dict.get(x) for x in solution_df.columns])
![](/assets/blank.gif)
![](/assets/blank.gif)
对比两张图,可以看到,在网络中,疫情发展更为迅猛,最高峰的感染人数也远远高于传统 SIR 模型的模拟结果。为什么呢?作为一个开放性的问题,留给大家自己去想吧。
上面的疫情模拟中展示了每一天不同人群的变化,那么在网络中每一天到底是哪些人感染了?我们可以通过 networkx
提供的网络可视化工具深入地分析。通过 networkx.draw_networkx
函数可以方便地将图画出来。
nx.draw_networkx(ba,with_labels = True, node_color = get_node_color(ba), pos = nx.spring_layout(ba,seed = 1))
![](/assets/blank.gif)
默认的展示并不是很清晰,我们需要一些额外的设置,例如标签颜色,去掉图片边框,去掉坐标轴刻度,重新设置图的大小等。
fig, ax = plt.subplots(figsize=(9, 6)) #将图的大小设置为 9X6
pos = nx.spring_layout(ba,seed = 1) #设置网络布局,将 seed 固定为 1
ax.axis("off") #关闭坐标轴
plt.box(False) #不显示框
nx.draw(ba,with_labels = True,font_color="white",node_color = get_node_color(ba), edge_color = "#D8D8D8",pos = pos, ax=ax)
![](/assets/blank.gif)
将每一天的节点状态以动画演示。我们借助 matplotlib
包中的动画模块 animation
来实现这一效果。
#运行出结果大概需要90秒
#初始化节点 state 属性
for node in ba:ba.nodes[node]["state"] = "S"
ba.nodes[55]["state"] = "I"#绘制网络节点颜色
fig, ax = plt.subplots(figsize=(9, 6))def graph_draw(i,G,pos,ax,beta,gamma): ## 实现动画中每一帧的绘制函数,i为第几帧ax.axis("off")ax.set_title("day " + str(i) + " 黄色(易感者),红色(感染者),绿色(恢复者)")plt.box(False)if i == 0: #第一帧,直接绘制网络nx.draw(G,with_labels = True, font_color="white", node_color = get_node_color(G), edge_color = "#D8D8D8",pos = pos, ax=ax)else: # 其余帧,先更新网络状态,再绘制网络updateNetworkState(G, beta, gamma)nx.draw_networkx_nodes(G, with_labels = True, font_color="white", node_color = get_node_color(G),pos = pos, ax=ax) plt.close()#演示网络动态变化
import matplotlib.animation as animation
from IPython.display import HTMLanimator = animation.FuncAnimation(fig, graph_draw, frames= range(0,days),fargs=(ba,pos,ax,beta,gamma), interval=200)
HTML(animator.to_jshtml())
![](/assets/blank.gif)
3.4 初始感染者的影响
初始感染者在网络中位置不一样,对疾病的传播也会造成影响。这里我们用两种策略进行对比分析:
- 随机选择初始感染者;
- 选择度数最高的点作为初始感染者。
首先,我们随机选择一个种子节点。
import numpy as np
seed_node = np.random.randint(0,N)
print(seed_node)
#初始化节点 state 属性
for node in ba:ba.nodes[node]["state"] = "S"
ba.nodes[seed_node]["state"] = "I"# 模拟days天,更新节点状态
SIR_list = []
for t in range(0,days):updateNetworkState(ba,beta, gamma)SIR_list.append(list(countSIR(ba)))
df = pd.DataFrame(SIR_list,columns=["S","I","R"])
df.plot(figsize=(9,6),color=[color_dict.get(x) for x in df.columns])
![](/assets/blank.gif)
![](/assets/blank.gif)
然后,计算图中节点的度,选择度最高的节点作为种子感染者。
# 使用networkx 的 degree_centrality 函数计算图中节点的度中心度
node_degree = nx.degree_centrality(ba)node_degree_df = pd.DataFrame.from_dict(node_degree, orient = "index", columns=["degree"])
node_degree_df = node_degree_df.reset_index().rename(columns = {"index":"node"})#查看度数最高的节点
node_degree_df.sort_values(by = "degree",inplace = True,ascending= False)
node_degree_df.head()
![](/assets/blank.gif)
下面我们将度数最高的节点设置为种子感染者,观察疫情变化曲线。
seed_node = int(node_degree_df.values[0,0])
print(seed_node)
#初始化节点 state 属性
for node in ba:ba.nodes[node]["state"] = "S"
ba.nodes[seed_node]["state"] = "I"# 模拟days天,更新节点状态
SIR_list = []
for t in range(0,days):updateNetworkState(ba,beta, gamma)SIR_list.append(list(countSIR(ba)))
df = pd.DataFrame(SIR_list,columns=["S","I","R"])
df.plot(figsize=(9,6),color=[color_dict.get(x) for x in df.columns])
![](/assets/blank.gif)
![](/assets/blank.gif)
4 真实人群网络中的 SIR 疫情模拟
我们这里选用一个真实的人与人之间的接触网络来模拟疫情的传播。Infectious 是 2009 年在都柏林举办的主题为 “INFECTIOUS: STAY AWAY” 的展览活动中人与人之间面对面接触的网络。图中一共包含 410 个节点,每个节点表示一个人。如果两个人之间有超过20秒以上的面对面接触,则它们之间存在一条边。原始数据集中两个节点之间可能存在多条边,为了简化分析我们只保留其中的一条边。数据集来源于网站KNOECT。原始图数据见文件out.sociopatterns-infectious.csv,去除了多余边的网络文件见infectious.csv。
首先,我们使用 networkx
的 read_edgelist
函数加载网络,并简单可视化。
infectious_network = nx.read_edgelist("./input/infectious.csv",delimiter=",")
fig, ax = plt.subplots(figsize=(24, 16)) #节点较多,将图片大小也调整大些
pos_infectious = nx.spring_layout(infectious_network,seed = 22)
ax.axis("off")
plt.box(False)
nx.draw(infectious_network,with_labels = True,font_color="white", node_color = "orange", edge_color = "#D8D8D8",pos = pos_infectious, ax=ax)
![](/assets/blank.gif)
选择最上图中左边的 397 号节点作为种子节点,观看疫情的变化。注意,398 号节点处于网络中边缘的位置,在模拟过程中有可能尚未将疾病传播出去 397 号就恢复了健康,因此疾病不会在网络中继续传播。
beta = 0.10 # 为了更好观察,我们减小传染率参数
gamma = 0.05
#初始化节点 state 属性
for node in infectious_network:infectious_network.nodes[node]["state"] = "S"
infectious_network.nodes["397"]["state"] = "I"# 模拟days天,更新节点状态
SIR_list = []
for t in range(0,days):updateNetworkState(infectious_network, beta, gamma)SIR_list.append(list(countSIR(infectious_network)))df = pd.DataFrame(SIR_list,columns=["S","I","R"])df.plot(figsize=(9,6),color=[color_dict.get(x) for x in df.columns])
![](/assets/blank.gif)
![](/assets/blank.gif)
选取 272 号节点作为种子节点,观察疫情发展。
#初始化节点 state 属性
for node in infectious_network:infectious_network.nodes[node]["state"] = "S"
infectious_network.nodes["272"]["state"] = "I"# 模拟days天,更新节点状态
SIR_list = []
for t in range(0,days):updateNetworkState(infectious_network, beta, gamma)SIR_list.append(list(countSIR(infectious_network)))df = pd.DataFrame(SIR_list,columns=["S","I","R"])
df.plot(figsize=(9,6),color=[color_dict.get(x) for x in df.columns])
![](/assets/blank.gif)
![](/assets/blank.gif)
观察上述两个图中疫情高峰到来的时间,可以看到网络中重要的节点可以更快地传播疾病。这些人也是疫情防控中的关键。最后我们通过一个动画来感受下在这个真实的接触网络中疫情的传播态势。
#初始化节点 state 属性
#plt.rcParams['figure.dpi'] = 150
plt.rcParams['animation.embed_limit'] = 217277560for node in infectious_network:infectious_network.nodes[node]["state"] = "S"
infectious_network.nodes["397"]["state"] = "I"fig, ax = plt.subplots(figsize=(16, 12))
animator = animation.FuncAnimation(fig, graph_draw, frames= range(0,days),fargs=(infectious_network,pos_infectious,ax,beta,gamma), interval=200)
HTML(animator.to_jshtml())
![](/assets/blank.gif)
5 总结
在本案例中,我们首先介绍了传染病中核心的 SIR 模型。然后使用 Scipy
中的 odeint
函数对其进行数值求解,模拟疫情的传播。
在基本的 SIR 模型中假设人之间的接触是随机的。而在真实情况中,人与人的接触以网络形式存在。为了探索在网络中SIR模型的传播。我们介绍了一个网络中的 SIR 模型,借助 networkx
工具,我们实现了该模型。
进一步地,我们使用 networkx
提供的随机图生成算法利用 BA 模型生成了一个无标度网络,并在该网络中对疫情的传播进行了模拟,同时与基本的 SIR 模型进行了对比分析。
最后我们利用一个包含 410 个节点的真实接触网络上对疫情进行了模拟分析。
本文介绍的方法,在更多领域中有很多有趣的应用,例如社交网络中的广告营销,专家发现。感兴趣的可以关注社交网络分析或复杂网络领域中的影响最大化方面的研究。
往期直播
数据科学人工智能:使用Python爬取COVID-19疫情数据zhuanlan.zhihu.com
![](/assets/blank.gif)
数据科学人工智能:Pandas疫情探索性分析zhuanlan.zhihu.com
![](/assets/blank.gif)
数据科学人工智能:基于PyEcharts的COVID-19疫情可视化分析zhuanlan.zhihu.com
![](/assets/blank.gif)
数据科学人工智能:直播案例 | 使用 SIR 模型进行疫情模拟预测zhuanlan.zhihu.com
![](/assets/blank.gif)
sir模型matlab案例_直播案例 | 人群接触网络中的 SIR 疫情模拟相关推荐
- python 复杂网络中的 SIR 模型
本文地址:https://goodgoodstudy.blog.csdn.net/article/details/110152880 如下图所示的传染病传播过程,绿色为正常.红色为染病.蓝色为康复(并 ...
- 一对一视频直播源码实现网络中一对一视频聊天
一对一视频直播源码实现网络中一对一视频聊天 代码实现步骤 概述 首先要通信那就得满足通信的基础,我选择和目标放通信,前提就是我通过一定的条件将自己和目标建立链接,然后再将自己的通信信息交给目标,目标也 ...
- sir模型初始值_传播模型(SIR)
#include#include#include#include#include#include #define MaxVertexNum 90000 #define RAND_MAX 0x7fff ...
- ahp层次分析法matlab代码_(案例)AHP层次决策分析Matlab编码计算
"层次聚类分析Matlab编码计算" 01 - AHP层次决策分析计算函数 求判断矩阵最大特征根和归一化特征向量: function [maxEigVal,w] = maxEigV ...
- 用python对人口模型数据拟合_万字案例 | 用Python建立客户流失预测模型(含源数据+代码)...
1.采集数据 本数据集来自DF ,数据源地址: https://www.datafountain.cn/dataSets/35/details# 本数据集描述了电信用户是否流失以及其相关信息,共包含7 ...
- fvdm 跟驰模型 matlab仿真_强大的系统级热流体仿真软件Flownex了解一下,还有大咖免费培训哦...
点击上方蓝字,关注并设为星标 \ 知识分享 · 精品课程 · 工程仿真 \ 随着仿真在工业领域的逐渐深入,企业对仿真的需求也越来越高,不仅需要对零件.部件进行详细的仿真.设计和优化,也需要对系 ...
- 异步调用案例_异步案例研究
异步调用案例 In March 2020 our work related world changed drastically due to COVID-19. Previously our team ...
- 词云分析案例_品牌案例中的案例研究词
词云分析案例 A close-up look at Friday's design process. 近看星期五的设计过程. WORDS IN THE WILD IS A BAY-AREA NONPR ...
- 基于python的数据建模与分析案例_基于案例详解Python数据分析与机器学习
课程概述: 使用数据领域最主流语言Python及其分析与建模库作为核心武器.对于机器学习经典算法给出完整的原理推导并基于实例进行讲解,基于案例演示如何应用机器学习算法解决实际问题. 课程特色: 通俗易 ...
- html5 项目案例_互动案例技术分析(3)
这是该系列文章的第3篇,这次我们选择了三个稍微高阶一点的案例,使用了相对复杂的 Canvas 来实现.我们的目标并非是推广技术,而是展示技术所能实现的效果. 技术不是互动营销的全部,但技术可以让互动营 ...
最新文章
- Alpha冲刺 - 事后诸葛亮
- acwing算法题--01背包问题
- fullcalendar5.X版本 显示自定义html内容
- ajax实现highchart与数据库数据结合完整案例分析(三)---柱状折线图
- npm WARN build `npm build` called with no arguments. Did you mean to `npm run-script build`
- storm消息可靠机制(ack)的原理和使用
- i5 1135g7什么水平_i7-10510U和i5-1135G7对比,该怎么选择呢?
- angular_directive动感超人
- Unique Binary Search Trees ll -深度优先遍历DFS
- 异常错误 - MySQL导入时错误
- The VMRC console has disconnected solution
- C语言动态规划——背包问题详解
- 此beta版已额满_《魔域口袋版》福利狂欢:现金红包天天送 魔石神器免费拿
- 计算机图形学教程动画实验报告,计算机图形学画圆实验报告.doc
- 无码科技发布第一款产品:Readhub
- 旧手机怎么当文件服务器,用旧手机做云存储服务器
- 单个正态总体均值的区间估计_总体均值的区间估计 (正态总体: σ2 已知实例).pdf...
- 英语作业介绍一项发明计算机,计算机专业英语第1次作业.doc
- sdfasfasdf
- 短时交通预测方法总结