SIR模型的应用 - Influence maximization in social networks based on TOPSIS(3)

在Influence maximization in social networks based on TOPSIS一文中，作者利用SIR模型如下：

即：种子节点处于感染阶段，其余节点处于易感，最终的评判以网络中恢复的人数来判断。
那么，不妨来了解下SIR模型是如何工作的。

`SIR`模型

经典的SIR模型提出比较早，我们都知道它存在三种状态，分别是：
Susceptible：易感人群，指未得病者，但缺乏免疫能力，与感病者接触后容易受到感染；
Infective：感染人群，指染上传染病的人，他可以传播给易感人群；
Removed：移除人群，指因病愈（具有免疫力）或死亡的人。这部分人不再参与感染和被感染过程。

随便百度下，我们都可以看见一些比较复杂的微分公式，如下图：

那么，接下来就来看看它是如何计算出来的。根据百度百科中介绍SIR模型的建立有几个基本假设：

基本假设

① 总人数N不变，即：不考虑人口的出生、死亡、流动等因素；即：
S+I+R=NS + I + R = N S+I+R=N
② 易感者、感染者和移出者的比例分别是s(t)、i(t)和r(t)，满足：
s(t)+i(t)+r(t)=1s(t) + i(t) + r(t) = 1 s(t)+i(t)+r(t)=1
③ 易感者数目S和感染者数目I成正比，比例系数为β(感染系数)；那么，单位时间(Δ(t)\Delta(t)Δ(t))内被感染的人数就是：
βNs(t)i(t)Δ(t)βNs(t)i(t)\Delta(t) βNs(t)i(t)Δ(t)
可以理解为：Ni(t)Ni(t)Ni(t)表示感染者的总人数，以感染系数βββ进行传播，而这种感染效应作用在以s(t)s(t)s(t)为比例的人群，而传播的时间为Δ(t)\Delta(t)Δ(t)。
也就是可以理解为，易感人群在时间Δ(t)\Delta(t)Δ(t)内，对易感人群进行的一个感染的结果。
④ 单位时间内从染病者中移出的人数与病人数量成正比，比例系数为γ(恢复系数)，单位时间(Δ(t)\Delta(t)Δ(t))内移出者的人数为：
Ni(t)γΔ(t)Ni(t)γ\Delta(t) Ni(t)γΔ(t)
和上面的区别在于没有再乘r(t)r(t)r(t)的比例，因为Δr\Delta rΔr和RRR实际是无关的（即：现在治愈/死亡多少人和将来能治愈/死亡多少人无直接关系）；而Δi\Delta iΔi的产生却和III是有关的（即：感染的人数一定是和原先就存在的感染人数是相关的）。

我们可以得出下面的一些结论：

相关结论

① 感染个体的增长率=期望的增长率-死亡/恢复率：
βs(t)i(t)−γi(t)βs(t)i(t) - γi(t) βs(t)i(t)−γi(t)
② 易感个体的下降率：
βs(t)i(t)βs(t)i(t) βs(t)i(t)
③ 移除个体的增长率：
γi(t)γi(t) γi(t)

易感者从患病到移出的过程可以用微分方程表示如下：
{ds(t)dt=−βs(t)i(t);di(t)dt=βs(t)i(t)−γi(t);di(t)dt=γi(t);i(0)=i0;s(0)=s0;通常初始取r(0)=0，i(0)+s(0)+r(0)=1;\begin{cases} \frac{ds(t)}{dt} = - βs(t)i(t) ;\\ \frac{di(t)}{dt} = βs(t)i(t) - γi(t) ;\\ \frac{di(t)}{dt} = γi(t) ;\\ i(0) = i_0 ;\\ s(0) = s_0 ;\\ 通常初始取r(0)=0，i(0) + s(0) + r(0) = 1; \end{cases}⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎧dtds(t)=−βs(t)i(t);dtdi(t)=βs(t)i(t)−γi(t);dtdi(t)=γi(t);i(0)=i0;s(0)=s0;通常初始取r(0)=0，i(0)+s(0)+r(0)=1;
关于i(t) , s(t) 的非线性微分方程组，没有解析解，只能通过数值计算得到s(t), i(t), r(t)的曲线。

那么，我们不妨用Python来模拟下：

# %matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeintclass SIR(object):def __init__(self, beta, gama, s0, i0, r0):"""初始化模型参数:param beta: 感染系数:param gama: 恢复系数"""self.beta = betaself.gama = gamaself.initalization(s0, i0, r0)def initalization(self, s0, i0, r0):"""初始化易感者、感染者和移除者在总体人数N中的比例:param s0::param i0::param r0::return:"""self.s0 = s0self.i0 = i0self.r0 = r0def sir_model(self, y, t, beta, gama):"""定义模型函数，以输入到odeint()函数odeint()函数是scipy库中一个数值求解微分方程的函数即：dy/dt = [ds_dt, di_dt, dr_dt]  的微分方程"""s, i, r = yds_dt = - beta * s * idi_dt = beta * s * i - gama * idr_dt = gama * ireturn np.array([ds_dt, di_dt, dr_dt])def calc(self):# 均匀生成1-100的100000个数的时间点time = np.linspace(1, 100, 100000)# odeint()函数是scipy库中一个数值求解微分方程的函数return time, odeint(self.sir_model, [self.s0, self.i0, self.r0], time, args=(self.beta, self.gama))def plot(self):time, res = self.calc()plt.figure()plt.plot(time, res[:, 0], label="S(t)")plt.plot(time, res[:, 1], label="I(t)")plt.plot(time, res[:, 2], label="R(t)")plt.legend()plt.grid()plt.xlabel("time")plt.ylabel("proportions")plt.title("SIR model simulation")plt.show()if __name__ == '__main__':SIR(0.55, 0.1, 0.8, 0.2, 0).plot()

在上面的过程中，虽然拟合了SIR模型，但是存在一个问题，就是文中是如何处理的，也就是如何处理：所选定的不同种子节点集合在SIR传播模型上，最终传播效果有所不同？

因为，在上面SIR模型中，并非是针对特定的节点进行模拟，而是根据初始S/I/R三种状态的比例，以及感染系数beta、恢复系数gama来进行，而没有考虑到不同的节点的因素。那么如何确定所选定的K个种子节点，就是选择的最优的呢？

这个问题我们在下篇进行思考。