PythonER随机网格&数据可视化实现"以SEIR传播动力学模型实验为例"

SEIR传播动力学模型

1. 什么是SEIR传播动力学模型？
   经典SEIR模型将人群分为S (Susceptible，易感者)，I (Infected，感染者)，E (Exposed，潜伏者)和R (Recovered，康复人群)。

2. 本文根据SEIR传播动力学模型，利用Python实现生成ER随机网格、并在网格间随机结边、模拟SEIR型传染病传播

3. 该实验可分为两个步骤进行：
   一、生成一个ER随机网络，生成网络节点的度分布图、画出节点分布图
   二、编写程序：编程SEIR传播动力学模型。随机选择一个传播节点，基于SEIR传播动力学模型，统计每个时间步长对应的S,E,I,R四个状态节点数，作图
   *注意：为防止概率传播误差，宜进行100次独立实验，取平均值

4. 关键词解析：
   • 网络是由节点和边组成。节点代表个体，边代表节点之间的关系，比如两个节点是朋友，那么他们之间有连边，网络可以用一个邻接矩阵表示，矩阵中的0表示没有接触，即没有连边，1表示有接触，即存在连边。下面是一个简单的例子：
   
   • 度（degree）: 邻居节点的个数，如节点a的度为4，对应于邻接矩阵中a行或列1的个数。

5. 生成ER随机网络—— 利用 numpy 函数
   初始: 网络总结点数 N, 没有边；
   网络构建过程：每对节点以概率 p 被选择，进行连边，不允许重复连边。代码如下：

def create_ER(N, p):                #利用numpy生成N*N[0，1]的随机二维数组，再用np（where）条件赋值
    ERmap = np.random.rand(N, N)
    ERmap = np.where(ERmap>p, 0, 1)
    for i in range(N):
        ERmap[i, i] = 0
        ERmap[i, :] = ERmap[:, i]   #ER邻边连接随机网格，应具有对称性、对角线为0
    return ERmap

这里继续将ER随机网格打印为文件，代码如下：

def save_ER_to_file(ERmap):         #将ER邻边连接随机网格写入txt文档
    file = open('/Users/apple/Documents/Class Folders/2019-2020 2/Python programming /Codes/Python/'
                'SEIR 传播模型/ER随机网格.txt', 'w+')
    for i in range(len(ERmap)):     #打开文档，若无将自动创建
        a = ERmap[i]
        for j in range(len(ERmap)):
            s = str(a[j])           #先取第一行，再取各列
            file.write(s)
            file.write(' ')
        file.write('\n')
    file.close()

接着将ER随机网格转为ER随机分布图，代码如下：

def showGraph(ERmap):                 #将生成的ER随机网格连接生成分布图
    G = nx.Graph()
    for i in range(len(ERmap)):
        for j in range(len(ERmap)):
            if ERmap[i][j] == 1:
                G.add_edge(i, j)    #将值为1的点相连接
    nx.draw(G)
    plt.savefig('/Users/apple/Documents/Class Folders/2019-2020 2/Python programming /Codes/Python/'
                'SEIR 传播模型/ER网格图.png')
    plt.show()                      #保存图片、显示图片

现在已经成功生成了ER随机分布图，我们接下来完成度分布概率图（并顺手把概率分布存入txt文件）：

def calculateDegreeDistribution(ERmap):
    avedegree = 0.0             #计算度分布
    identify = 0.0                  #若算法正确，度概率分布总和应为0
    p_degree = np.zeros((len(ERmap)), dtype=float)
                                    #statistic下标为度值
                                    #(1)先计数该度值的量
                                    #(2)再除以总节点数N得到比例
    degree = np.zeros(len(ERmap), dtype=int)
                                    #degree用于存放各个节点的度之和
    for i in range(len(ERmap)):
        for j in range(len(ERmap)):
            degree[i] = degree[i] + ERmap[i][j]
                                    #汇总各个节点的度之和
    for i in range(len(ERmap)):
        avedegree += degree[i]
                                    #汇总每个节点的度之和
    print('该模型平均度为\t' + str(avedegree /len(ERmap)))
                                    # 计算平均度
    for i in range(len(ERmap)):
        p_degree[degree[i]] = p_degree[degree[i]] + 1
                                    #先计数该度值的量
    for i in range(len(ERmap)):     #再除以总节点数N得到比例
        p_degree[i] = p_degree[i] / len(ERmap)
        identify = identify + p_degree[i]
                                    #将所有比例相加，应为1
    identify = int(identify)
    plt.figure(figsize=(10, 4), dpi=120)
                                    #绘制度分布图
    plt.xlabel('$Degree$', fontsize=21)
                                    #横坐标标注——Degrees
    plt.ylabel('$P$', fontsize=26)
                                    #纵坐标标注——P
    plt.plot(list(range(len(ERmap))),list(p_degree),'-*',markersize=15,label='度',color = "#ff9c00")
                                    #自变量为list(range(N)),因变量为list(p_degree)
                                    #图形标注选用星星*与线条-，大小为15，标注图例为度，颜色是水果橘

    plt.xlim([0, 12])               #给x轴设限制值，这里是我为了美观加入的，也可以不写
    plt.ylim([-0.05, 0.5])          #给y轴设限制值，也是为了美观加入的，不在意可以跳过
    plt.xticks(fontsize=20)         #设置x轴的字体大小为21
    plt.yticks(fontsize=20)         #设置y轴的字体大小为21
    plt.legend