目前整个世界同新冠肺炎的斗争仍在继续，从新冠肺炎爆发以来，很多人都提出了各种各样的模型来预测肺炎疫情的发展，其中比较常见的就是SIR模型。这是因为SIR是目前疾病防控领域最经典也是最常用的一个模型，而今天笔者就用图论来讲述一下SIR模型。
首先对SIR模型做一个简单的介绍，SIR全称就是Susceptible—Infected—Recovered，翻译过来就是易感—感染—康复，即易感人群Susceptible有α概率被某种疾病感染，成为感染人群Infected，而感染人群Infected又有β概率康复，成为康复人群Recovered。这个过程可以是一次或者是多次的，而且还衍生出不少类似的模型，比如SIS、SIRS和SI等。整个SIR过程可以用下图来描述。
图1. SIR模型示意图
而笔者这次用的是图论，也就是我们平常所说的网络分析来说明SIR模型，大概过程如下图所示。
图2. SIR网络模型示意图
这个过程其实很简单，就是由少数人感染多数人，然后最终所有感染人群都康复的一个过程，这个模型简化了我们很多实际情况所遇到的部分问题，比如不考虑感染人群的死亡问题，不考虑感染人群的再感染等等。这次用到的图论方面的库是networkx，这是Python中最常用的图论分析的库。下面还是直接用代码来具体说明一下。
首先导入各种包。
<code>import matplotlib.pyplot as plt
import numpy.random as rdm
import networkx as nx/<code>
接下来定义两个变量n和g。n是总人数，而nx.erdos_renyi_graph则是一种随机网络图，Erdos和Renyi是两位匈牙利数学家，他们的主要贡献就是建立了著名的E-R随机图理论（Random graph theory），这被公认为在数学上开创了复杂网络拓扑结构的系统性分析。而这个E-R图也是本文应用的基础。而nx.erdos_renyi_graph(n, 0.01)的意思就是以概率0.01来连接100个节点。这里100和0.01这两个数字笔者用的比较随意，大家可以根据自己的研究来任意设定相关数值，这里主要还是用于简单说明。然后是三个string类型的变量susceptible、infected和recovered，分别代表“易感”、“感染”和“康复”人群。
<code>n = 100
g = nx.erdos_renyi_graph(n, 0.01)
susceptible = 'S'
infected = 'I'
recovered = 'R'/<code>
接下来我们要定义几个函数。
<code>def onset(g): #初始设置
for i in g.node.keys():
g.node[i]['state'] = susceptible
def infect_prop(g, proportion): #设置感染比例
for i in g.node.keys():
if(rdm.random() <= proportion):
g.node[i]['state'] = infected
def build_model(pInfect, pRecover): #模型构建
def model(g, i):
if g.node[i]['state'] == infected:
for m in g.neighbors(i):
if g.node[m]['state'] == susceptible:
if rdm.random() <= pInfect:
g.node[m]['state'] = infected
if rdm.random() <= pRecover:
g.node[i]['state'] = recovered
return model
def model_run(g, model): #单次模型运行
for i in g.node.keys():
model(g, i)
def model_iter(g, model, iter): #多次模型循环
for i in range(iter):
model_run(g, model)/<code>
第一个函数onset，也就是刚开始的状态，在这里我们把每个节点的状态都设置为“S”，也就是易感的意思，这里我们假设所有人都属于易感人群中的。然后就开始有部分人群出现感染，而这个感染函数就是infect_prop，prop意思就是proportion，这个函数的意思就是感染比例，这里我们用到了numpy.random.random方法，也就是返回一个均匀分布，数值大小在0-1之间，不包含1，当rdm.random()<= proportion时，我们就让这样的人群变为感染人群，这样让感染人群更加均匀一些。第三个函数是build_model，是一个嵌套函数，有两个参数，pInfect和pRecover，分别代表感染概率和康复概率，而这里的g.neighbors(i)的意思是节点i的每个相邻的节点，而g是我们生成的一个图（graph），是一个class，后面会有说明。在这些感染人群中，当他们的相邻节点的状态是“susceptible”时，我们让这些节点（人群）以概率pInfect来进行感染；而在这些感染人群中，以pRecover的概率进行康复。这就是一个简单的感染—康复过程。而后面的两个函数model_run和model_iter则是将这个模型运行一次和多次，分别用来模拟一个循环和多个循环。
接下来就是画图。
<code>fig, ax= plt.subplots(figsize=(12, 10))
ax.set_xticks([])
ax.set_yticks([])
pos = nx.spring_layout(g, k=0.2)
nx.draw_networkx_edges(g, pos, alpha=0.5, width = 1)
nx.draw_networkx_nodes(g, pos, node_size=80)
plt.show()/<code>
首先设置图片的大小，并去掉坐标轴，然后设置网络图的位置，nx.spring_layout就是设置网络图的位置的方法，k是节点间的最佳距离，这个可以随意设置，值越大节点越分散，接下来绘制节点和连线，nx.draw_networkx_nodes(g, pos, node_size=80)用来绘制节点，pos就是刚才设置的位置参数，再设置一下节点的大小，而nx.draw_networkx_edges(g, pos, alpha=0.5, width = 1)用来绘制连线，同样要传入位置参数，再设置透明度和线宽。最后生成的图如下。
图3. 模型生成的网络图
networkx使用的绘图算法是随机的，同时我们使用的参数也是随机的，所以这个图每次生成的结果都不同，但大体相似，我们可以看到这里面已经有部分节点相连，疾病也就是通过他们开始传播。
最后就是计算感染率。
<code>onset(g)
infect_prop(g, 0.05)
model = build_model(0.2, 0.8)
model_iter(g, model, 10)
infected = [ v for (v, attr) in g.nodes(data = True) if attr['state'] == recovered ]
infection_rate = len(infected)/n
print(infection_rate)/<code>
这当中onset(g)中的g就是前面我们用g = nx.erdos_renyi_graph(n, 0.01)生成的图，这是一个类的实例，infect_prop(g, 0.05)中的0.05就是设置人群初始感染率为0.05，model = build_model(0.2, 0.8)中的感染几率和康复几率分别设置为0.2和0.8，至于为何这么设置，主要是根据常用的“二八理论”，这个数字可以根据使用者的模拟情况随意设置，而model_iter(g, model, 10)中我们让这个模拟过程重复10次，最后计算出感染人数，并得出最终稳定的感染率infection_rate，因为这些参数不少都是随机的，所以得出的结果理论上每次都是不同的，笔者得出的结果从0.03到0.12不等，这个结果意义不大，主要还是理解模拟的方法。
閱讀更多 愛學習的美少女 的文章
