数学建模社区-数学中国
标题:
Python实现简单的SI传播模型
[打印本页]
作者:
杨利霞
时间:
2020-4-18 16:16
标题:
Python实现简单的SI传播模型
Python实现简单的SI传播模型
9 l) _1 o: A2 z [, h
#SI疾病传播模型的原理
9 ^! X: `7 n4 F4 Q1 G# ], \
在经典的传染病模型中,种群(Population)内N个个体的状态可分为如下几类
/ \, O" W$ i4 ]
6 ~2 [$ e; C) `2 O3 H; p x
易感状态(Susceptible)。一个个体在感染前是处于易感状态的,即该个体有可能被邻居个体感染。
% d: a: { Q/ X+ Q0 A3 X7 f( o
易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。,即该个体还会以一定概率感染其邻居个体。
9 Z' {: P3 X/ D9 B
移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态,当一个个体经历过一个完整的感染周期后,该个体就不再被感染,因此就可以不再考虑改革提。
* t* r0 w. ]- A3 ~: U
SI传播模型是最简单的疾病传播模型,模型中的所有个体都只可能处于两个状态中的一个
5 ~0 b, ] D' _. A, L" q
即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
$ ^" p. }6 _6 R$ l- W
在给定时刻t,令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目,显然有
; e3 s ]# T( U. ~" y
S(t)+I(t)恒等于N,这里,N是个体总数。随着时间t的增长,易感个体与感染个体的接触
, y# f/ m) W/ s, b B5 a. w
会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β,疾病仅在
0 h2 P( t2 ^- w6 h. r* [2 R( C
感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t,易感个体的比例为S(t)/N,感染个体的数量为I(t),一次,易感个体的数量将以如下变化率减少
- q2 ]# R$ g* O. r' k
ds/dt = -β*S(t)I(t)/N
! y1 x% q8 a: h: S
同时,感染个体的数量会以与易感个体相反的变化率增加,
' q2 c. U8 n! M- j' Y) _
ds/dt = βS(t)*I(t)/N
" E- H" h/ O& g
分别将时刻t处于易感状态和感染状态的个体所占比例记为,
4 Y8 {: J6 K/ w! v* z
s(t)=S(t)/N
8 h G! {# P; R% B2 Z7 E
i(t)=I(t)/N
/ |5 D4 i6 [$ i/ i% w$ x1 N/ Q
显然有,
# F5 K' j; d: z
s(t)+i(t)恒等于1,此时之前的公式可以记做
! S# R2 U b: p# V; U2 }
ds/dt=-βsi
: O: l% \ i# S( @
di/dt=βsi
, o: p0 J" N* V2 H( {8 p
即
" K) ~1 w$ M9 ?0 J& B9 [) K5 J
di/dt=βi(1-i)
$ W; p- I! F9 F: x+ T
上式也成为Logistic增长方程式(Logistic growth equation),
) F; v! g& C7 S
方程的解和图像如图
! S) r! t0 `- Q, j# ]- q
2020-4-18 16:14 上传
下载附件
(213.53 KB)
. l' k8 A5 Z) {
代码和相关文件以及环境链接:链接:https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ
/ H( w; e2 u1 r) C( V( g, d! ~
提取码:z448
! g1 J: k3 ~, |9 t0 B2 L+ s
; s5 u0 t8 G. v, |5 Y
7 M# U; b X. a/ l" J5 r: ]0 j
'''
1 o9 o- F1 Q8 Q% y. {$ e3 c: H
实验环境Python2.7.13,igraph包,cairo包,numpy包
4 N" ?# b! a8 L1 U9 l3 S7 J
'''
( [0 ]4 n* E* C, ~
# -*- coding:utf8 -*
3 P y* E' u+ J$ M- S" h
from igraph import *
3 f1 O6 j: h4 k% g; D' o, J
import numpy as numpy
0 W6 M. \- q/ U3 E
from numpy import *
- i4 k: C v8 S: {! V2 T; E4 h
import random
* _* c: C; f3 ^; p# J& x
* H; F) m. {) L& x# p+ B, c( m
def len_arr(infected_array,nodes_num):#获取感染数组长度
4 c, {- w3 J! t
len_value=0#初始化长度
! n+ f5 v# M4 p& Y. U/ h
len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
! v6 g& ?4 U4 x* i2 L* q. T4 P7 @
return len_value
3 U% b/ c% i' J% v0 E9 m
2 M1 ]" ~7 g! k( u
g=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
9 L, m$ Q8 @; N. S" `
summary(g)
9 g' d: U+ A6 S$ |; B& I# E
nodes_num=g.vcount()#统计图中的结点个数
" z4 y) g: f$ p& T! W) @! b! G
net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
% d) f6 {8 d: e2 B0 f0 ~
g.vs["color"]=["white"]#给图的顶点序列颜色赋值白色
4 f/ [6 U& X, b+ `/ u4 r9 @
a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a
" N; @( J. L. }/ u; I
nodes_state=matrix(a).T#nodes_state通过转置a矩阵创建,用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重!!!)
2 C& t% D. h8 O; `
#第一列是节点编号,第二列是节点状态,感染状态用-2表示,第三列是节点感染的时间
9 v9 P! I# L2 \9 S- C% Q, \) J
print(nodes_state)
% V: z8 }4 z4 S2 ]; E8 R E$ f
infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染 34代表网络节点数
2 d/ s3 F3 a6 V$ _$ X* Q6 P9 n
print(infected_array)
9 j5 o0 `$ z2 j
' E( P$ ^" b% E( F
infe_rate=1#传播率(感染率) 1代表邻接点100%被感染
5 P t0 I7 @4 y, W
set_time=2#传播次数(感染次数) 2次
: O. T* c7 J9 v
source_seed=1#感染源位置
& N5 |0 i8 B: }/ Z$ q
nodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
* J* V3 R8 S) O6 j
nodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
5 V1 k8 @6 _; r; z4 |0 P
nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
, q! {! r# N( f+ d! q
g.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
& H$ t; c) b7 e: h9 C9 q
infected_array[0]=source_seed#将感染源的位置存入被感染节点列表
1 J% _" @1 G) k) l5 {8 u) D+ U& ?
plot(g)#绘制
: c6 [) ?" ?6 t' C6 T" I+ _
4 _8 K7 z/ G! Y! `2 K
stop=False#感染过程结束的标记
% F' _/ ]- @0 F; {
temp_time=0#第几次感染
, ]8 K" D8 g1 m$ v: H$ S
temp_len=0#本轮的感染源数量初始化
' N/ [% r# T* \0 L! V
3 P% J# A- s. p" n4 p" i
while not stop:
r5 [; K; g: M/ ]- _, b& W( z6 K
i=0#记录让每个感染源都传播一次
5 ~* ?; \, D, I6 D7 B5 Y
if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行
: E6 K4 D D+ X F8 A7 K: @' q1 L
temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量
4 V4 s" q" D2 x
while i<temp_len:
! ?+ l8 F3 l9 V; Y
temp_time=nodes_state[infected_array
-1,2]#获取每一个节点的感染时间
5 M/ d! v+ c5 l G3 Y
nei_count=0#下一轮可以被感染到的节点数量
) f1 R" J4 g: r2 F
#生成下一轮可能被感染的节点的集合nei_arr
5 Z6 R+ t& r/ D! {, h% o- e
for j in range(nodes_num):#遍历节点
3 k p! U& m: i0 p8 S2 N- B# C
if net_mat[infected_array
-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
B; f g- T8 R- u( Q' E5 L
nei_count=nei_count+1#下一轮可以被感染到的节点数量++
* L# `+ C, h Y g. Y9 P" [ D
nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
0 A4 b+ e# t. {; R0 C
t=0
* G3 A: v g! d$ B$ i8 r
for j in range(nodes_num):
; r. m! b( u8 R k; r, {
if net_mat[infected_array
-1,j]==1 and nodes_state[j,1]!=-2:
8 a3 t% Z! {* E Y4 l+ N; D, ^
nei_arr[t]=j+1
7 f$ L" e. u2 {; f8 [
t=t+1
3 H) F7 F. J& u2 L- [
ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
: n0 j* Q( P2 x, N
#random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象
1 u. i; \: ~+ o# }4 M6 R
if len(ran_infe_arr)>0:#存在需要被感染的节点
6 r/ \( v7 h9 I7 S0 g6 f; f
t=0#让ran_infe_arr内每个感染源都被感染
/ A7 \- j+ W2 g/ I; e; ]
while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
. O6 t% L- i5 [; L
nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
- k" x; i) V5 C& V) E
nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间
: }: f3 v7 O: Z1 @8 h4 i
infected_array[len_arr(infected_array,nodes_num)]=nei_arr[ran_infe_arr[t]]#将此次感染节点放入总的感染节点数组中
7 `9 ^0 j0 O7 @ R
g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色
" h% s( ~4 O, \( l
plot(g)#绘制
. x# w. |& a$ h2 q ]- p) }6 x
t=t+1
5 A4 j0 m1 A" L% i7 h- M
i=i+1
$ s5 N* f8 h# w# _
if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
/ g# ~9 |2 F% W2 M/ \6 W+ w
stop=True
! [. k4 ~% q7 }$ t( a% B5 y3 D
/ W* H* S) {) y4 I+ _
U% E/ ~/ `( _. o# Y: v
视频演示bilibili传送门
- V& [ F1 G9 Y' ^& O
效果图
3 S/ o) Q0 x: M7 O J4 \
9 R+ s0 m9 f. X, m2 j
* x) |) B: f) }+ o% G2 n$ w
2020-4-18 16:15 上传
下载附件
(174.22 KB)
S z8 R* F) e [3 i
6 B/ a& x6 Q! s7 M) H1 v
2020-4-18 16:15 上传
下载附件
(22.44 KB)
+ T# q- z: Z4 v7 L3 J
# Y' `6 }4 B! g- u+ o: P# x- b1 `
2020-4-18 16:16 上传
下载附件
(42.85 KB)
* m) F6 q% X x
, l6 j. Z: F+ N7 L$ q& O
2020-4-18 16:16 上传
下载附件
(17.44 KB)
& J! ?! J3 P% `2 w$ @8 s. R: u
2 f, Y, U# U/ u" i3 v
2020-4-18 16:16 上传
下载附件
(19.78 KB)
1 H$ l( b' y d7 i' k0 C' \
, ?: m( x) r. Q' q
2020-4-18 16:16 上传
下载附件
(21.16 KB)
————————————————
7 B; Y- j, b+ h4 w6 T
版权声明:本文为CSDN博主「eck_燃」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。
9 J( Q. p1 j( f& V4 C. n* P
原文链接:https://blog.csdn.net/wdays83892469/article/details/80878862
: h- w4 u& b5 A, n4 `( j: U; j9 F
" f* |$ |; c- W. K6 ~& a- o, J* b: C
7 n2 C% H) _% J m; k L4 J$ c
作者:
尔雅
时间:
2020-4-18 17:50
发表回复棒!
- {( x" N/ q2 p8 a/ S; ]! @: y
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5