SIR模型及python复现

SIR模型

SIR模型是传染病模型中的经典模型，可以用在传染病过程中的模拟预测，也可以用作抽象表达社交网络中的信息传播过程，本文讲述的是后者。
SIR模型将社交网络中的节点分为三类：$S$类，易感染者，指未接收信息且不具备传播能力的节点；$I$类，感染者，指接收信息且具备传播能力的节点；$R$类，免疫者，指接收信息后退出传播过程的节点。
SIR模型的传播机制如下图所示：

其中，$\beta$代表感染率，即时刻 $t$ 时单位时间内$I$类节点新增的感染个体与$S$类节点个数成比例，比例系数为$\beta$；$\gamma$代表免疫率(恢复率)，即时刻 $t$ 时单位时间内$R$类节点新增的免疫个体与$R$类节点个数成比例，比例系数为$\gamma$。

SIR数学模型

传播动力学方程

我们将上述传播模型表达为动力学方程。这里设置节点总数为$N$，其他三类节点初始数目分别为$S、I、R$且$N=S+I+R$。则有下列方程
$$\{\begin{array}{l}\frac{dS}{dt}=-\beta S\frac{I}{N},\\ \frac{dI}{dt}=\beta S\frac{I}{N}-\gamma I,\\ \frac{dR}{dt}=\gamma I\end{array}$$
其中，$\frac{I}{N}$代表$S$类与$I$类的接触率。

python实现

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as pltdef sir_model(input_v, t):s, i, r = input_vn = s + i + rds_dt = - beta * s * (i / n)di_dt = beta * s * (i / n) - gamma * idr_dt = gamma * ireturn [ds_dt, di_dt, dr_dt]if __name__ == '__main__':# 模型输入值N = 10000  # 设定总节点数量S0 = 9980  # S类节点数量I0 = 18  # I类节点数量R0 = N - S0 - I0  # R类节点数量input_value = (S0, I0, R0)beta = 0.5  # 设定感染率gamma = 0.1  # 设定免疫率# 设定步长 即传播次数times = np.linspace(0, 60, 60)# 求解微分方程result = odeint(sir_model, input_value, t=times)# 画图plt.plot(result[:, 0], '-ro', label='S')plt.plot(result[:, 1], '-b^', label='I')plt.plot(result[:, 2], '-gs', label='R')plt.legend(loc=0)plt.xlabel('times')plt.ylabel('number')plt.show()

实现结果如下图所示。

模拟社交网络中SIR模型的信息传播过程

在这里由于没有数据集，所以只能用python中的Networkx库，这里我们使用BA无标度网络模型(个人认为比较接近于现实社交网络)。

import networkx as nx
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import random
# 确保 中文 和 -
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = Falsedef update_node_status(G, node, beta, gamma):"""更改节点状态:param G: 输入图:param node: 节点序数:param beta: 感染率:param gamma: 免疫率"""# 如果当前节点状态为 感染者(I) 有概率gamma变为 免疫者(R)if G.nodes[node]['status'] == 'I':p = random.random()if p < gamma:G.nodes[node]['status'] = 'R'# 如果当前节点状态为 易感染者(S) 有概率beta变为 感染者(I)if G.nodes[node]['status'] == 'S':for adj_node in G[node]:if G.nodes[adj_node]['status'] == 'I':p = random.random()if p < beta:G.nodes[node]['status'] = 'I'breakdef update_network_data(G, beta, gamma):"""更改图数据:param G: 输入图:param beta: 感染率:param gamma: 免疫率"""for node in G:update_node_status(G, node, beta, gamma)def initial_network_data(G, i_num, r_num):"""初始化图数据:param G: 输入图:param i_num: 感染者数量:param r_num: 免疫者数量"""# 感染节点集i_set = set(random.sample(G.nodes, i_num))# 免疫节点集r_set = set(random.sample(G.nodes, r_num))# 两个集合不能重复while r_set & i_set:r_set = set(random.sample(G.nodes, r_num))# 初始化节点状态for node in G:if node in i_set:G.nodes[node]['status'] = 'I'elif node in r_set:G.nodes[node]['status'] = 'R'else:G.nodes[node]['status'] = 'S'def count_node(G):"""计算当前图内各个节点的数目:param G: 输入图:return: 各个节点数目"""s_num, i_num, r_num = 0, 0, 0for node in G:if G.nodes[node]['status'] == 'S':s_num += 1elif G.nodes[node]['status'] == 'I':i_num += 1else:r_num += 1return s_num, i_num, r_numdef draw_network(G):"""输出初始网络节点分布:param G: 输入图"""# 设置图大小fig, ax = plt.subplots(figsize=(12, 8))ax.set_title("易感染者-感染者-免疫者节点初始分布")pos = nx.spring_layout(G, scale=1)nx.draw(G, pos=pos, with_labels=True, font_color='white', edge_color='grey',node_color=[color_dict[ba.nodes[node]['status']] for node in G])def draw_node_trend(G, beta, gamma):"""输出各类节点趋势:param G: 输入图:param beta: 感染率:param gamma: 免疫率"""# 设定传播步长t_list = np.linspace(1, 24, 24)# 开始模拟传播for t in range(len(t_list)):# 计算并存储当前各个节点数目node_list.append(count_node(G))update_network_data(G, beta, gamma)# 整理数据df = pd.DataFrame(data=node_list, index=t_list, columns=['S', 'I', 'R'])# 显示数据曲线df.plot(figsize=(8, 6), color=[color_dict.get(x) for x in df.columns])plt.ylabel('nodes(节点数)')plt.xlabel('days(天数)')plt.title('易感染者-感染者-免疫者节点趋势')plt.savefig('SIR_model')plt.show()def update_graph(i, G, ax, pos, beta, gamma):"""动态更新节点:param i: 输入帧:param ax: 输入图参数:param G: 输入图:param beta: 感染率:param gamma: 免疫率"""i = int(i)ax.set_title("第" + str(i + 1) + "天 易感染者-感染者-免疫者节点分布")ax.axis('off')plt.box(False)if i == 1:# 第一天  初始节点分布  直接画出nx.draw(G, with_labels=True, font_color='white', edge_color='grey',node_color=[color_dict[ba.nodes[node]['status']] for node in G], pos=pos)else:# 以后变化 需要演变节点update_network_data(G, beta, gamma)nx.draw_networkx_nodes(G, with_labels=True, font_color='white', edge_color='grey',node_color=[color_dict[ba.nodes[node]['status']] for node in G], pos=pos)def draw_network_trend(G, beta, gamma, days):"""输出网络动态变化视频:param G: 输入图:param beta: 感染率:param gamma: 免疫率:param days: 需要的时间(迭代次数)"""fig, ax = plt.subplots(figsize=(12, 8))pos = nx.spring_layout(G, scale=1)ani = animation.FuncAnimation(fig, update_graph, frames=days,fargs=(G, ax, pos, beta, gamma), interval=300)writer = animation.FFMpegWriter()ani.save('network_trend.mp4', writer=writer)if __name__ == '__main__':# 总人数N = 10000# 易感染者人数s = 9980# 感染者人数i = 18# 免疫者人数r = N - s - i# 各个节点数目列表node_list = []# 节点颜色color_dict = {'S': 'blue', 'I': 'red', 'R': 'green'}# 创建BA无标度网络ba = nx.barabasi_albert_graph(N, 3, seed=1)# 初始化网络节点initial_network_data(ba, i, r)# 输出节点趋势图draw_node_trend(ba, 0.4, 0.1)# 初始节点分布图# 输出初始节点网络图# draw_network(ba)# 输出网络动态变化图# draw_network_trend(ba, 0.2, 0.15, 50)

上述代码得到的结果如下图所示。

这里说一下，如果有数据集可以修改一下这个代码，将数据集中的边集和节点集输入Networkx中的空白图中，剩余可以不改变，也可以模拟SIR模型在真实社交网络下的传播过程。
从SIR模型中我们可以看到，$\beta$和$\gamma$在很大程度会影响信息传播过程。如果利用真实数据集，也要注意不同的$\beta$和$\gamma$值，会有不同的结果。同时，利用好真实数据集，我们可以求出近似真实值$\beta$和$\gamma$值。
另外在上述主程序代码的最后几行，我增加了生成网络初始分布图和动态实现社交网络中节点变化的函数。动态实现那部分需要一个ffmpeg的软件才可以生成。