|
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)
这样就可以避免传入多个值报错的问题了。
|