De-Sim示例分析（三）SIR传染病模型

De-Sim前两个示例介绍了一类message、多类message、一个仿真对象、对个仿真对象的使用方法，De-Sim的第三个示例传染病SIR模型介绍了检查点对象的使用方法，通过给检查点对象设置记录频率和记录内容，可以实时记录仿真对象的情况。

SIR模型将传染病传播中人的状态分为三类Susceptible、Infectious、Recovered，分别表示易受感染、患病、痊愈。SIR通过S向I的转变、I向R的转变，具体表现为数值的变化，来预测传染病的传播情况。
该过程在具有马尔科夫性质，设在时间内的转移概率为

则由ODE model给出的转移概率为

以上过程如何推导可以不了解，只需要知道该模型的结论，定义以下符号

每经过τ的时间，P1的概率有1个S转化为I，P2的概率有1个I转化为R

一个SIR过程用折线图可以表示

该示例由枚举（Enum）类型设置了两类message，表示S转化为I和I转化为R。
普通类SIR，该类的属性s/i/N/beta/gamma是SIR模型的参数，random_state为随机数种子，可以用于仿真过程的再现
该类有两个方法
·handle_state_transition(self,transition_type)，用于根据传入的transition_type更新SIR的s和i的数量
·plan_next_transition(self)，根据当前SIR的参数，计算转化的概率并确定转化的类型和时间tau，返回tuple类型
仿真对象类SIRSimutor 该类在构造函数内示例化了SIR类该类有四个方法 ·init_before_run(self),初始化该类
·handle_event(self,event),该方法调用sir的handle_state_transition方法，通过更新sir属性的方法记录了event的信息
·plan_and_create_event(self),该方法调用sir的plan_next_transition方法，确定了下一次转化的时间和类型
·schedule_next_event(self,planed_event),该方法执行plan_and_create_event中确定的转化时间和类型
类内方法外的属性定义和前两个示例一致
访问状态类该类是AccessStateObjectInterface的实现
该类在构造函数内传入了sir对象，也就是在SIRSimutor的构造函数内创建的对象，通过记录该对象的属性实时记录仿真状态两个方法
·get_checkpoint_state(self,time)，返回要记录的信息
·get_random_state(self)，返回随机数种子，用于仿真过程的复现
普通类RunSIR 该类整合了启动仿真和仿真结果展示的方法该类有三个方法
·simulate（self，recording_period,max_time,**sir_args）
该方法创建了simulator，为simulator添加了仿真对象和检查点对象，并运行了max_time
·last_checkpoint(self),用于返回最后一个检查点的信息
·history2dataframe(self),将检查点中的信息转成dataframe的格式

添加中文注释的代码如下:
(原始代码使用tempfile暂存了输出，我直接plt.show())

import de_sim
import numpy
import enum
import pandas
import matplotlib.pylab as plt
import sysfrom de_sim.checkpoint import AccessCheckpoints
from de_sim.simulation_checkpoint_object import (AccessStateObjectInterface,CheckpointSimulationObject)#设置转换类型
class StateTransitionType(enum.Enum):s2i='Transition from Susceptible to Infectious'i2r='Transition from Infectious to Recovered'#通过转换类型设置事件message
class TransitionMessage(de_sim.EventMessage):#指示应该发生哪个转换，属性将存储一个实例"Messge for all model transitions"transition_type:StateTransitionType#创建SIR类，该类将个体在SIR状态下的转换表示为离散事件，该类将循环访问各事件
class SIR(object):"""SIR每次迭代执行以下操作接收事件->执行事件message中的状态转换->规划下一个message->安排下一次事件"""def __init__(self,s,i,N,beta,gamma):"""Args:s:int,易感染者S初始数量i: int,感染者I初始数量N: int,总数beta: SIR beta参数gamma: SIR gamma参数"""self.s=sself.i=iself.N=Nself.beta=betaself.gamma=gammaself.random_state=numpy.random.RandomState()#设置随机数种子def handle_state_transition(self,transition_type):"""更新SIR状态Args:transition_type: 接收的信息"""if transition_type is StateTransitionType.s2i:self.s-=1self.i+=1elif transition_type is StateTransitionType.i2r:self.i-=1def plan_next_transition(self):"""计划下一次事件的时间及messageReturns:tuple:(delay,event message) if an event is scheduled,otherwise None"""rates={'s2i':self.beta*self.s*self.i/self.N,'i2r':self.gamma*self.i}lambda_val=rates['s2i']+rates['i2r']sys.exit()if lambda_val==0:return Nonetau=self.random_state.exponential(1.0/lambda_val)prob_s2i=rates['s2i']/lambda_valprint(prob_s2i,'****')if self.random_state.random_sample()<prob_s2i:return (tau, TransitionMessage(StateTransitionType.s2i))else:return (tau, TransitionMessage(StateTransitionType.i2r))#setattr(SIR, 'handle_state_transition', handle_state_transition)#创建和调度事件，接收并执行事件
class SIRSimulator(de_sim.SimulationObject):#通过继承SimulationObject来使用功能，obtain SIR logic by composing a instancedef __init__(self,name,s=98,i=2,N=100,beta=0.3,gamma=0.15):self.sir=SIR(s,i,N,beta,gamma)super().__init__(name)def init_before_run(self):#Initialize before a simulation run.Send the first eventself.plan_and_create_event()def handle_event(self,event):#调用SIR的方法根据传来的信息更新SIR的属性self.sir.handle_state_transition(event.message.transition_type)self.plan_and_create_event()def plan_and_create_event(self):#Plan and create the next event#调用SIR的方法按给定参数确定的概率确定传递出去的messageplanned_event=self.sir.plan_next_transition()if planned_event is not None:self.schedule_next_event(planned_event)def schedule_next_event(self,planned_event):#发送messagetau,transition_message=planned_eventself.send_event(tau,self,transition_message)#该主体收到TransitionMessage消息，就调用handle_event函数处理event_handlers=[(TransitionMessage,'handle_event')]messages_sent=[TransitionMessage]#创建检查点定期保存模拟的状态，以便后续的分析
#模拟的状态被概念化为一组值及其随机数生成器的状态
class AccessSIRObjectState(AccessStateObjectInterface):"""Get the State of a SIR objectAttributes:sir:a sir Objectrandom_state:a random state"""def __init__(self,sir):self.sir=sirself.random_state=sir.random_state#这两个函数是必要的，一个返回要记录的信息，一个返回此时的随机数种子def get_checkpoint_state(self, time):"""Get the SIR simulation stateArgs:time:current time;ignoredReturns:dict:sir state"""return dict(s=self.sir.s,i=self.sir.i)def get_random_state(self):"""Get the SIR object random stateReturns:random state"""return self.random_state.get_state()#运行SIR模型并可视化其预测的代码
#需要给定名称、检查点周期、保存检查点的目录以及可以访问模拟状态的对象
class RunSIR(object):def __init__(self,checkpoint_dir):self.checkpoint_dir=checkpoint_dirdef simulate(self,recording_period,max_time,**sir_args):"""Create and run a SIR simulationArgs:recording_period: float,interval between state checkpointsmax_time: float,simulation end time**sir_args: arguments for a SIR object"""#create a simulatorsimulator=de_sim.Simulator()#create a SIRSimulator instanceself.sir_simulator=SIRSimulator(**sir_args)simulator.add_object(self.sir_simulator)#create a checkpoint simulation objectaccess_state_object=AccessSIRObjectState(self.sir_simulator.sir)#CheckpointSimulationObject的四个参数#name/记录周期recording_period/checkpoint文件的保存地址/记录内容及格式access_state_objectcheckpointing_obj=CheckpointSimulationObject('checkpointing_obj',recording_period,self.checkpoint_dir,access_state_object)#为simulator添加检查点对象simulator.add_object(checkpointing_obj)#initialize simulation,which sends the SIR instance an initial event messagesimulator.initialize()#run the simulationevent_num=simulator.simulate(max_time).num_eventsdef last_checkpoint(self):"""Get the last checkpoint of the last simulation runReturns:obj checkpoint:the last checkpoint of the simulation"""#根据地址读取checkpoint文件access_checkpoints = AccessCheckpoints(self.checkpoint_dir)last_checkpoint_time=access_checkpoints.list_checkpoints()[-1]return access_checkpoints.get_checkpoint(time=last_checkpoint_time)def history2dataframe(self):"""Convert the checkpoint history to a Pandas dataframeReturns:pandas.DataFrame:columns are SIR values,rows are timepoints"""fields=('s','i','r')hist=[]index=[]#读取所有检查点添加的时间，返回形如[0,10,20,30,...]样式的列表access_checkpoints=AccessCheckpoints(self.checkpoint_dir)for checkpoint_time in access_checkpoints.list_checkpoints():#对于每个检查点，读取其状态，其类型在AccessSIRObjectState内被指定为字典state=access_checkpoints.get_checkpoint(time=checkpoint_time).state#转换成列表state_as_list=[state['s'],state['i'],self.sir_simulator.sir.N-state['s']-state['i']]hist.append(dict(zip(fields,state_as_list)))index.append(checkpoint_time)return index,pandas.DataFrame(hist)#转换成dataframeif __name__=='__main__':sir_args=dict(name='sir',s=98,i=2,N=100,beta=0.3,gamma=0.15)run_sir = RunSIR('./checkpoint_dir')for i in range(10):run_sir.simulate(10,100,**sir_args)#print and plot an epidemic predicted trajectoryindex,sir_data_frame=run_sir.history2dataframe()plt.plot(index,sir_data_frame['s'],label='s')plt.plot(index,sir_data_frame['i'],label='i')plt.plot(index, sir_data_frame['r'],label='r')plt.ylabel('population')plt.xlabel('time')plt.legend()plt.title('epidemic model%d'%i)#plt.savefig('epidemic %d.jpg'%i)plt.show()

以下是我保留基础功能的简化版本

import de_sim
import numpy
import matplotlib.pylab as plt
import sysfrom de_sim.checkpoint import AccessCheckpoints
from de_sim.simulation_checkpoint_object import (AccessStateObjectInterface,CheckpointSimulationObject)
class TransformMessages2i(de_sim.EventMessage):"Transition from Susceptible to Infectious"class TransformMessagei2r(de_sim.EventMessage):"Transition from Infectious to Recovered"class AccessStateObject(AccessStateObjectInterface):def __init__(self, sir):self.sir = sirself.random_state = sir.random_statedef get_checkpoint_state(self, time):return dict(s=self.sir.s, i=self.sir.i,r=self.sir.r)def get_random_state(self):return self.random_state.get_state()class SIR(de_sim.SimulationObject):def __init__(self,args):self.random_state = numpy.random.RandomState()self.beta=args['beta']self.gamma=args['gamma']self.s=args['s']self.i=args['i']self.r=args['N']-args['s']-args['i']self.N = args['N']super().__init__(args['name'])def init_before_run(self):self.schedule_next_event()def handle_step_event(self,event):if isinstance(event.message,TransformMessages2i):self.s-=1self.i+=1elif isinstance(event.message,TransformMessagei2r):self.i-=1self.r+=1self.schedule_next_event()def schedule_next_event(self):lambda_val = self.beta * self.s * self.i / self.N + self.gamma* self.iif lambda_val==0:print('None')return Noneprob_s2i = (lambda_val -self.gamma*self.i)/ (lambda_val)tau=self.random_state.exponential(1.0/lambda_val)print(prob_s2i,'***')if self.random_state.random()<prob_s2i:self.send_event(tau,self,TransformMessages2i())else:self.send_event(tau, self, TransformMessagei2r())event_handlers=[(TransformMessages2i,'handle_step_event'),(TransformMessagei2r,'handle_step_event')]messages_sent=[TransformMessages2i,TransformMessagei2r]class RunSir(object):@staticmethoddef simulate(max_time,sir_args):simulator = de_sim.Simulator()sir_obj = SIR(sir_args)simulator.add_object(sir_obj)access_state_object = AccessStateObject(sir_obj)che_obj = CheckpointSimulationObject(name='checkpoint_obj',checkpoint_dir='./checkpoint_dir',checkpoint_period=10,access_state_obj=access_state_object)simulator.add_object(che_obj)simulator.initialize()simulator.simulate(max_time)@staticmethoddef plot_result():index = []s_hist, i_hist, r_hist = [], [], []access_checkpoints = AccessCheckpoints('./checkpoint_dir')for checkpoint_time in access_checkpoints.list_checkpoints():# 对于每个检查点，读取其状态，其类型在AccessSIRObjectState内被指定为字典state = access_checkpoints.get_checkpoint(time=checkpoint_time).state# 转换成列表s_hist.append(state['s'])i_hist.append(state['i'])r_hist.append(state['r'])index.append(checkpoint_time)plt.plot(index, s_hist, label='s')plt.plot(index, i_hist, label='i')plt.plot(index, r_hist, label='r')plt.ylabel('population')plt.xlabel('time')plt.legend()plt.title('epidemic model')# plt.savefig('epidemic %d.jpg'%i)plt.show()if __name__=='__main__':sir_args = dict(name='sir',s=98,i=2,N=100,beta=0.3,gamma=0.15)for i in range(4):RunSir.simulate(100,sir_args)RunSir.plot_result()

De-Sim示例分析（三）SIR传染病模型相关推荐

最近疯传的SIR传染病模型是什么？
总第188篇/张俊红最近看到在网上传的一张SIR传染病模型的图,很多人应该对这个模型不是很了解,今天就讲一下这个模型.这一篇只讲学术,不讨论别的. SIR模型是传染病模型中最经典的一个,类似的还有S ...
SIR传染模型Matlab代码,sir传染病模型 MATLAB代码运行不了,
问题描述: sir传染病模型 MATLAB代码运行不了, function y=ill(t,x) a=1;b=0.3; y=[a*x(1)*x(2)-b*x(1),-a*x(1)*x(2)]'; ts ...
SIR传染病模型介绍+R语言简单应用
此文用于整理回顾写论文时看的文献资料和学到的知识,也希望能带来一些参考. 什么是SIR模型? 1927年,Kermack 和 McKendrick 为了研究17世纪肆虐伦敦的黑死病和20世纪席卷孟买的 ...
Python学习：SIR传染病模型
SIR模型是传染病模型中最经典的一个,类似的还有SI和SIS两种.SIR是三个单词首字母的缩写,其中S是Susceptible的缩写,表示易感者:I是Infective的缩写,表示感染者,R是Remo ...
微分方程（人口预测与传染病模型）
一.定义微分方程:含导数或微分的方程微分方程的阶数:所含导数或微分的最高阶数,如y'''+2y''-2x=0是三阶微分方程微分方程的解:使得微分方程成立的函数例如y'-2x=0的解可以为x²或 ...
NetLogo学习笔记5 —— 物种与传染病模型
NetLogo学习笔记5 -- 物种与传染病模型 (模型经过一些修改,与标准SIR模型有些出入) 在上一篇文章,我们学习了随机选择.伪并发.ifelse和of语法.实现了用于解释隔离现象的谢林模型这 ...
TVM开发三个示例分析
TVM开发三个示例分析把自主生成的代码生成TVM 把自主生成的代码生成TVM 目录简介要生成C代码. 要生成任何其它图形表示. 实现一个C代码生成器实现[CodegenC] 运算符代码生成输 ...
【组合数学】排列组合 ( 多重集组合数示例 | 三个计数模型 | 选取问题 | 多重集组合问题 | 不定方程非负整数解问题 )
文章目录一.多重集组合示例二.三个计数模型排列组合参考博客 : [组合数学]基本计数原则 ( 加法原则 | 乘法原则 ) [组合数学]集合的排列组合问题示例 ( 排列 | 组合 | 圆排列 | ...
matlab求解常微分方程组/传染病模型并绘制SIR曲线
看了很多关于传染病模型的matlab程序,大都是绘制出两条曲线(I.S)的,本文最大的不同是绘出SIR三条曲线. 先给出SIR微分方程组函数文件: run的程序:
数模学习第三天--微分方程（传染病模型）
微分方程定义包含关于未知变量.未知变量的导数以及自变量的方程第一类为可分离变量的微分方程. 在高数中我们学过相关内容,如: 建模准备 1.运用已知规律列方程利用数学.力学.物理.化学等学科中的 ...

De-Sim示例分析（三）SIR传染病模型

De-Sim示例分析（三）SIR传染病模型相关推荐

最新文章

热门文章