拓端tecdat|R語言SIR模型（Susceptible Infected Recovered Model）代碼sir模型實例

原文鏈接：http://tecdat.cn/?p=14593 

 

SIR模型定義

SIR模型是一種傳播模型，是信息傳播過程的抽象描述。
SIR模型是傳染病模型中最經典的模型，其中S表示易感者，I表示感染者，R表示移除者。

S：Susceptible，易感者
I：Infective，感染者
R：Removal，移除者

---

SIR模型的應用

SIR模型應用於信息傳播的研究。

傳播過程大致如下：最初，所有的節點都處於易感染狀態。然后，部分節點接觸到信息后，變成感染狀態，這些感染狀態的節點試着去感染其他易感染狀態的節點，或者進入恢復狀態。感染一個節點即傳遞信息或者對某事的態度。恢復狀態，即免疫，處於恢復狀態的節點不再參與信息的傳播。

SIR的微分方程

a為感染率、b恢復率

注意：

t為某個時刻，例如t=1，S(1)為第一天易感人群的人數。
無論t為什么時刻，總人數是不變的，即N(t)=S(t)+I(t)+R(t)。
人口總數總保持一個常數，即N(t)=K，不考慮人口的出生、死亡、遷移等因素。

這里介紹一個使用R模擬網絡擴散的例子。

 

第一步，生成網絡。

規則網

1.  
   g =graph.tree(size, children =2); plot(g)
2.  
    

1.  
   g =graph.star(size); plot(g)
2.  
    

1.  
   g =graph.full(size); plot(g)
2.  
    

1.  
   g =graph.ring(size); plot(g)
2.  
    

第二步，隨機選取一個或n個隨機種子。
 

1.  
   # initiate the diffusers
2.  
    
3.  
   seeds_num =1 diffusers =sample(V(g),seeds_num) ;
4.  
    
5.  
   diffusers
6.  
    
7.  
   ## + 1/50 vertex:
8.  
    
9.  
   ## [1] 43
10.  
     
11.  
    infected =list()
12.  
     
13.  
    infected[[1]]=diffusers#

 

第三步，傳染能力

在這個簡單的例子中，每個節點的傳染能力是0.5，即與其相連的節點以0.5的概率被其感染,每個節點的回復能力是0.5，即其以0.5的概率被其回復。在R中的實現是通過拋硬幣的方式來實現的。

1.  
   ## [1] 0
2.  
    

顯然，這很容易擴展到更一般的情況，比如節點的平均感染能力是0.128，那么可以這么寫： 節點的平均回復能力是0.1，那么可以這么寫

1.  
   p =0.128
2.  
    
3.  
   coins =c(rep(1, p*1000), rep(0,(1-p)*1000))
4.  
    
5.  
   sample(coins, 1, replace=TRUE, prob=rep(1/n, n))
6.  
    
7.  
   ## [1] 0
8.  
    
9.  
   n =length(coins2)
10.  
     
11.  
    sample(coins2, 1, replace=TRUE, prob=rep(1/n, n))
12.  
     
13.  
    ## [1] 0

當然最重要的一步是要能按照“時間”更新網絡節點被感染的信息。
 

1.  
   keep =unlist(lapply(nearest_neighbors[,2], toss))
2.  
    
3.  
   new_infected =as.numeric(as.character(nearest_neighbors[,1][keep >=1]))
4.  
    
5.  
   diffusers =unique(c(as.numeric(diffusers), new_infected))
6.  
    
7.  
   return(diffusers)}
8.  
    
9.  
   set.seed(1);

 

開啟擴散過程！

先看看S曲線吧：

為了可視化這個擴散的過程，我們用紅色來標記被感染者。

1.  
   # generate a palette#
2.  
    
3.  
   plot(g, layout =layout.old)
4.  
    
5.  
   set.seed(1)#
6.  
    
7.  
   library(animation)# start the plot
8.  
    
9.  
   m =1

如同在Netlogo里一樣，我們可以把網絡擴散與增長曲線同時展示出來：

1.  
   set.seed(1)
2.  
    
3.  
   # start the plot
4.  
    
5.  
   m =1
6.  
    
7.  
   p_cum=numeric(0)
8.  
    
9.  
   h_cum=numeric(0)
10.  
     
11.  
    i_cum=numeric(0)
12.  
     
13.  
    while( m<50 ) {# start the plot
14.  
     
15.  
    layout(matrix(c(1, 2, 1, 3), 2,2, byrow =TRUE), widths=c(3,1), heights=c(1, 1))
16.  
     
17.  
    V(g)$color = "white"
18.  
     
19.  
    V(g)$color[V(g)%in%infected[[m ]] ] = "red"
20.  
     
21.  
    V(g)$color[V(g)%in%health[[m ]]] = "green"
22.  
     
23.  
    if(m<=length(infected))
24.  
     
25.  
     
26.  
     
27.  
     
28.  
     
29.  
     
30.  
     
31.  
     
32.  
     
33.  
     
34.  
     
35.  
    plot(pp~time, type ="h", ylab ="PDF", xlab ="Time",xlim =c(0,i), ylim =c(0,1), frame.plot =FALSE)
36.  
     
37.  
    m =m +1
38.  
     
39.  
    }

 

---

參考文獻

1.R語言泊松Poisson回歸模型分析案例

2.R語言進行數值模擬：模擬泊松回歸模型

3.r語言泊松回歸分析

4.R語言對布豐投針（蒲豐投針）實驗進行模擬和動態可視化

5.用R語言模擬混合制排隊隨機服務排隊系統

6.GARCH（1,1），MA以及歷史模擬法的VaR比較

7.R語言做復雜金融產品的幾何布朗運動的模擬

8.R語言進行數值模擬：模擬泊松回歸模型

9.R語言對巨災風險下的再保險合同定價研究案例：廣義線性模型和帕累托分布Pareto distributions