QQ登录

只需要一步,快速开始

 注册地址  找回密码
查看: 3071|回复: 1
打印 上一主题 下一主题

Python实现简单的SI传播模型

[复制链接]
字体大小: 正常 放大
杨利霞        

5273

主题

82

听众

17万

积分

  • TA的每日心情
    开心
    2021-8-11 17:59
  • 签到天数: 17 天

    [LV.4]偶尔看看III

    网络挑战赛参赛者

    网络挑战赛参赛者

    自我介绍
    本人女,毕业于内蒙古科技大学,担任文职专业,毕业专业英语。

    群组2018美赛大象算法课程

    群组2018美赛护航培训课程

    群组2019年 数学中国站长建

    群组2019年数据分析师课程

    群组2018年大象老师国赛优

    跳转到指定楼层
    1#
    发表于 2020-4-18 16:16 |只看该作者 |倒序浏览
    |招呼Ta 关注Ta
    Python实现简单的SI传播模型1 |! Y8 @6 ]) H5 T6 \0 n. ^; O
    #SI疾病传播模型的原理
    ; N1 v/ Q8 l& w在经典的传染病模型中,种群(Population)内N个个体的状态可分为如下几类
    0 D: g2 a9 s+ O: j) ?$ p+ K$ [. W5 x4 F! e' s- L
    易感状态(Susceptible)。一个个体在感染前是处于易感状态的,即该个体有可能被邻居个体感染。* u1 l. J" v$ Z" \" @
    易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。,即该个体还会以一定概率感染其邻居个体。
    ; H- c6 U" h9 s移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态,当一个个体经历过一个完整的感染周期后,该个体就不再被感染,因此就可以不再考虑改革提。
    ( ]' M( w* C' @6 b' HSI传播模型是最简单的疾病传播模型,模型中的所有个体都只可能处于两个状态中的一个& j3 F7 p" U/ ~% H& f
    即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。
    7 Z  T( e* T: _: L' J$ A在给定时刻t,令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目,显然有
    % z' W. V0 e: K! C; n, {$ ]$ M5 cS(t)+I(t)恒等于N,这里,N是个体总数。随着时间t的增长,易感个体与感染个体的接触
    0 G6 Y& z4 f, B  B会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β,疾病仅在
    4 A: O5 q' l8 Q, U& m+ E8 K" m感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t,易感个体的比例为S(t)/N,感染个体的数量为I(t),一次,易感个体的数量将以如下变化率减少$ t- Q- l6 @9 Y+ G2 y/ E
    ds/dt = -β*S(t)I(t)/N
    , |, A( J# B9 Q  L2 s/ ?0 T3 B% ^同时,感染个体的数量会以与易感个体相反的变化率增加,
    0 q6 x" i4 z" i- Sds/dt = βS(t)*I(t)/N6 C9 i% _/ S4 h( e3 d5 f/ `
    分别将时刻t处于易感状态和感染状态的个体所占比例记为,
    4 w/ {( R( I8 p% |% Z9 rs(t)=S(t)/N
    6 A, h% Z, _9 b$ P4 J. Z: [i(t)=I(t)/N" P" J$ P+ X2 w$ f; `
    显然有," D; Q3 O' ]' t6 d
    s(t)+i(t)恒等于1,此时之前的公式可以记做
    * u0 x0 I5 C  b1 ~: i) I: Z& qds/dt=-βsi' U! W7 Z! D- J( r9 ]6 \
    di/dt=βsi
    ; y3 B& z. U  q8 O) u
    $ Q7 t) W8 ?5 k( H( n! Z) U, Xdi/dt=βi(1-i)
    7 B; N: U* |" m4 _$ Y0 ?0 _上式也成为Logistic增长方程式(Logistic growth equation),
    5 `9 j& q! U  _, T& X- d方程的解和图像如图) ~, x: p' k" s+ X
    1.jpg 9 K& a# w* j. `
    代码和相关文件以及环境链接:链接:https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ9 ?' i! O* D, c1 C$ o
    提取码:z4484 ~& p* \4 L& Z
    : M  g# w$ r: H. P  e) L; ]) D

    1 K, I) E; t& r; ]1 E'''
    + L' g! y/ a/ n$ ~, T. {实验环境Python2.7.13,igraph包,cairo包,numpy包1 z, Z! [, v- \5 }2 x8 }# E
    '''
    * B- k" ]) V# {  [0 L# -*- coding:utf8 -*- m; ?; T4 @0 L; C  b1 r
    from igraph import *# O5 {5 c: d' Z) {; M. J3 c: t9 j9 L
    import numpy as numpy
    / I, I9 y% [: V  t: Dfrom  numpy import *. s; r6 O, @- j- ]: J9 n9 M
    import random
    . q7 A9 q( K7 a% v4 U+ d. R' V& i4 y) v6 k- {8 M' g! i9 c& L0 K: t) [
    def len_arr(infected_array,nodes_num):#获取感染数组长度) `6 @. L: U- m! b/ c# ~+ ]  [9 z
        len_value=0#初始化长度
    & r3 F* {  e3 I4 M# H' E4 Y    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)
    8 h; t/ P7 K) V  m' U1 \# |, r( a9 j    return len_value- O" m- m& d; e9 b6 M
    ( K& Z" O0 C# j+ C
    g=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)/ v' P3 M+ d0 E) m3 H- c
    summary(g)
    * C9 \+ b9 U! a4 c) \, g/ |nodes_num=g.vcount()#统计图中的结点个数
    9 Z) ^* q7 K  y2 n  a) P3 t5 ~net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
    ( d* W4 Z. N, {5 [) R% K& fg.vs["color"]=["white"]#给图的顶点序列颜色赋值白色7 w1 r: [) \8 k
    a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a  {' A% x+ G# O/ d0 O4 W
    nodes_state=matrix(a).T#nodes_state通过转置a矩阵创建,用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重!!!), z4 e* N$ e( S4 p5 G4 c  K
                                                        #第一列是节点编号,第二列是节点状态,感染状态用-2表示,第三列是节点感染的时间! D# o0 @. o9 s$ @8 ~" V
    print(nodes_state); U1 \4 {, x6 I. v0 B) v
    infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数3 V' ]5 p, x; h5 P  H
    print(infected_array)
    + p$ E' O# s7 b. E& E8 g/ D6 ^) ^9 E5 H
    infe_rate=1#传播率(感染率) 1代表邻接点100%被感染+ b  q1 _) z, O  ?$ U( I
    set_time=2#传播次数(感染次数) 2次
    8 b/ b8 j6 z- t" ^9 Xsource_seed=1#感染源位置0 L( Z2 x- T/ L
    nodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1
    / g* B6 r" }( F5 f2 I: V" h! ]3 Onodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态
    ; H8 c0 |  ]: Z. j1 T% Z( @nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为1
    # |, ^) `, _7 |6 ^9 K  n' {+ }g.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红
    , t- i6 i  i: I, r: Q! q! Kinfected_array[0]=source_seed#将感染源的位置存入被感染节点列表$ A# D8 O# d+ h& g7 a3 P" }
    plot(g)#绘制
    ( ^, w% U" E$ F: J% @. [6 H, y4 R
    0 J! O) O5 a# h' O  B4 G& \stop=False#感染过程结束的标记
    " g9 B! A  W# N5 K- gtemp_time=0#第几次感染  u  L; q* B" e4 m, w' R
    temp_len=0#本轮的感染源数量初始化. v! Q; A8 V/ G& ?2 y9 m
    ! M5 \/ h3 B% m6 P& r  L: F2 f  p
    while not stop:0 m& R/ ^- f8 y/ G9 t
        i=0#记录让每个感染源都传播一次! P; O( o. i! E) I7 B
        if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行
      ^6 s5 K8 z" T( t% a. [' o/ J        temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量
      C0 f. \2 o) ]3 B, p3 e! V        while i<temp_len:9 W: \" Z, B7 I
                temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间7 [. J1 F, R/ P2 ]
                nei_count=0#下一轮可以被感染到的节点数量
    ( d' I9 |- ?; B2 W/ l) X4 U            #生成下一轮可能被感染的节点的集合nei_arr0 O0 ]2 j9 ~# R/ X9 h7 n
                for j in range(nodes_num):#遍历节点9 O. f4 K" K7 d& h
                    if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染
    2 }( o4 B) n  s' m* ?* e                    nei_count=nei_count+1#下一轮可以被感染到的节点数量++4 }3 g0 N( M- a+ S0 B5 s
                nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
    / z% Y! V/ k! a$ ^6 U' E            t=0
    6 B5 d. @/ W- y- U$ I0 o1 i9 X: V            for j in range(nodes_num):6 A% r+ w6 R* `, _3 v+ M; Y
                    if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:; j) l# X( G/ M' D
                        nei_arr[t]=j+1" F5 ?, d+ A/ j/ B
                        t=t+1
    % @9 g1 y/ L( Q            ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组8 ^3 k% W+ V5 Y+ {( E6 q5 W5 j0 o
                                            #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象. O1 X' ?+ A) M5 n6 E( ?. [, \9 W
                if len(ran_infe_arr)>0:#存在需要被感染的节点
    5 k, j9 R6 s+ B5 I- l! y* U+ }                t=0#让ran_infe_arr内每个感染源都被感染
    9 T  r4 D9 X+ D, o# v8 R                while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染# a# W& R+ c) j9 B8 B& w& c/ k
                        nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
    4 J6 a0 x$ A( J" l                    nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间
    9 n$ [/ M- ?  e8 l  i2 G                    infected_array[len_arr(infected_array,nodes_num)]=nei_arr[ran_infe_arr[t]]#将此次感染节点放入总的感染节点数组中
    0 `5 l: o# k9 G$ x* A. S                    g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色
    * E; a* b) F+ j1 q                    plot(g)#绘制
    . B) v0 @1 B# o& U                    t=t+1
    ; n- n7 \& o! K1 Z            i=i+1
    . w/ f6 l+ h- E; ^7 v8 s    if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染
    ! v$ h. M3 Y: a) e$ Y8 f        stop=True / W% P$ ^7 m3 r! D' G( V7 h1 O6 E
    & r0 j6 m2 L' k1 S1 S
    ) n' R7 G! r6 H$ ~/ H% K6 O
    视频演示bilibili传送门
    1 L  _' Q# p) b( j( Z; F效果图3 V% N* c  ?- G3 R7 {+ ?5 V- c

    / g: T. t. y! p
    " Q6 s3 c2 E/ N8 z% M 2.jpg # c$ ]1 V7 s* g( \' n' L( t0 ^

    4 I; }; V3 ]& z6 [: S 3.png 5 b; F* j" W1 c, l& e6 K0 ^
    . U* `1 z$ {( F7 `! K* ]/ W
    4.jpg 1 G) Y5 h) _* }2 P# M* i
    6 W$ D+ e" A- a+ _+ @; v
    5.png   x8 `3 L2 i( Y: @7 D5 r
    0 u8 J6 X+ {$ i+ m) x9 A8 p- u, ]
    6.png
    " ?2 n, m' i4 o: J: b+ f( G3 e$ A+ ?, }& s* m" R
    7.png ————————————————
    % M8 {. ?$ e. J% X6 k7 P版权声明:本文为CSDN博主「eck_燃」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。* k/ ]' m  a' A' T, ]1 \
    原文链接:https://blog.csdn.net/wdays83892469/article/details/808788620 ]" _; B( V' ^& c/ I
    5 ?% C+ r. _6 A) M

    3 n8 b2 {0 k3 i  [7 B
    zan
    转播转播0 分享淘帖0 分享分享0 收藏收藏0 支持支持0 反对反对0 微信微信
    尔雅 实名认证       

    0

    主题

    1

    听众

    208

    积分

    升级  54%

  • TA的每日心情
    开心
    2022-2-18 09:12
  • 签到天数: 17 天

    [LV.4]偶尔看看III

    网络挑战赛参赛者

    国际赛参赛者

    回复

    使用道具 举报

    您需要登录后才可以回帖 登录 | 注册地址

    qq
    收缩
    • 电话咨询

    • 04714969085
    fastpost

    关于我们| 联系我们| 诚征英才| 对外合作| 产品服务| QQ

    手机版|Archiver| |繁體中文 手机客户端  

    蒙公网安备 15010502000194号

    Powered by Discuz! X2.5   © 2001-2013 数学建模网-数学中国 ( 蒙ICP备14002410号-3 蒙BBS备-0002号 )     论坛法律顾问:王兆丰

    GMT+8, 2026-4-21 04:18 , Processed in 0.456979 second(s), 59 queries .

    回顶部