QQ登录

只需要一步,快速开始

 注册地址  找回密码
查看: 3073|回复: 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传播模型' O+ b3 [9 P% K, X, D. G
    #SI疾病传播模型的原理
    8 p1 y7 T1 y. j3 }$ `5 E在经典的传染病模型中,种群(Population)内N个个体的状态可分为如下几类3 l* z( q+ G6 w- S
    4 S5 A; d& O- u( d: h
    易感状态(Susceptible)。一个个体在感染前是处于易感状态的,即该个体有可能被邻居个体感染。
      l. L$ _- K+ ?/ N- v+ |易感状态I(Infected)。一个感染上某种病毒的个体就称为是处于感染状态。,即该个体还会以一定概率感染其邻居个体。
    4 L/ ]* r, K% N1 T3 _4 P" \- L+ S移除状态(Remove,Refractory或者Recovered)。也成为免疫状态或恢复状态,当一个个体经历过一个完整的感染周期后,该个体就不再被感染,因此就可以不再考虑改革提。+ Z2 @& R3 K2 U! f/ B. O
    SI传播模型是最简单的疾病传播模型,模型中的所有个体都只可能处于两个状态中的一个  V0 g: E( @4 n, m# G0 {0 V
    即易感(S)状态或感染(I)状态。SI模型中的个体一旦被感染后就永远处于感染状态。: ^# P( p/ B, e* x5 [8 U
    在给定时刻t,令S(t)与I(t)分别代表该时刻处于易感和感染状态的个体数目,显然有0 w1 r. I& G' i7 n
    S(t)+I(t)恒等于N,这里,N是个体总数。随着时间t的增长,易感个体与感染个体的接触
    , h% p5 s! C5 y) |! k会导致感染个体数量的增加。加入由于个体之间的接触而导致疾病传播的概率为β,疾病仅在& {! U9 R6 u7 K6 f* I* J7 X
    感染个体和易感个体之间进行接触时才会以概率β将疾病传染给易感个体。在时刻t,易感个体的比例为S(t)/N,感染个体的数量为I(t),一次,易感个体的数量将以如下变化率减少0 y; h0 b# |% Q) j$ m# N/ J
    ds/dt = -β*S(t)I(t)/N; R8 k0 @1 e7 e8 J5 K; x
    同时,感染个体的数量会以与易感个体相反的变化率增加,
    9 g3 n; W1 n. J  eds/dt = βS(t)*I(t)/N
    3 d- x6 p8 |) c* M" v# m4 ~分别将时刻t处于易感状态和感染状态的个体所占比例记为,# J" F7 \  H1 A& D8 E
    s(t)=S(t)/N% W' l9 ?$ Q; z
    i(t)=I(t)/N
    3 b* H& `4 ]+ I1 x( ~' }7 \显然有,+ {7 ~# R! w1 S4 E  }/ s; h
    s(t)+i(t)恒等于1,此时之前的公式可以记做
    ( O9 C, F& ^% F) H* F) N0 l2 ~- lds/dt=-βsi
    % Q1 O3 M/ o- A( o4 r  idi/dt=βsi7 T" E* K( \' X1 w. @  _* b

    1 L% o8 H' c' K9 @4 Hdi/dt=βi(1-i)0 W8 {6 x8 L) L7 ]$ R  b
    上式也成为Logistic增长方程式(Logistic growth equation),. J' x0 ~8 d. ~- w' j
    方程的解和图像如图
    # R/ K: O3 q) U- x8 Q" w 1.jpg $ y0 @, R$ s8 m! @
    代码和相关文件以及环境链接:链接:https://pan.baidu.com/s/1JSfHuTPaglFimeEBLdSDyQ9 P, o! i. x1 g& Z9 S( _
    提取码:z448
    ) W& T: f* }& P& ^
    ( q9 {, ]% o5 c0 N5 h0 x. b3 T. z" c3 G: {* F" _" `- `6 I
    '''
    0 F( \6 e, C9 ]8 b9 y实验环境Python2.7.13,igraph包,cairo包,numpy包7 f) U% f1 D3 u, P/ d0 o! Q1 Y3 {' O
    '''
    ) U) s' g$ w8 V1 I( ~, P9 M# -*- coding:utf8 -*( X2 e: W- ^6 W/ b6 T; M3 w& J$ \
    from igraph import *9 ^5 w* X- x  }4 h* N% @' N
    import numpy as numpy& D3 L! i/ S/ T! m
    from  numpy import *
    & D& @3 c4 ]4 l0 m% `1 yimport random
      j( H, z# D* S: }$ H" O6 Z* A' p! H: S9 Q3 Z8 q
    def len_arr(infected_array,nodes_num):#获取感染数组长度
    % n# L% m2 h, K3 r    len_value=0#初始化长度
    1 H* o. {6 w( A3 Q    len_value=nodes_num-infected_array.count(-1)#被感染数量是结点总数减去未感染节点数(未感染的结点被标记为-1)1 ?7 N; _% H+ I
        return len_value% G5 U5 H" H3 z) \

      \* x0 _/ H( c' Q4 pg=Graph.Read_GML("C:\python27\e1.gml")#将本地保存的网络数据读入变量g(生成图)
      u9 z6 Z$ `  s0 T8 i% D  ^/ Fsummary(g)5 H  \- E7 y+ ]! y: l$ P5 H6 y% N
    nodes_num=g.vcount()#统计图中的结点个数
    & p7 a% k, L4 M4 z3 ^net_mat=g.get_adjacency(type=GET_ADJACENCY_BOTH)#将网络数据转换为邻接矩阵存储在变量net_mat
    * m/ [' n6 g; L) @( R0 o, m! N' |g.vs["color"]=["white"]#给图的顶点序列颜色赋值白色6 B+ O4 Z. F. k! S( |
    a=[arange(nodes_num)+1]*3#声明一个N行3列的数组a. h1 S4 y$ V. U, S* H; u# E3 j
    nodes_state=matrix(a).T#nodes_state通过转置a矩阵创建,用于存放每个节点的状态信息以及其被感染的时间(这个是理解算法的重中之重!!!)
    9 s: D4 X" P9 R                                                    #第一列是节点编号,第二列是节点状态,感染状态用-2表示,第三列是节点感染的时间" a( O- g" l5 P6 d; T& w) s) m1 H& @
    print(nodes_state)& O0 Y; j. L/ ^2 a& _
    infected_array=[-1]*34#用于存放本轮被感染的结点, 这些结点将参与下一次感染   34代表网络节点数  M4 z. Y& t( n2 j9 k! S3 N+ O
    print(infected_array)
    / ^7 E, s: k% B. i% a$ d# y$ C  \- e5 r6 r& d9 _" @
    infe_rate=1#传播率(感染率) 1代表邻接点100%被感染
    7 U+ r  b6 M/ H: M5 Oset_time=2#传播次数(感染次数) 2次1 D3 t5 a+ C. \4 O
    source_seed=1#感染源位置, T+ t& o9 K( S$ n, t
    nodes_state[0:nodes_num,2]=-1#给所有节点初始化感染时间为-1' h6 _& B# k! g6 c" T
    nodes_state[source_seed-1,1]=-2#设置第一个感染源感染状态 -2代表感染状态! b5 G4 L+ f4 ^5 h8 H
    nodes_state[source_seed-1,2]=1#设置第一个感染源的感染时间为11 J6 y1 q- p* P6 I% M  b* _. c
    g.vs[source_seed-1]["color"]="red"#将感染的顶点颜色标红+ g8 w% G- n7 a/ j6 J6 L
    infected_array[0]=source_seed#将感染源的位置存入被感染节点列表
      d- r- j9 ~" a# `$ w6 Jplot(g)#绘制
    % B( F6 Z) {, E
    9 n! m- v1 x" l) W- b5 _6 tstop=False#感染过程结束的标记$ B) B( h0 j5 h5 V3 c$ s
    temp_time=0#第几次感染$ r2 U' Y( c4 i1 X
    temp_len=0#本轮的感染源数量初始化8 c* A/ L3 A% J" F5 q+ I

    3 F+ B9 f% j1 S& o" G4 ^& I  |# K! Qwhile not stop:
    & w) q5 Y) U+ v/ `0 m    i=0#记录让每个感染源都传播一次
    9 v+ f9 ^+ d0 W' `    if len_arr(infected_array,nodes_num)>0 and len_arr(infected_array,nodes_num)<=nodes_num:#感染可以进行3 Q! E# A3 W' B+ q  l# Y" |
            temp_len=len_arr(infected_array,nodes_num)#获取本轮的感染源数量8 b4 h  c7 b4 \3 b
            while i<temp_len:
    . v7 ?  A, V' I+ I( K# p0 C- F            temp_time=nodes_state[infected_array-1,2]#获取每一个节点的感染时间
    6 ^6 S0 r; o( I' ^            nei_count=0#下一轮可以被感染到的节点数量& P* s6 ]# ?7 T- V
                #生成下一轮可能被感染的节点的集合nei_arr  p, v5 b& G. T
                for j in range(nodes_num):#遍历节点4 w1 ?5 Z) @) V2 f9 O+ R
                    if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:#是邻接节点而且未被感染* n# u* z* ~9 A1 Y- F# d2 u
                        nei_count=nei_count+1#下一轮可以被感染到的节点数量++: m$ I, E' V& [
                nei_arr=[-1]*nei_count#用于临时存放本轮被感染的结点, 这些结点将参与下一次感染
    + I7 }5 h% j6 B3 {8 E            t=0; M9 T$ _: ~+ w- a
                for j in range(nodes_num):( s- E. \5 A8 M! |! w
                    if net_mat[infected_array-1,j]==1 and nodes_state[j,1]!=-2:
    ' Z- f. Y3 L6 f2 ^1 P% s                    nei_arr[t]=j+1
    3 V/ g2 @& [' k5 k5 h                    t=t+17 Z( W' E: @6 Z& N  [2 U
                ran_infe_arr=random.sample(range(nei_count),int(nei_count*infe_rate))#随机生成会被感染的节点的数组
    - E7 f$ u; T& a                                        #random.simple(arg1,num) 从arg1集合中随机取num个数据生成一个对象& Q* k% n- x( D& ]& Y
                if len(ran_infe_arr)>0:#存在需要被感染的节点% `- v2 j+ l3 E+ x5 b* A
                    t=0#让ran_infe_arr内每个感染源都被感染' N$ B$ x9 c# \0 f9 f
                    while t<len(ran_infe_arr):#对刚才生成的会被感染的数组内的节点进行感染
    ' ~. G9 s( i' N) v9 L1 K4 N                    nodes_state[nei_arr[ran_infe_arr[t]]-1,1]=-2#标记为感染状态
    ) h6 D' B' v) @. R% ^- D; ?                    nodes_state[nei_arr[ran_infe_arr[t]]-1,2]=temp_time+1#记录感染时间7 L/ L# ?3 t. G
                        infected_array[len_arr(infected_array,nodes_num)]=nei_arr[ran_infe_arr[t]]#将此次感染节点放入总的感染节点数组中
    7 m% k& A% U2 m                    g.vs[nei_arr[ran_infe_arr[t]]-1]["color"]="pink"#将此次感染的节点集的所有节点颜色置为粉色6 E7 s8 C* V4 p8 e
                        plot(g)#绘制
    1 |) T; B) |0 r% w                    t=t+1/ j' z5 {5 y( V  k: [3 t4 {
                i=i+1* m6 K& v- r$ p
        if temp_time>set_time-1:#当执行感染的次数等于设置的次数结束感染# w0 R# f( O# _0 C. M4 ~) c$ j
            stop=True
    9 u# j2 y2 P# U
    & G9 v$ [6 e2 D' z5 K% m$ D/ E1 }7 W% B  v$ n( u' H! R. J
    视频演示bilibili传送门
    8 i+ v! y- w2 a效果图
    & U. R0 e) Y; ]8 c+ l9 D1 V+ W8 b( k8 p8 P1 Q
    $ @4 H- ^1 P% F, \, G
    2.jpg 6 D# N+ S- @9 {6 q1 i; ^; }
    8 W- b5 q1 n, c  j$ ~& j
    3.png
    0 G7 y8 d- m8 x; V/ y0 a
    7 w$ t: F. s" J6 u: T( w 4.jpg & {( l* P# [. C! p/ e' `% d2 f! f
    & ~% t9 f" H4 x- ^
    5.png
    / }: G$ n9 F6 \2 g) o! T& w. e2 r+ B; v- Z: r: k+ W5 i1 {" F
    6.png
    , T# l- B) e0 O9 O  F3 q
    7 W' \9 g9 |0 c/ z$ Y- o 7.png ————————————————
    % m/ k. u% {6 a, s- A4 ~7 X7 ]1 j- y版权声明:本文为CSDN博主「eck_燃」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。9 _3 t$ a. X$ r
    原文链接:https://blog.csdn.net/wdays83892469/article/details/80878862
    " q- g4 U; N! c6 j& S5 u3 T. o8 d' h0 e$ Y! I, _/ E

    - t+ b, u  i+ @3 g8 S# h
    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 05:07 , Processed in 0.438488 second(s), 59 queries .

    回顶部