|
25鱼币
- #导入模块
- import numpy as np
- import random
- import math
- import matplotlib.pyplot as plt
- #使用蒙特卡洛方法模拟传染
- def transmission(p):
- '''p: float, 传染率
- 如果生成的随机数比p小,则被感染
- 反之则不被感染'''
- if random.random() < p:
- return True
- else:
- return False
- def SIS(N, infecteds, alpha_number, beta, gamma, tmin=0, tmax=100):
- '''
- N: 总人口数, int
- infetceds: 初始感染人数, int
- alpha_number: 传染人数,int, 每个人能接触到的人数
- beta: 感染率, float
- gamma: 恢复率, float
- tmin: 初始时间, int
- tmax: 最大时间步长, int
- '''
- # 初始状态所有人均为感染,用数组表示,为全零数组
- state = [0] * N # 创建一个长度为N的全零数组
- # 随机从数组选infected个人,将他们的状态改为1
- for i in range(infecteds): # 循环infecteds次
- index = random.randint(0, N - 1) # 随机选择一个数组下标
- state[index] = 1 # 将该下标对应的元素设为1
- # 创建两个空列表,分别统计当前0和1的人数
- S_list = [] # 创建一个空列表,用于存储易感者的人数
- I_list = [] # 创建一个空列表,用于存储感染者的人数
- # 循环从tmin到tmax,每次增加1
- for t in range(tmin, tmax + 1):
- # 统计当前0和1的人数,并添加到对应的列表中
- S = state.count(0) # 计算数组中0的个数,即易感者的人数
- I = state.count(1) # 计算数组中1的个数,即感染者的人数
- S_list.append(S) # 将S添加到S_list中
- I_list.append(I) # 将I添加到I_list中
- # 对每个感染者进行操作
- for i in range(N): # 循环遍历数组
- if state[i] == 1: # 如果当前元素为1,即当前人为感染者
- # 随机选择alpha_number个人接触,并判断是否传染或恢复
- for j in range(alpha_number): # 循环alpha_number次
- k = random.randint(0, N - 1) # 随机选择一个数组下标
- if state[k] == 0 and transmission(beta): # 如果该下标对应的元素为0,即易感者,并且以beta的概率返回True
- state[k] = 1 # 将该元素设为1,即感染
- if transmission(gamma): # 如果以gamma的概率返回True
- state[i] = 0 # 将当前元素设为0,即恢复
- return S_list, I_list # 返回两个列表
- def PlotAxes(ax, xlabel, ylabel, title, mode=False):
- '''
- Decorate the axes
- Parameters
- ----------
- ax : axes
- axes.
- xlabel : str
- set the xlabel.
- ylabel : str
- set the ylabel.
- title : str
- set the title.
- mode : bool, optional
- whether to show the legend. The default is False.
- Returns
- -------
- None.
- '''
- fontsize = 16
- font_label = {'family': "Arial", 'size': fontsize}
- n_legend = 12
- ax.set_xlabel(xlabel, fontdict=font_label)
- ax.set_ylabel(ylabel, fontdict=font_label)
- ax.set_title(title, loc='left', fontdict={'family': "Arial", 'size': fontsize})
- ax.tick_params(direction='out', which='both', length=4, width=1, pad=1, labelsize=n_legend)
- # ax.minorticks_on()
- if mode == True:
- ax.legend(loc='best', framealpha=0, fontsize=n_legend)
- # 调用SIS函数,传入参数,得到两个列表
- N = 1000 # 总人口数
- beta = 0.3 # 感染率
- gamma = 0.1 # 恢复率
- tmin = 0 # 初始时间
- tmax = 100 # 最大时间步长
- S_list, I_list = SIS(N, beta, gamma, tmin, tmax, infecteds=10, alpha_number=5)
- # 创建一个图形对象
- fig = plt.figure()
- # 添加一个子图,设置为3D模式
- ax = fig.add_subplot(111)
- # 调用PlotAxes函数,为子图添加标签、标题和图例
- PlotAxes(ax, 'Time', 'Population', 'SIS model', mode=True)
- # 生成一个时间序列,与两个列表对应
- t = np.arange(tmin, tmax + 1)
- # 在子图上绘制两条曲线,分别表示易感者和感染者的人数随时间的变化
- ax.plot(t, S_list, label='Susceptible')
- ax.plot(t, I_list, label='Infected')
- # 显示图形
- plt.show()
复制代码
在调用SIS函数时,传入的参数顺序有误。应该按照函数定义的顺序传入参数。正确的调用方式如下:
S_list, I_list = SIS(N, infecteds=10, alpha_number=5, beta=0.3, gamma=0.1, tmin=0, tmax=100)
这样就可以避免传入多个值报错的问题了。
|
最佳答案
查看完整内容
在调用SIS函数时,传入的参数顺序有误。应该按照函数定义的顺序传入参数。正确的调用方式如下:
S_list, I_list = SIS(N, infecteds=10, alpha_number=5, beta=0.3, gamma=0.1, tmin=0, tmax=100)
这样就可以避免传入多个值报错的问题了。
|