信息传播学习笔记（2）——SIS模型python代码实现（基于全局与基于节点）_sis 模型代码

本文是讲解python代码实现SIS模型，共有两部分代码：
第一部分是传统的SIS模型，是获得全局感染者（I）和易感者（S）的比率。
第二部分是基于节点的SIS模型，是获得每个节点的受感染率。

对了宝贝儿们，卑微小李的公众号【野指针小李】已开通，期待与你一起探讨学术哟~摸摸大！

1 传统SIS模型

传统的SIS模型如果不熟悉的同学，可以参考我的上一篇文章《信息传播学习笔记（1）——SIS模型原理与公式推导》。

1.1 初始值设定

本文设定的初始值如下：

nodes = np.arange(0, 5, 1)
lamb = 0.5  # 传染率
mu = 0.2  # 治愈率
t_range = np.arange(0, 50, 0.1)  # 时间步长
init_i = 2  # 初始感染者
1
2
3
4
5

即$\lambda =0.5,\mu =0.2$，迭代次数500次，共5个节点，2个感染者（I），3个易感者（S）。

1.2 微分方程

def sis_i(i, t):
    """
    SIS模型
    di(t)/dt = lambda * i(t) * (1 - i(t)) - mu * i(t)

    :param i: 感染者
    :param t: 时间
    :return i: t时刻的i
    """
    i = lamb * i * (1 - i) - mu * i

    return i


def sis_s(s, t):
    """
    SIS模型
    ds(t)/dt = -lambda * s(t) * (1 - s(t)) - mu * (1 - s(t))

    :param s: 易感者
    :param t: 时间
    :return s: t时刻的s
    """
    s = lamb * s * (1 - s) - mu * (1 - s)

    return s
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

sis_i是计算感染者比率的微分方程，其公式为：
$$N\frac{di(t)}{dt}=Ni(t)\lambda s(t)-\mu Ni(t)$$

sis_s是计算易感者比率的微分方程，其公式为：
$$N\frac{ds(t)}{dt}=-Ni(t)\lambda s(t)+\mu Ni(t)$$

由于$s(t)+i(t)=1$，所以上面的公式转换为：
$$\frac{di(t)}{dt}=\lambda i(t)(1-i(t))-\mu i(t)$$

$$\frac{ds(t)}{dt}=-\lambda s(t)(1-s(t))+\mu s(t)$$

1.2 微分方程调用与绘图

通过scipy计算微分方程：

import scipy.integrate as spi

def calc_ODE(nodes):
    """
    计算常微分方程

    :param nodes: 节点
    """
    i = calc_init_ratio(nodes)
    result_i = spi.odeint(func=sis_i, y0=i, t=t_range)
    result_s = spi.odeint(func=sis_s, y0=i, t=t_range)
    draw_pic(result_i, result_s)
1
2
3
4
5
6
7
8
9
10
11
12

其中calc_init_ratio(nodes)是计算$t(0)$时刻感染率：

def calc_init_ratio(nodes):
    """
    计算初始时刻的i(0)

    :param nodes: 节点
    :return: 第0时刻的感染率
    """
    n = len(nodes)
    return init_i / n
1
2
3
4
5
6
7
8
9

其中通过draw_pic来绘图：

def draw_pic(result_i, result_s):
    """
    绘制感染情况图
    :param result_i: 感染者的比率构成的向量
    :param result_s: 易感者的比率构成的向量
    """
    plt.plot(t_range, result_i)
    plt.plot(t_range, result_s)
    plt.xlabel("time")
    plt.ylabel("ratio")
    plt.show()
1
2
3
4
5
6
7
8
9
10
11

1.3 结果展示

其中蓝色的线是感染者I，橘色的线是易感者S。

1.4 完整代码

import scipy.integrate as spi
import matplotlib.pyplot as plt
import numpy as np

nodes = np.arange(0, 5, 1)
lamb = 0.5  # 传染率
mu = 0.2  # 治愈率
t_range = np.arange(0, 50, 0.1)  # 时间步长
init_i = 2  # 初始感染者


def sis_i(i, t):
    """
    SIS模型
    di(t)/dt = lambda * i(t) * (1 - i(t)) - mu * i(t)

    :param i: 感染者
    :param t: 时间
    :return i: t时刻的i
    """
    i = lamb * i * (1 - i) - mu * i

    return i


def sis_s(s, t):
    """
    SIS模型
    ds(t)/dt = -lambda * s(t) * (1 - s(t)) - mu * (1 - s(t))

    :param s: 易感者
    :param t: 时间
    :return s: t时刻的s
    """
    s = lamb * s * (1 - s) - mu * (1 - s)

    return s


def calc_init_ratio(nodes):
    """
    计算初始时刻的i(0)

    :param nodes: 节点
    :return: 第0时刻的感染率
    """
    n = len(nodes)
    return init_i / n


def calc_ODE(nodes):
    """
    计算常微分方程

    :param nodes: 节点
    """
    i = calc_init_ratio(nodes)
    result_i = spi.odeint(func=sis_i, y0=i, t=t_range)
    result_s = spi.odeint(func=sis_s, y0=i, t=t_range)
    draw_pic(result_i, result_s)


def draw_pic(result_i, result_s):
    """
    绘制感染情况图
    :param result_i: 感染者的比率构成的向量
    :param result_s: 易感者的比率构成的向量
    """
    plt.plot(t_range, result_i)
    plt.plot(t_range, result_s)
    plt.xlabel("time")
    plt.ylabel("ratio")
    plt.show()


if __name__ == '__main__':
    calc_ODE(nodes)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77

2. 基于节点的SIS模型

以下内容参考的是《节点（个体）级SIS模型python代码》

代码和上文的基本一致，我这里只是讲解下这位大佬的代码的含义，再加上了一条曲线。

2.1 微分方程

def sis_i(i, t):
    """
    节点级的计算sis
    这里计算了每一轮的i(t)值, 将i(t)代入到方程中计算出t这个时刻的值, 并赋值到I向量中,
    记录这一个时刻的节点情况
    di/dt = (1-i)λn - μi
    :param i: 感染者(向量)
    :param t: 迭代次数(时间)
    :return i: 每一时刻的节点的感染情况
    """
    I = np.zeros(n)
    for j in range(n):
        # 计算每个节点的邻居的感染情况
        neighbor = 0
        for k in range(n):
            neighbor += adj[j, k] * i[k]
        I[j] = (1 - i[j]) * lamb * neighbor - mu * i[j]

    return I
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

这串代码最开始对于我这个没有接触过scipy的人来说有点头痛，但是我查了下资料后自我总结了一番，可能有问题，欢迎大家指正。

首先我们要明确两点：

1. 这个微分方程传入的参数$i$是一个向量，该向量每一维代表的是整个网络中的某一个节点；
2. 这个微分方程是针对节点的，那么说明该向量中每一维的值代表着这个节点受感染的概率。

整个微分方程我们就着重于一行代码：

I[j] = (1 - i[j]) * lamb * neighbor - mu * i[j]
1

这个I和neighbor是什么我们一会说，等号后面这一坨就是微分方程，将其转换为公式为：
$$\frac{di(t)}{dt}=(1-i(t))\lambda n-\mu i(t)$$

我们回忆下上面的公式：$s(t)=1-i(t)$，那么公式可以改为：
$$\frac{di(t)}{dt}=\lambda ns(t)-\mu i(t)$$

那么再说neighbor，我们先不考虑neighbor具体代表什么，我们就当做是邻居节点，这里我们用$s(t{)}_j$表示该事件该节点是易感者的概率，那么$\lambda ns(t)$就代表着该节点是易感者且被邻居感染的概率。

$\mu i(t)$就很好理解，就是该节点自行治愈的概率。

那么整个微分方程代表的变化率就是说该节点被其他节点传染的概率，再减去自己治愈的概率，得到该节点最终的概率。

接着我们再来说neighbor，neighbor可以看做是这个公式：${\sum }_{j=1}^n{\sum }_{k=1}^{n}ad{j}_{jk}\ast i(t{)}_k$，这里$adj$指的是邻接矩阵，那么就是说求的是$j$节点在$t$时刻的邻居的状态。

而这里至于$I$是什么，我查了很多资料，发现好像是只要返回向量等，都需要通过这种重新赋值的方式来进行操作。

2.2 另一条曲线

根据上面的代码，我重新构建了感染者情况的微分方程：

def sis_s(s, t):
    """
    节点级的计算sis
    di/dt = -sλn - μ(1-s)
    :param s: 易感者(向量)
    :param t: 迭代次数(时间)
    :return i: 每一时刻的节点的易感情况
    """
    S = np.zeros(n)
    for j in range(n):
        # 计算每个节点的邻居的感染情况
        neighbor = 0
        for k in range(n):
            neighbor += adj[j, k] * (1 - s[k])
        S[j] = -s[j] * lamb * neighbor + mu * (1 - s[j])

    return S
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

这里为什么neighbor还是${\sum }_{j=1}^n{\sum }_{k=1}^{n}ad{j}_{jk}\ast i(t{)}_k$这个公式，是因为我们微分方程中前面这一坨始终是针对感染者的，所以这里也需要针对感染情况。

2.3 结果展示

蓝色代表的是每个时刻每个节点被感染的平均概率，橘色代表的是每个时刻每个节点被治愈（变为易感者）的平均概率。

2.4 完整代码

import scipy.integrate as spi
import numpy as np
import matplotlib.pyplot as plt
import random as rd
from data import static_para as ds
import networkx as nx
import itertools


def create_direct_graph(nodes):
    dg = nx.DiGraph()
    edges = itertools.permutations(nodes, 2)

    dg.add_edges_from(edges)

    return dg


def create_matrix(dg):
    """
    生成邻接矩阵
    :param dg: 有向图
    :return: 邻接矩阵
    """
    return nx.to_numpy_matrix(dg)


nodes = np.arange(0, 5, 1)
lamb = ds.lamb
mu = ds.mu
t_range = ds.t_range
dg = create_direct_graph(nodes)
adj = create_matrix(dg)
n = len(nodes)

i = np.zeros(n)

i_1 = rd.randint(0, n-1)
i_2 = rd.randint(0, n-1)
while i_1 == i_2:
    i_2 = rd.randint(0, n-1)
i[i_1] = 1
i[i_2] = 1


def sis_i(i, t):
    """
    节点级的计算sis
    这里计算了每一轮的i(t)值, 将i(t)代入到方程中计算出t这个时刻的值, 并赋值到I向量中,
    记录这一个时刻的节点情况
    di/dt = (1-i)λn - μi
    :param i: 感染者(向量)
    :param t: 迭代次数(时间)
    :return i: 每一时刻的节点的感染情况
    """
    I = np.zeros(n)
    for j in range(n):
        # 计算每个节点的邻居的感染情况
        neighbor = 0
        for k in range(n):
            neighbor += adj[j, k] * i[k]
        I[j] = (1 - i[j]) * lamb * neighbor - mu * i[j]

    return I


def sis_s(s, t):
    """
    节点级的计算sis
    di/dt = -sλn - μ(1-s)
    :param s: 易感者(向量)
    :param t: 迭代次数(时间)
    :return i: 每一时刻的节点的易感情况
    """
    S = np.zeros(n)
    for j in range(n):
        # 计算每个节点的邻居的感染情况
        neighbor = 0
        for k in range(n):
            neighbor += adj[j, k] * (1 - s[k])
        S[j] = -s[j] * lamb * neighbor + mu * (1 - s[j])

    return S


def calc_ODE():
    """
    计算微分方程
    :return:
    """
    result_1 = spi.odeint(func=sis_i, y0=i, t=t_range)
    result_1_mean = np.mean(result_1, axis=1)
    result_2 = spi.odeint(func=sis_s, y0=1-i, t=t_range)
    result_2_mean = np.mean(result_2, axis=1)
    draw_result(result_1_mean, result_2_mean)


def draw_result(result_1_mean, result_2_mean):
    """
    绘制SIS曲线
    :param result_mean: 结果向量
    """
    plt.plot(t_range, result_1_mean)
    plt.plot(t_range, result_2_mean)
    plt.xlabel("time")
    plt.ylabel("ratio")
    plt.show()


if __name__ == '__main__':
    calc_ODE()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111

3 参考

[1]7im0thyZhang.节点（个体）级SIS模型python代码[EB/OL].https://blog.csdn.net/timothytt/article/details/79718114,2018-3-27.